diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index 59c6cc69d..1e47215ae 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -1,518 +1,531 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Pieter in 't Veld (SNL) ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_colloid.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neigh_list.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) enum{SMALL_SMALL,SMALL_LARGE,LARGE_LARGE}; /* ---------------------------------------------------------------------- */ PairColloid::PairColloid(LAMMPS *lmp) : Pair(lmp) {} /* ---------------------------------------------------------------------- */ PairColloid::~PairColloid() { if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); memory->destroy_2d_int_array(form); memory->destroy_2d_double_array(a12); memory->destroy_2d_double_array(sigma); memory->destroy_2d_double_array(d1); memory->destroy_2d_double_array(d2); memory->destroy_2d_double_array(a1); memory->destroy_2d_double_array(a2); + memory->destroy_2d_double_array(diameter); memory->destroy_2d_double_array(cut); memory->destroy_2d_double_array(offset); memory->destroy_2d_double_array(sigma3); memory->destroy_2d_double_array(sigma6); memory->destroy_2d_double_array(lj1); memory->destroy_2d_double_array(lj2); memory->destroy_2d_double_array(lj3); memory->destroy_2d_double_array(lj4); } } /* ---------------------------------------------------------------------- */ void PairColloid::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r,forcelj,factor_lj; double r2inv,r6inv,c1,c2,fR,dUR,dUA; double K[9],h[4],g[4]; int *ilist,*jlist,*numneigh,**firstneigh; evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq >= cutsq[itype][jtype]) continue; switch (form[itype][jtype]) { case SMALL_SMALL: r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; if (eflag) evdwl = r6inv*(r6inv*lj3[itype][jtype]-lj4[itype][jtype]) - offset[itype][jtype]; break; case SMALL_LARGE: c2 = a2[itype][jtype]; K[1] = c2*c2; K[2] = rsq; K[0] = K[1] - rsq; K[4] = rsq*rsq; K[3] = K[1] - K[2]; K[3] *= K[3]*K[3]; K[6] = K[3]*K[3]; fR = sigma3[itype][jtype]*a12[itype][jtype]*c2*K[1]/K[3]; fpair = 4.0/15.0*sqrt(rsq)*fR*factor_lj * (2.0*(K[1]+K[2]) * (K[1]*(5.0*K[1]+22.0*K[2])+5.0*K[4]) * sigma6[itype][jtype]/K[6]-5.0) / K[0]; if (eflag) evdwl = 2.0/9.0*fR * (1.0-(K[1]*(K[1]*(K[1]/3.0+3.0*K[2])+4.2*K[4])+K[2]*K[4]) * sigma6[itype][jtype]/K[6]) - offset[itype][jtype]; break; case LARGE_LARGE: r = sqrt(rsq); c1 = a1[itype][jtype]; c2 = a2[itype][jtype]; K[0] = c1*c2; K[1] = c1+c2; K[2] = c1-c2; K[3] = K[1]+r; K[4] = K[1]-r; K[5] = K[2]+r; K[6] = K[2]-r; K[7] = 1.0/(K[3]*K[4]); K[8] = 1.0/(K[5]*K[6]); g[0] = pow(K[3],-7.0); g[1] = pow(K[4],-7.0); g[2] = pow(K[5],-7.0); g[3] = pow(K[6],-7.0); h[0] = ((K[3]+5.0*K[1])*K[3]+30.0*K[0])*g[0]; h[1] = ((K[4]+5.0*K[1])*K[4]+30.0*K[0])*g[1]; h[2] = ((K[5]+5.0*K[2])*K[5]-30.0*K[0])*g[2]; h[3] = ((K[6]+5.0*K[2])*K[6]-30.0*K[0])*g[3]; g[0] *= 42.0*K[0]/K[3]+6.0*K[1]+K[3]; g[1] *= 42.0*K[0]/K[4]+6.0*K[1]+K[4]; g[2] *= -42.0*K[0]/K[5]+6.0*K[2]+K[5]; g[3] *= -42.0*K[0]/K[6]+6.0*K[2]+K[6]; fR = a12[itype][jtype]*sigma6[itype][jtype]/r/37800.0; evdwl = fR * (h[0]-h[1]-h[2]+h[3]); dUR = evdwl/r + 5.0*fR*(g[0]+g[1]-g[2]-g[3]); dUA = -a12[itype][jtype]/3.0*r*((2.0*K[0]*K[7]+1.0)*K[7] + (2.0*K[0]*K[8]-1.0)*K[8]); fpair = factor_lj * (dUR+dUA)/r; if (eflag) evdwl += a12[itype][jtype]/6.0 * (2.0*K[0]*(K[7]+K[8])-log(K[8]/K[7])) - offset[itype][jtype]; break; } if (eflag) evdwl *= factor_lj; f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,0.0,fpair,delx,dely,delz); } } if (vflag_fdotr) virial_compute(); } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairColloid::allocate() { allocated = 1; int n = atom->ntypes; setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); form = memory->create_2d_int_array(n+1,n+1,"pair:form"); a12 = memory->create_2d_double_array(n+1,n+1,"pair:a12"); sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); d1 = memory->create_2d_double_array(n+1,n+1,"pair:d1"); d2 = memory->create_2d_double_array(n+1,n+1,"pair:d2"); a1 = memory->create_2d_double_array(n+1,n+1,"pair:a1"); a2 = memory->create_2d_double_array(n+1,n+1,"pair:a2"); + diameter = memory->create_2d_double_array(n+1,n+1,"pair:diameter"); cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); offset = memory->create_2d_double_array(n+1,n+1,"pair:offset"); sigma3 = memory->create_2d_double_array(n+1,n+1,"pair:sigma3"); sigma6 = memory->create_2d_double_array(n+1,n+1,"pair:sigma6"); lj1 = memory->create_2d_double_array(n+1,n+1,"pair:lj1"); lj2 = memory->create_2d_double_array(n+1,n+1,"pair:lj2"); lj3 = memory->create_2d_double_array(n+1,n+1,"pair:lj3"); lj4 = memory->create_2d_double_array(n+1,n+1,"pair:lj4"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairColloid::settings(int narg, char **arg) { if (narg != 1) error->all("Illegal pair_style command"); cut_global = atof(arg[0]); // reset cutoffs that have been explicitly set if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; } } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairColloid::coeff(int narg, char **arg) { if (narg < 6 || narg > 7) error->all("Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); double a12_one = atof(arg[2]); double sigma_one = atof(arg[3]); double d1_one = atof(arg[4]); double d2_one = atof(arg[5]); double cut_one = cut_global; if (narg == 7) cut_one = atof(arg[6]); if (d1_one < 0.0 || d2_one < 0.0) error->all("Invalid d1 or d2 value for pair colloid coeff"); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { a12[i][j] = a12_one; sigma[i][j] = sigma_one; if (i == j && d1_one != d2_one) error->all("Invalid d1 or d2 value for pair colloid coeff"); d1[i][j] = d1_one; d2[i][j] = d2_one; + diameter[i][j] = 0.5*(d1_one+d2_one); cut[i][j] = cut_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all("Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairColloid::init_one(int i, int j) { if (setflag[i][j] == 0) { a12[i][j] = mix_energy(a12[i][i],a12[j][j],sigma[i][i],sigma[j][j]); sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); d1[i][j] = mix_distance(d1[i][i],d1[j][j]); d2[i][j] = mix_distance(d2[i][i],d2[j][j]); + diameter[i][j] = 0.5 * (d1[i][j] + d2[i][j]); cut[i][j] = mix_distance(cut[i][i],cut[j][j]); } sigma3[i][j] = sigma[i][j]*sigma[i][j]*sigma[i][j]; sigma6[i][j] = sigma3[i][j]*sigma3[i][j]; if (d1[i][j] == 0.0 && d2[i][j] == 0.0) form[i][j] = SMALL_SMALL; else if (d1[i][j] == 0.0 || d2[i][j] == 0.0) form[i][j] = SMALL_LARGE; else form[i][j] = LARGE_LARGE; // for SMALL_SMALL, a1/a2 do not need to be set // for SMALL_LARGE, a1 does not need to be set, a2 = radius for i,j and j,i // for LARGE_LARGE, a1/a2 are radii, swap them for j,i if (form[i][j] == SMALL_LARGE) { if (d1[i][j] > 0.0) a2[i][j] = 0.5*d1[i][j]; else a2[i][j] = 0.5*d2[i][j]; a2[j][i] = a2[i][j]; } else if (form[i][j] == LARGE_LARGE) { a2[j][i] = a1[i][j] = 0.5*d1[i][j]; a1[j][i] = a2[i][j] = 0.5*d2[i][j]; } form[j][i] = form[i][j]; a12[j][i] = a12[i][j]; sigma[j][i] = sigma[i][j]; sigma3[j][i] = sigma3[i][j]; sigma6[j][i] = sigma6[i][j]; + diameter[j][i] = diameter[i][j]; cut[j][i] = cut[i][j]; cutsq[j][i] = cutsq[i][j] = cut[i][j] * cut[i][j]; double epsilon = a12[i][j]/144.0; lj1[j][i] = lj1[i][j] = 48.0 * epsilon * sigma6[i][j] * sigma6[i][j]; lj2[j][i] = lj2[i][j] = 24.0 * epsilon * sigma6[i][j]; lj3[j][i] = lj3[i][j] = 4.0 * epsilon * sigma6[i][j] * sigma6[i][j]; lj4[j][i] = lj4[i][j] = 4.0 * epsilon * sigma6[i][j]; offset[j][i] = offset[i][j] = 0.0; if (offset_flag) { double tmp; offset[j][i] = offset[i][j] = single(0,0,i,j,cutsq[i][j],0.0,1.0,tmp); } return cut[i][j]; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairColloid::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&a12[i][j],sizeof(double),1,fp); fwrite(&sigma[i][j],sizeof(double),1,fp); fwrite(&d1[i][j],sizeof(double),1,fp); fwrite(&d2[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairColloid::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (comm->me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (comm->me == 0) { fread(&a12[i][j],sizeof(double),1,fp); fread(&sigma[i][j],sizeof(double),1,fp); fread(&d1[i][j],sizeof(double),1,fp); fread(&d2[i][j],sizeof(double),1,fp); fread(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&a12[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&d1[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&d2[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairColloid::write_restart_settings(FILE *fp) { fwrite(&cut_global,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairColloid::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { fread(&cut_global,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- */ double PairColloid::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double K[9],h[4],g[4]; double r,r2inv,r6inv,forcelj,c1,c2,phi,fR,dUR,dUA; switch (form[itype][jtype]) { case SMALL_SMALL: r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fforce = factor_lj*forcelj*r2inv; phi = r6inv*(r6inv*lj3[itype][jtype]-lj4[itype][jtype]) - offset[itype][jtype]; break; case SMALL_LARGE: r = sqrt(rsq); c2 = a2[itype][jtype]; K[1] = c2*c2; K[2] = rsq; K[4] = rsq*rsq; K[3] = K[1] - K[2]; K[3] *= K[3]*K[3]; K[6] = K[3]*K[3]; fR = sigma3[itype][jtype]*a12[itype][jtype]*c2*K[1]/K[3]; fforce = 4.0/15.0*r*fR*factor_lj * (2.0*(K[1]+K[2])*(K[1]*(5.0*K[1]+22.0*K[2])+5.0*K[4]) * sigma6[itype][jtype]/K[6] - 5.0)/K[0]; phi = 2.0/9.0*fR * (1.0-(K[1]*(K[1]*(K[1]/3.0+3.0*K[2])+4.2*K[4])+K[2]*K[4]) * sigma6[itype][jtype]/K[6]) - offset[itype][jtype]; break; case LARGE_LARGE: r = sqrt(rsq); c1 = a1[itype][jtype]; c2 = a2[itype][jtype]; K[0] = c1*c2; K[1] = c1+c2; K[2] = c1-c2; K[3] = K[1]+r; K[4] = K[1]-r; K[5] = K[2]+r; K[6] = K[2]-r; K[7] = 1.0/(K[3]*K[4]); K[8] = 1.0/(K[5]*K[6]); g[0] = pow(K[3],-7.0); g[1] = pow(K[4],-7.0); g[2] = pow(K[5],-7.0); g[3] = pow(K[6],-7.0); h[0] = ((K[3]+5.0*K[1])*K[3]+30.0*K[0])*g[0]; h[1] = ((K[4]+5.0*K[1])*K[4]+30.0*K[0])*g[1]; h[2] = ((K[5]+5.0*K[2])*K[5]-30.0*K[0])*g[2]; h[3] = ((K[6]+5.0*K[2])*K[6]-30.0*K[0])*g[3]; g[0] *= 42.0*K[0]/K[3]+6.0*K[1]+K[3]; g[1] *= 42.0*K[0]/K[4]+6.0*K[1]+K[4]; g[2] *= -42.0*K[0]/K[5]+6.0*K[2]+K[5]; g[3] *= -42.0*K[0]/K[6]+6.0*K[2]+K[6]; fR = a12[itype][jtype]*sigma6[itype][jtype]/r/37800.0; phi = fR * (h[0]-h[1]-h[2]+h[3]); dUR = phi/r + 5.0*fR*(g[0]+g[1]-g[2]-g[3]); dUA = -a12[itype][jtype]/3.0*r*((2.0*K[0]*K[7]+1.0)*K[7] + (2.0*K[0]*K[8]-1.0)*K[8]); fforce = factor_lj*(dUR+dUA)/r; phi += a12[itype][jtype]/6.0*(2.0*K[0]*(K[7]+K[8])-log(K[8]/K[7])) - offset[itype][jtype]; break; } return factor_lj*phi; } + +/* ---------------------------------------------------------------------- */ + +void *PairColloid::extract(char *str) +{ + if (strcmp(str,"diameter") == 0) return (void *) diameter; + return NULL; +} diff --git a/src/COLLOID/pair_colloid.h b/src/COLLOID/pair_colloid.h index 87e10ef86..00f41b4b3 100644 --- a/src/COLLOID/pair_colloid.h +++ b/src/COLLOID/pair_colloid.h @@ -1,48 +1,49 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator www.cs.sandia.gov/~sjplimp/lammps.html Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories 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 PAIR_COLLOID_H #define PAIR_COLLOID_H #include "pair.h" namespace LAMMPS_NS { class PairColloid : public Pair { public: PairColloid(class LAMMPS *); ~PairColloid(); void compute(int, int); void settings(int, char **); void coeff(int, char **); double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); double single(int, int, int, int, double, double, double, double &); + void *extract(char *); private: double cut_global; double **cut; - double **a12,**d1,**d2,**a1,**a2,**offset; + double **a12,**d1,**d2,**diameter,**a1,**a2,**offset; double **sigma,**sigma3,**sigma6; double **lj1,**lj2,**lj3,**lj4; int **form; void allocate(); }; } #endif diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp index b2323d2c9..78c7088fa 100644 --- a/src/COLLOID/pair_lubricate.cpp +++ b/src/COLLOID/pair_lubricate.cpp @@ -1,480 +1,478 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Randy Schunk (SNL) ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_lubricate.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "update.h" #include "memory.h" #include "error.h" #include "domain.h" using namespace LAMMPS_NS; #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) /* ---------------------------------------------------------------------- */ PairLubricate::PairLubricate(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; } /* ---------------------------------------------------------------------- */ PairLubricate::~PairLubricate() { if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); memory->destroy_2d_double_array(cut); memory->destroy_2d_double_array(cut_inner); } } /* ---------------------------------------------------------------------- */ void PairLubricate::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz; double rsq,r,h_sep,radi; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; double vt1,vt2,vt3,w1,w2,w3,v_shear1,v_shear2,v_shear3; double omega_t_1,omega_t_2,omega_t_3; double n_cross_omega_t_1,n_cross_omega_t_2,n_cross_omega_t_3; double wr1,wr2,wr3,wnnr,wn1,wn2,wn3,inv_inertia; double P_dot_wrel_1,P_dot_wrel_2,P_dot_wrel_3; double a_squeeze,a_shear,a_pump,a_twist; int *ilist,*jlist,*numneigh,**firstneigh; double PI = 4.0*atan(1.0); if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **v = atom->v; double **f = atom->f; double **omega = atom->omega; double **angmom = atom->angmom; double **torque = atom->torque; int *type = atom->type; int nlocal = atom->nlocal; int newton_pair = force->newton_pair; int omega_flag = atom->omega_flag; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; a_squeeze = a_shear = a_pump = a_twist = 0.0; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; radi = atom->shape[itype][0]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { r = sqrt(rsq); // relative translational velocity vr1 = v[i][0] - v[j][0]; vr2 = v[i][1] - v[j][1]; vr3 = v[i][2] - v[j][2]; // normal component N.(v1-v2) = nn.(v1-v2) vnnr = vr1*delx + vr2*dely + vr3*delz; vnnr /= r; vn1 = delx*vnnr / r; vn2 = dely*vnnr / r; vn3 = delz*vnnr / r; // tangential component -P.(v1-v2) // P = (I - nn) where n is vector between centers vt1 = vr1 - vn1; vt2 = vr2 - vn2; vt3 = vr3 - vn3; // additive rotational velocity = omega_1 + omega_2 // use omega directly if it exists, else angmom // angular momentum = I*omega = 2/5 * M*R^2 * omega if (omega_flag) { w1 = omega[i][0] + omega[j][0]; w2 = omega[i][1] + omega[j][1]; w3 = omega[i][2] + omega[j][2]; } else { inv_inertia = 1.0 / (0.4*atom->mass[itype]*radi*radi); w1 = inv_inertia * (angmom[i][0] + angmom[j][0]); w1 = inv_inertia * (angmom[i][0] + angmom[j][0]); w1 = inv_inertia * (angmom[i][0] + angmom[j][0]); } // relative velocities n X P . (v1-v2) = n X (I-nn) . (v1-v2) v_shear1 = (dely*vt3 - delz*vt2) / r; v_shear2 = -(delx*vt3 - delz*vt1) / r; v_shear3 = (delx*vt2 - dely*vt1) / r; // relative rotation rate P.(omega1 + omega2) omega_t_1 = w1 - delx*(delx*w1) / rsq; omega_t_2 = w2 - dely*(dely*w2) / rsq; omega_t_3 = w3 - delz*(delz*w3) / rsq; // n x omega_t n_cross_omega_t_1 = (dely*omega_t_3 - delz*omega_t_2) / r; n_cross_omega_t_2 = -(delx*omega_t_3 - delz*omega_t_1) / r; n_cross_omega_t_3 = (delx*omega_t_2 - dely*omega_t_1) / r; // N.(w1-w2) and P.(w1-w2) if (omega_flag) { wr1 = omega[i][0] - omega[j][0]; wr2 = omega[i][1] - omega[j][1]; wr3 = omega[i][2] - omega[j][2]; } else { wr1 = inv_inertia * (angmom[i][0] - angmom[j][0]); wr1 = inv_inertia * (angmom[i][0] - angmom[j][0]); wr1 = inv_inertia * (angmom[i][0] - angmom[j][0]); } wnnr = wr1*delx + wr2*dely + wr3*delz; wn1 = delx*wnnr / rsq; wn2 = dely*wnnr / rsq; wn3 = delz*wnnr / rsq; P_dot_wrel_1 = wr1 - delx*(delx*wr1)/rsq; P_dot_wrel_2 = wr2 - dely*(dely*wr2)/rsq; P_dot_wrel_3 = wr3 - delz*(delz*wr3)/rsq; // compute components of pair-hydro h_sep = r - 2.0*radi; if (flag1) a_squeeze = (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep); if (flag2) a_shear = (PI*mu*2.*radi/2.0) * log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0; if (flag3) a_pump = (PI*mu*pow(2.0*radi,4)/8.0) * ((3.0/20.0) * log(2.0*radi/2.0/h_sep) + (63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep)); if (flag4) a_twist = (PI*mu*pow(2.0*radi,4)/4.0) * (h_sep/2.0/radi) * log(2.0/(2.0*h_sep)); if (h_sep >= cut_inner[itype][jtype]) { fx = -a_squeeze*vn1 - a_shear*(2.0/r)*(2.0/r)*vt1 + (2.0/r)*a_shear*n_cross_omega_t_1; fy = -a_squeeze*vn2 - a_shear*(2.0/r)*(2.0/r)*vt2 + (2.0/r)*a_shear*n_cross_omega_t_2; fz = -a_squeeze*vn3 - a_shear*(2.0/r)*(2.0/r)*vt3 + (2.0/r)*a_shear*n_cross_omega_t_3; torque[i][0] += -(2.0/r)*a_shear*v_shear1 - a_shear*omega_t_1 - a_pump*P_dot_wrel_1 - a_twist*wn1; torque[i][1] += -(2.0/r)*a_shear*v_shear2 - a_shear*omega_t_2 - a_pump*P_dot_wrel_2 - a_twist*wn2; torque[i][2] += -(2.0/r)*a_shear*v_shear3 - a_shear*omega_t_3 - a_pump*P_dot_wrel_3 - a_twist*wn3; } else { fpair = -vnnr*(3.0*PI*mu*radi/2.0)*radi/4.0/cut_inner[itype][jtype]; fx = fpair*delx/r; fy = fpair*dely/r; fz = fpair*delz/r; } f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; if (newton_pair || j < nlocal) { f[j][0] -= fx; f[j][1] -= fy; f[j][2] -= fz; if (h_sep >= cut_inner[itype][jtype]) { torque[j][0] += -(2.0/r)*a_shear*v_shear1 - a_shear*omega_t_1 + a_pump*P_dot_wrel_1 + a_twist*wn1; torque[j][1] += -(2.0/r)*a_shear*v_shear2 - a_shear*omega_t_2 + a_pump*P_dot_wrel_2 + a_twist*wn2; torque[j][2] += -(2.0/r)*a_shear*v_shear3 - a_shear*omega_t_3 + a_pump*P_dot_wrel_3 + a_twist*wn3; } } if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, 0.0,0.0,fx,fy,fz,delx,dely,delz); } } } if (vflag_fdotr) virial_compute(); } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairLubricate::allocate() { allocated = 1; int n = atom->ntypes; setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); cut_inner = memory->create_2d_double_array(n+1,n+1,"pair:cut_inner"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairLubricate::settings(int narg, char **arg) { if (narg != 7) error->all("Illegal pair_style command"); mu = atof(arg[0]); flag1 = atoi(arg[1]); flag2 = atoi(arg[2]); flag3 = atoi(arg[3]); flag4 = atoi(arg[4]); cut_inner_global = atof(arg[5]); cut_global = atof(arg[6]); // reset cutoffs that have been explicitly set if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) { cut_inner[i][j] = cut_inner_global; cut[i][j] = cut_global; } } } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairLubricate::coeff(int narg, char **arg) { if (narg != 2 && narg != 4) error->all("Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); double cut_inner_one = cut_inner_global; double cut_one = cut_global; if (narg == 4) { cut_inner_one = atof(arg[2]); cut_one = atof(arg[3]); } int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { cut_inner[i][j] = cut_inner_one; cut[i][j] = cut_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all("Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairLubricate::init_style() { if (!atom->torque_flag || atom->shape == NULL) error->all("Pair lubricate requires atom attributes torque, shape"); - if (domain->dimension != 3) - error->all("Pair lubricate only available for 3d"); // insure all particle shapes are spherical for (int i = 1; i <= atom->ntypes; i++) if ((atom->shape[i][0] != atom->shape[i][1]) || (atom->shape[i][0] != atom->shape[i][2]) || (atom->shape[i][1] != atom->shape[i][2]) ) error->all("Pair lubricate requires spherical particles"); // insure mono-dispersity for (int i = 2; i <= atom->ntypes; i++) if (atom->shape[i][0] != atom->shape[1][1]) error->all("Pair lubricate requires mono-disperse particles"); int irequest = neighbor->request(this); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairLubricate::init_one(int i, int j) { if (setflag[i][j] == 0) { cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); cut[i][j] = mix_distance(cut[i][i],cut[j][j]); } cut_inner[j][i] = cut_inner[i][j]; return cut[i][j]; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLubricate::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&cut_inner[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLubricate::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) { fread(&cut_inner[i][j],sizeof(double),1,fp); fread(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLubricate::write_restart_settings(FILE *fp) { fwrite(&mu,sizeof(double),1,fp); fwrite(&flag1,sizeof(int),1,fp); fwrite(&flag2,sizeof(int),1,fp); fwrite(&flag3,sizeof(int),1,fp); fwrite(&flag4,sizeof(int),1,fp); fwrite(&cut_inner_global,sizeof(double),1,fp); fwrite(&cut_global,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLubricate::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { fread(&mu,sizeof(double),1,fp); fread(&flag1,sizeof(int),1,fp); fread(&flag2,sizeof(int),1,fp); fread(&flag3,sizeof(int),1,fp); fread(&flag4,sizeof(int),1,fp); fread(&cut_inner_global,sizeof(double),1,fp); fread(&cut_global,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&mu,1,MPI_DOUBLE,0,world); MPI_Bcast(&flag1,1,MPI_INT,0,world); MPI_Bcast(&flag2,1,MPI_INT,0,world); MPI_Bcast(&flag3,1,MPI_INT,0,world); MPI_Bcast(&flag4,1,MPI_INT,0,world); MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } diff --git a/src/comm.cpp b/src/comm.cpp index 414f70c15..51f74e9b0 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -1,1489 +1,1528 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author (triclinic) : Pieter in 't Veld (SNL) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "stdio.h" #include "stdlib.h" #include "comm.h" #include "atom.h" #include "atom_vec.h" #include "force.h" #include "pair.h" #include "domain.h" #include "neighbor.h" #include "modify.h" #include "fix.h" +#include "group.h" #include "compute.h" #include "error.h" #include "memory.h" using namespace LAMMPS_NS; #define BUFFACTOR 1.5 #define BUFMIN 1000 #define BUFEXTRA 1000 #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define BIG 1.0e20 enum{SINGLE,MULTI}; /* ---------------------------------------------------------------------- setup MPI and allocate buffer space ------------------------------------------------------------------------- */ Comm::Comm(LAMMPS *lmp) : Pointers(lmp) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); user_procgrid[0] = user_procgrid[1] = user_procgrid[2] = 0; grid2proc = NULL; + igroup = 0; style = SINGLE; multilo = multihi = NULL; cutghostmulti = NULL; // initialize comm buffers & exchange memory maxsend = BUFMIN; buf_send = (double *) memory->smalloc((maxsend+BUFEXTRA)*sizeof(double),"comm:buf_send"); maxrecv = BUFMIN; buf_recv = (double *) memory->smalloc(maxrecv*sizeof(double),"comm:buf_recv"); maxswap = 6; allocate_swap(maxswap); sendlist = (int **) memory->smalloc(maxswap*sizeof(int *),"sendlist"); maxsendlist = (int *) memory->smalloc(maxswap*sizeof(int),"maxsendlist"); for (int i = 0; i < maxswap; i++) { maxsendlist[i] = BUFMIN; sendlist[i] = (int *) memory->smalloc(BUFMIN*sizeof(int),"sendlist[i]"); } maxforward_fix = maxreverse_fix = 0; maxforward_pair = maxreverse_pair = 0; } /* ---------------------------------------------------------------------- */ Comm::~Comm() { if (grid2proc) memory->destroy_3d_int_array(grid2proc); free_swap(); if (style == MULTI) { free_multi(); memory->destroy_2d_double_array(cutghostmulti); } memory->sfree(maxsendlist); if (sendlist) for (int i = 0; i < maxswap; i++) memory->sfree(sendlist[i]); memory->sfree(sendlist); memory->sfree(buf_send); memory->sfree(buf_recv); } /* ---------------------------------------------------------------------- setup 3d grid of procs based on box size ------------------------------------------------------------------------- */ void Comm::set_procs() { if (user_procgrid[0] == 0) procs2box(); else { procgrid[0] = user_procgrid[0]; procgrid[1] = user_procgrid[1]; procgrid[2] = user_procgrid[2]; } if (procgrid[0]*procgrid[1]*procgrid[2] != nprocs) error->all("Bad grid of processors"); if (domain->dimension == 2 && procgrid[2] != 1) error->all("Proc grid in z != 1 for 2d simulation"); if (grid2proc) memory->destroy_3d_int_array(grid2proc); grid2proc = memory->create_3d_int_array(procgrid[0],procgrid[1],procgrid[2], "comm:grid2proc"); // use MPI Cartesian routines to setup 3d grid of procs // grid2proc[i][j][k] = proc that owns i,j,k location in grid // let MPI compute it instead of LAMMPS in case it is machine optimized int reorder = 0; int periods[3]; periods[0] = periods[1] = periods[2] = 1; MPI_Comm cartesian; MPI_Cart_create(world,3,procgrid,periods,reorder,&cartesian); MPI_Cart_get(cartesian,3,procgrid,periods,myloc); MPI_Cart_shift(cartesian,0,1,&procneigh[0][0],&procneigh[0][1]); MPI_Cart_shift(cartesian,1,1,&procneigh[1][0],&procneigh[1][1]); MPI_Cart_shift(cartesian,2,1,&procneigh[2][0],&procneigh[2][1]); int coords[3]; int i,j,k; for (i = 0; i < procgrid[0]; i++) for (j = 0; j < procgrid[1]; j++) for (k = 0; k < procgrid[2]; k++) { coords[0] = i; coords[1] = j; coords[2] = k; MPI_Cart_rank(cartesian,coords,&grid2proc[i][j][k]); } MPI_Comm_free(&cartesian); // set lamda box params after procs are assigned if (domain->triclinic) domain->set_lamda_box(); if (me == 0) { if (screen) fprintf(screen," %d by %d by %d processor grid\n", procgrid[0],procgrid[1],procgrid[2]); if (logfile) fprintf(logfile," %d by %d by %d processor grid\n", procgrid[0],procgrid[1],procgrid[2]); } } /* ---------------------------------------------------------------------- */ void Comm::init() { triclinic = domain->triclinic; map_style = atom->map_style; + groupbit = group->bitmask[igroup]; // comm_only = 1 if only x,f are exchanged in forward/reverse comm comm_x_only = atom->avec->comm_x_only; comm_f_only = atom->avec->comm_f_only; // maxforward = # of datums in largest forward communication // maxreverse = # of datums in largest reverse communication // query pair,fix,compute for their requirements maxforward = MAX(atom->avec->size_comm,atom->avec->size_border); maxreverse = atom->avec->size_reverse; if (force->pair) maxforward = MAX(maxforward,force->pair->comm_forward); if (force->pair) maxreverse = MAX(maxreverse,force->pair->comm_reverse); for (int i = 0; i < modify->nfix; i++) { maxforward = MAX(maxforward,modify->fix[i]->comm_forward); maxreverse = MAX(maxreverse,modify->fix[i]->comm_reverse); } for (int i = 0; i < modify->ncompute; i++) { maxforward = MAX(maxforward,modify->compute[i]->comm_forward); maxreverse = MAX(maxreverse,modify->compute[i]->comm_reverse); } if (force->newton == 0) maxreverse = 0; // memory for multi-style communication if (style == MULTI && multilo == NULL) { allocate_multi(maxswap); cutghostmulti = memory->create_2d_double_array(atom->ntypes+1,3,"comm:cutghostmulti"); } if (style == SINGLE && multilo) { free_multi(); memory->destroy_2d_double_array(cutghostmulti); } } /* ---------------------------------------------------------------------- setup spatial-decomposition communication patterns function of neighbor cutoff(s) and current box size single style sets slab boundaries (slablo,slabhi) based on max cutoff multi style sets type-dependent slab boundaries (multilo,multihi) ------------------------------------------------------------------------- */ void Comm::setup() { // cutghost = max distance at which ghost atoms need to be acquired // for orthogonal: // cutghost is in box coords = neigh->cutghost in all 3 dims // for triclinic: // neigh->cutghost = distance between tilted planes in box coords // cutghost is in lamda coords = distance between those planes // for multi: // cutghostmulti = same as cutghost, only for each atom type int i; int ntypes = atom->ntypes; double *prd,*sublo,*subhi; if (triclinic == 0) { prd = domain->prd; sublo = domain->sublo; subhi = domain->subhi; cutghost[0] = cutghost[1] = cutghost[2] = neighbor->cutghost; if (style == MULTI) { double *cuttype = neighbor->cuttype; for (i = 1; i <= ntypes; i++) cutghostmulti[i][0] = cutghostmulti[i][1] = cutghostmulti[i][2] = cuttype[i]; } } else { prd = domain->prd_lamda; sublo = domain->sublo_lamda; subhi = domain->subhi_lamda; double *h_inv = domain->h_inv; double length0,length1,length2; length0 = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]); cutghost[0] = neighbor->cutghost * length0; length1 = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]); cutghost[1] = neighbor->cutghost * length1; length2 = h_inv[2]; cutghost[2] = neighbor->cutghost * length2; if (style == MULTI) { double *cuttype = neighbor->cuttype; for (i = 1; i <= ntypes; i++) { cutghostmulti[i][0] = cuttype[i] * length0; cutghostmulti[i][1] = cuttype[i] * length1; cutghostmulti[i][2] = cuttype[i] * length2; } } } // need = # of procs I need atoms from in each dim based on max cutoff // for 2d, don't communicate in z need[0] = static_cast (cutghost[0] * procgrid[0] / prd[0]) + 1; need[1] = static_cast (cutghost[1] * procgrid[1] / prd[1]) + 1; need[2] = static_cast (cutghost[2] * procgrid[2] / prd[2]) + 1; if (domain->dimension == 2) need[2] = 0; // if non-periodic, do not communicate further than procgrid-1 away // this enables very large cutoffs in non-periodic systems int *periodicity = domain->periodicity; if (periodicity[0] == 0) need[0] = MIN(need[0],procgrid[0]-1); if (periodicity[1] == 0) need[1] = MIN(need[1],procgrid[1]-1); if (periodicity[2] == 0) need[2] = MIN(need[2],procgrid[2]-1); // allocate comm memory nswap = 2 * (need[0]+need[1]+need[2]); if (nswap > maxswap) grow_swap(nswap); // setup parameters for each exchange: // sendproc = proc to send to at each swap // recvproc = proc to recv from at each swap // for style SINGLE: // slablo/slabhi = boundaries for slab of atoms to send at each swap // use -BIG/midpt/BIG to insure all atoms included even if round-off occurs // if round-off, atoms recvd across PBC can be < or > than subbox boundary // note that borders() only loops over subset of atoms during each swap // set slablo > slabhi for swaps across non-periodic boundaries // this insures no atoms are swapped // only for procs owning sub-box at non-periodic end of global box // for style MULTI: // multilo/multihi is same as slablo/slabhi, only for each atom type // pbc_flag: 0 = nothing across a boundary, 1 = something across a boundary // pbc = -1/0/1 for PBC factor in each of 3/6 orthog/triclinic dirs // for triclinic, slablo/hi and pbc_border will be used in lamda (0-1) coords // 1st part of if statement is sending to the west/south/down // 2nd part of if statement is sending to the east/north/up int dim,ineed; int iswap = 0; for (dim = 0; dim < 3; dim++) { for (ineed = 0; ineed < 2*need[dim]; ineed++) { pbc_flag[iswap] = 0; pbc[iswap][0] = pbc[iswap][1] = pbc[iswap][2] = pbc[iswap][3] = pbc[iswap][4] = pbc[iswap][5] = 0; if (ineed % 2 == 0) { sendproc[iswap] = procneigh[dim][0]; recvproc[iswap] = procneigh[dim][1]; if (style == SINGLE) { if (ineed < 2) slablo[iswap] = -BIG; else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]); slabhi[iswap] = sublo[dim] + cutghost[dim]; } else { for (i = 1; i <= ntypes; i++) { if (ineed < 2) multilo[iswap][i] = -BIG; else multilo[iswap][i] = 0.5 * (sublo[dim] + subhi[dim]); multihi[iswap][i] = sublo[dim] + cutghostmulti[i][dim]; } } if (myloc[dim] == 0) { if (periodicity[dim] == 0) { if (style == SINGLE) slabhi[iswap] = slablo[iswap] - 1.0; else for (i = 1; i <= ntypes; i++) multihi[iswap][i] = multilo[iswap][i] - 1.0; } else { pbc_flag[iswap] = 1; pbc[iswap][dim] = 1; if (triclinic) { if (dim == 1) pbc[iswap][5] = 1; else if (dim == 2) pbc[iswap][4] = pbc[iswap][3] = 1; } } } } else { sendproc[iswap] = procneigh[dim][1]; recvproc[iswap] = procneigh[dim][0]; if (style == SINGLE) { slablo[iswap] = subhi[dim] - cutghost[dim]; if (ineed < 2) slabhi[iswap] = BIG; else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]); } else { for (i = 1; i <= ntypes; i++) { multilo[iswap][i] = subhi[dim] - cutghostmulti[i][dim]; if (ineed < 2) multihi[iswap][i] = BIG; else multihi[iswap][i] = 0.5 * (sublo[dim] + subhi[dim]); } } if (myloc[dim] == procgrid[dim]-1) { if (periodicity[dim] == 0) { if (style == SINGLE) slabhi[iswap] = slablo[iswap] - 1.0; else for (i = 1; i <= ntypes; i++) multihi[iswap][i] = multilo[iswap][i] - 1.0; } else { pbc_flag[iswap] = 1; pbc[iswap][dim] = -1; if (triclinic) { if (dim == 1) pbc[iswap][5] = -1; else if (dim == 2) pbc[iswap][4] = pbc[iswap][3] = -1; } } } } iswap++; } } } /* ---------------------------------------------------------------------- communication of atom coords every timestep other stuff may also be sent via pack/unpack routines ------------------------------------------------------------------------- */ void Comm::communicate() { int n; MPI_Request request; MPI_Status status; AtomVec *avec = atom->avec; double **x = atom->x; double *buf; // exchange data with another proc // if other proc is self, just copy // if comm_x_only set, exchange or copy directly to x, don't unpack for (int iswap = 0; iswap < nswap; iswap++) { if (sendproc[iswap] != me) { if (comm_x_only) { if (size_comm_recv[iswap]) buf = x[firstrecv[iswap]]; else buf = NULL; MPI_Irecv(buf,size_comm_recv[iswap],MPI_DOUBLE, recvproc[iswap],0,world,&request); n = avec->pack_comm(sendnum[iswap],sendlist[iswap], buf_send,pbc_flag[iswap],pbc[iswap]); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); } else { MPI_Irecv(buf_recv,size_comm_recv[iswap],MPI_DOUBLE, recvproc[iswap],0,world,&request); n = avec->pack_comm(sendnum[iswap],sendlist[iswap], buf_send,pbc_flag[iswap],pbc[iswap]); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_recv); } } else { if (comm_x_only) { if (sendnum[iswap]) n = avec->pack_comm(sendnum[iswap],sendlist[iswap], x[firstrecv[iswap]],pbc_flag[iswap],pbc[iswap]); } else { n = avec->pack_comm(sendnum[iswap],sendlist[iswap], buf_send,pbc_flag[iswap],pbc[iswap]); avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send); } } } } /* ---------------------------------------------------------------------- reverse communication of forces on atoms every timestep other stuff can also be sent via pack/unpack routines ------------------------------------------------------------------------- */ void Comm::reverse_communicate() { int n; MPI_Request request; MPI_Status status; AtomVec *avec = atom->avec; double **f = atom->f; double *buf; // exchange data with another proc // if other proc is self, just copy // if comm_f_only set, exchange or copy directly from f, don't pack for (int iswap = nswap-1; iswap >= 0; iswap--) { if (sendproc[iswap] != me) { if (comm_f_only) { MPI_Irecv(buf_recv,size_reverse_recv[iswap],MPI_DOUBLE, sendproc[iswap],0,world,&request); if (size_reverse_send[iswap]) buf = f[firstrecv[iswap]]; else buf = NULL; MPI_Send(buf,size_reverse_send[iswap],MPI_DOUBLE, recvproc[iswap],0,world); MPI_Wait(&request,&status); } else { MPI_Irecv(buf_recv,size_reverse_recv[iswap],MPI_DOUBLE, sendproc[iswap],0,world,&request); n = avec->pack_reverse(recvnum[iswap],firstrecv[iswap],buf_send); MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap],0,world); MPI_Wait(&request,&status); } avec->unpack_reverse(sendnum[iswap],sendlist[iswap],buf_recv); } else { if (comm_f_only) { if (sendnum[iswap]) avec->unpack_reverse(sendnum[iswap],sendlist[iswap], f[firstrecv[iswap]]); } else { n = avec->pack_reverse(recvnum[iswap],firstrecv[iswap],buf_send); avec->unpack_reverse(sendnum[iswap],sendlist[iswap],buf_send); } } } } /* ---------------------------------------------------------------------- exchange: move atoms to correct processors atoms exchanged with all 6 stencil neighbors send out atoms that have left my box, receive ones entering my box atoms will be lost if not inside some proc's box can happen if atom moves outside of non-periodic bounary or if atom moves more than one proc away this routine called before every reneighboring for triclinic, atoms must be in lamda coords (0-1) before exchange is called ------------------------------------------------------------------------- */ void Comm::exchange() { int i,m,nsend,nrecv,nrecv1,nrecv2,nlocal; double lo,hi,value; double **x; double *sublo,*subhi,*buf; MPI_Request request; MPI_Status status; AtomVec *avec = atom->avec; // clear global->local map since atoms move & new ghosts are created if (map_style) atom->map_clear(); // subbox bounds for orthogonal or triclinic if (triclinic == 0) { sublo = domain->sublo; subhi = domain->subhi; } else { sublo = domain->sublo_lamda; subhi = domain->subhi_lamda; } // loop over dimensions for (int dim = 0; dim < 3; dim++) { // fill buffer with atoms leaving my box, using < and >= // when atom is deleted, fill it in with last atom x = atom->x; lo = sublo[dim]; hi = subhi[dim]; nlocal = atom->nlocal; i = nsend = 0; while (i < nlocal) { if (x[i][dim] < lo || x[i][dim] >= hi) { if (nsend > maxsend) grow_send(nsend,1); nsend += avec->pack_exchange(i,&buf_send[nsend]); avec->copy(nlocal-1,i); nlocal--; } else i++; } atom->nlocal = nlocal; // send/recv atoms in both directions // if 1 proc in dimension, no send/recv, set recv buf to send buf // if 2 procs in dimension, single send/recv // if more than 2 procs in dimension, send/recv to both neighbors if (procgrid[dim] == 1) { nrecv = nsend; buf = buf_send; } else { MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][0],0, &nrecv1,1,MPI_INT,procneigh[dim][1],0,world,&status); nrecv = nrecv1; if (procgrid[dim] > 2) { MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][1],0, &nrecv2,1,MPI_INT,procneigh[dim][0],0,world,&status); nrecv += nrecv2; } if (nrecv > maxrecv) grow_recv(nrecv); MPI_Irecv(buf_recv,nrecv1,MPI_DOUBLE,procneigh[dim][1],0, world,&request); MPI_Send(buf_send,nsend,MPI_DOUBLE,procneigh[dim][0],0,world); MPI_Wait(&request,&status); if (procgrid[dim] > 2) { MPI_Irecv(&buf_recv[nrecv1],nrecv2,MPI_DOUBLE,procneigh[dim][0],0, world,&request); MPI_Send(buf_send,nsend,MPI_DOUBLE,procneigh[dim][1],0,world); MPI_Wait(&request,&status); } buf = buf_recv; } // check incoming atoms to see if they are in my box // if so, add to my list m = 0; while (m < nrecv) { value = buf[m+dim+1]; if (value >= lo && value < hi) m += avec->unpack_exchange(&buf[m]); else m += static_cast (buf[m]); } } } /* ---------------------------------------------------------------------- borders: make lists of nearby atoms to send to neighboring procs at every timestep one list is created for every swap that will be made as list is made, actually do swaps this does equivalent of a communicate (so don't need to explicitly call communicate routine on reneighboring timestep) this routine is called before every reneighboring for triclinic, atoms must be in lamda coords (0-1) before borders is called ------------------------------------------------------------------------- */ void Comm::borders() { int i,n,itype,iswap,dim,ineed,maxneed,nsend,nrecv,nfirst,nlast,smax,rmax; double lo,hi; - int *type; + int *type,*mask; double **x; double *buf,*mlo,*mhi; MPI_Request request; MPI_Status status; AtomVec *avec = atom->avec; int size_border = avec->size_border; // clear old ghosts atom->nghost = 0; // do swaps over all 3 dimensions iswap = 0; smax = rmax = 0; for (dim = 0; dim < 3; dim++) { nlast = 0; maxneed = 2*need[dim]; for (ineed = 0; ineed < maxneed; ineed++) { // find atoms within slab boundaries lo/hi using <= and >= // check atoms between nfirst and nlast // for first swaps in a dim, check owned and ghost // for later swaps in a dim, only check newly arrived ghosts // store sent atom indices in list for use in future timesteps x = atom->x; if (style == SINGLE) { lo = slablo[iswap]; hi = slabhi[iswap]; } else { type = atom->type; mlo = multilo[iswap]; mhi = multihi[iswap]; } if (ineed % 2 == 0) { nfirst = nlast; nlast = atom->nlocal + atom->nghost; } nsend = 0; - if (style == SINGLE) { - for (i = nfirst; i < nlast; i++) - if (x[i][dim] >= lo && x[i][dim] <= hi) { - if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend); - sendlist[iswap][nsend++] = i; - } + + // SINGLE vs MULTI, all atoms versus only atoms in group + + if (igroup) { + mask = atom->mask; + if (style == SINGLE) { + for (i = nfirst; i < nlast; i++) + if (mask[i] & groupbit) { + if (x[i][dim] >= lo && x[i][dim] <= hi) { + if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend); + sendlist[iswap][nsend++] = i; + } + } + } else { + for (i = nfirst; i < nlast; i++) + if (mask[i] & groupbit) { + itype = type[i]; + if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) { + if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend); + sendlist[iswap][nsend++] = i; + } + } + } } else { - for (i = nfirst; i < nlast; i++) { - itype = type[i]; - if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) { - if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend); - sendlist[iswap][nsend++] = i; + if (style == SINGLE) { + for (i = nfirst; i < nlast; i++) + if (x[i][dim] >= lo && x[i][dim] <= hi) { + if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend); + sendlist[iswap][nsend++] = i; + } + } else { + for (i = nfirst; i < nlast; i++) { + itype = type[i]; + if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) { + if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend); + sendlist[iswap][nsend++] = i; + } } } } - + // pack up list of border atoms if (nsend*size_border > maxsend) grow_send(nsend*size_border,0); n = avec->pack_border(nsend,sendlist[iswap],buf_send, pbc_flag[iswap],pbc[iswap]); // swap atoms with other proc // put incoming ghosts at end of my atom arrays // if swapping with self, simply copy, no messages if (sendproc[iswap] != me) { MPI_Sendrecv(&nsend,1,MPI_INT,sendproc[iswap],0, &nrecv,1,MPI_INT,recvproc[iswap],0,world,&status); if (nrecv*size_border > maxrecv) grow_recv(nrecv*size_border); MPI_Irecv(buf_recv,nrecv*size_border,MPI_DOUBLE, recvproc[iswap],0,world,&request); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); buf = buf_recv; } else { nrecv = nsend; buf = buf_send; } // unpack buffer avec->unpack_border(nrecv,atom->nlocal+atom->nghost,buf); // set all pointers & counters smax = MAX(smax,nsend); rmax = MAX(rmax,nrecv); sendnum[iswap] = nsend; recvnum[iswap] = nrecv; size_comm_recv[iswap] = nrecv * avec->size_comm; size_reverse_send[iswap] = nrecv * avec->size_reverse; size_reverse_recv[iswap] = nsend * avec->size_reverse; firstrecv[iswap] = atom->nlocal + atom->nghost; atom->nghost += nrecv; iswap++; } } // insure send/recv buffers are long enough for all forward & reverse comm int max = MAX(maxforward*smax,maxreverse*rmax); if (max > maxsend) grow_send(max,0); max = MAX(maxforward*rmax,maxreverse*smax); if (max > maxrecv) grow_recv(max); // reset global->local map if (map_style) atom->map_set(); } /* ---------------------------------------------------------------------- forward communication invoked by a Pair ------------------------------------------------------------------------- */ void Comm::comm_pair(Pair *pair) { int iswap,n; double *buf; MPI_Request request; MPI_Status status; for (iswap = 0; iswap < nswap; iswap++) { // pack buffer n = pair->pack_comm(sendnum[iswap],sendlist[iswap], buf_send,pbc_flag[iswap],pbc[iswap]); // exchange with another proc // if self, set recv buffer to send buffer if (sendproc[iswap] != me) { MPI_Irecv(buf_recv,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0, world,&request); MPI_Send(buf_send,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); buf = buf_recv; } else buf = buf_send; // unpack buffer pair->unpack_comm(recvnum[iswap],firstrecv[iswap],buf); } } /* ---------------------------------------------------------------------- reverse communication invoked by a Pair ------------------------------------------------------------------------- */ void Comm::reverse_comm_pair(Pair *pair) { int iswap,n; double *buf; MPI_Request request; MPI_Status status; for (iswap = nswap-1; iswap >= 0; iswap--) { // pack buffer n = pair->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send); // exchange with another proc // if self, set recv buffer to send buffer if (sendproc[iswap] != me) { MPI_Irecv(buf_recv,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0, world,&request); MPI_Send(buf_send,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,world); MPI_Wait(&request,&status); buf = buf_recv; } else buf = buf_send; // unpack buffer pair->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf); } } /* ---------------------------------------------------------------------- forward communication invoked by a Fix ------------------------------------------------------------------------- */ void Comm::comm_fix(Fix *fix) { int iswap,n; double *buf; MPI_Request request; MPI_Status status; for (iswap = 0; iswap < nswap; iswap++) { // pack buffer n = fix->pack_comm(sendnum[iswap],sendlist[iswap], buf_send,pbc_flag[iswap],pbc[iswap]); // exchange with another proc // if self, set recv buffer to send buffer if (sendproc[iswap] != me) { MPI_Irecv(buf_recv,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0, world,&request); MPI_Send(buf_send,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); buf = buf_recv; } else buf = buf_send; // unpack buffer fix->unpack_comm(recvnum[iswap],firstrecv[iswap],buf); } } /* ---------------------------------------------------------------------- reverse communication invoked by a Fix ------------------------------------------------------------------------- */ void Comm::reverse_comm_fix(Fix *fix) { int iswap,n; double *buf; MPI_Request request; MPI_Status status; for (iswap = nswap-1; iswap >= 0; iswap--) { // pack buffer n = fix->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send); // exchange with another proc // if self, set recv buffer to send buffer if (sendproc[iswap] != me) { MPI_Irecv(buf_recv,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0, world,&request); MPI_Send(buf_send,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,world); MPI_Wait(&request,&status); buf = buf_recv; } else buf = buf_send; // unpack buffer fix->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf); } } /* ---------------------------------------------------------------------- forward communication invoked by a Compute ------------------------------------------------------------------------- */ void Comm::comm_compute(Compute *compute) { int iswap,n; double *buf; MPI_Request request; MPI_Status status; for (iswap = 0; iswap < nswap; iswap++) { // pack buffer n = compute->pack_comm(sendnum[iswap],sendlist[iswap], buf_send,pbc_flag[iswap],pbc[iswap]); // exchange with another proc // if self, set recv buffer to send buffer if (sendproc[iswap] != me) { MPI_Irecv(buf_recv,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0, world,&request); MPI_Send(buf_send,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); buf = buf_recv; } else buf = buf_send; // unpack buffer compute->unpack_comm(recvnum[iswap],firstrecv[iswap],buf); } } /* ---------------------------------------------------------------------- reverse communication invoked by a Compute ------------------------------------------------------------------------- */ void Comm::reverse_comm_compute(Compute *compute) { int iswap,n; double *buf; MPI_Request request; MPI_Status status; for (iswap = nswap-1; iswap >= 0; iswap--) { // pack buffer n = compute->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send); // exchange with another proc // if self, set recv buffer to send buffer if (sendproc[iswap] != me) { MPI_Irecv(buf_recv,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0, world,&request); MPI_Send(buf_send,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,world); MPI_Wait(&request,&status); buf = buf_recv; } else buf = buf_send; // unpack buffer compute->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf); } } /* ---------------------------------------------------------------------- communicate atoms to new owning procs via irregular communication can be used in place of exchange() unlike exchange(), allows atoms to have moved arbitrarily long distances first setup irregular comm pattern, then invoke it for triclinic, atoms must be in lamda coords (0-1) before irregular is called ------------------------------------------------------------------------- */ void Comm::irregular() { // clear global->local map since atoms move to new procs if (map_style) atom->map_clear(); // subbox bounds for orthogonal or triclinic double *sublo,*subhi; if (triclinic == 0) { sublo = domain->sublo; subhi = domain->subhi; } else { sublo = domain->sublo_lamda; subhi = domain->subhi_lamda; } // loop over atoms, flag any that are not in my sub-box // fill buffer with atoms leaving my box, using < and >= // assign which proc it belongs to via irregular_lookup() // when atom is deleted, fill it in with last atom AtomVec *avec = atom->avec; double **x = atom->x; int nlocal = atom->nlocal; int nsend = 0; int nsendatom = 0; int *sizes = new int[nlocal]; int *proclist = new int[nlocal]; int i = 0; while (i < nlocal) { if (x[i][0] < sublo[0] || x[i][0] >= subhi[0] || x[i][1] < sublo[1] || x[i][1] >= subhi[1] || x[i][2] < sublo[2] || x[i][2] >= subhi[2]) { if (nsend > maxsend) grow_send(nsend,1); sizes[nsendatom] = avec->pack_exchange(i,&buf_send[nsend]); nsend += sizes[nsendatom]; proclist[nsendatom] = irregular_lookup(x[i]); nsendatom++; avec->copy(nlocal-1,i); nlocal--; } else i++; } atom->nlocal = nlocal; // create irregular communication plan, perform comm, destroy plan // returned nrecv = size of buffer needed for incoming atoms int nrecv; Plan *plan = irregular_create(nsendatom,sizes,proclist,&nrecv); if (nrecv > maxrecv) grow_recv(nrecv); irregular_perform(plan,buf_send,sizes,buf_recv); irregular_destroy(plan); delete [] sizes; delete [] proclist; // add received atoms to my list int m = 0; while (m < nrecv) m += avec->unpack_exchange(&buf_recv[m]); // reset global->local map if (map_style) atom->map_set(); } /* ---------------------------------------------------------------------- create an irregular communication plan n = # of atoms to send sizes = # of doubles for each atom proclist = proc to send each atom to (none to me) return nrecvsize = total # of doubles I will recv ------------------------------------------------------------------------- */ Comm::Plan *Comm::irregular_create(int n, int *sizes, int *proclist, int *nrecvsize) { int i; // allocate plan and work vectors Plan *plan = (struct Plan *) memory->smalloc(sizeof(Plan),"comm:plan"); int *list = new int[nprocs]; int *count = new int[nprocs]; // nrecv = # of messages I receive for (i = 0; i < nprocs; i++) { list[i] = 0; count[i] = 1; } for (i = 0; i < n; i++) list[proclist[i]] = 1; int nrecv; MPI_Reduce_scatter(list,&nrecv,count,MPI_INT,MPI_SUM,world); // allocate receive arrays int *proc_recv = new int[nrecv]; int *length_recv = new int[nrecv]; MPI_Request *request = new MPI_Request[nrecv]; MPI_Status *status = new MPI_Status[nrecv]; // nsend = # of messages I send for (i = 0; i < nprocs; i++) list[i] = 0; for (i = 0; i < n; i++) list[proclist[i]] += sizes[i]; int nsend = 0; for (i = 0; i < nprocs; i++) if (list[i]) nsend++; // allocate send arrays int *proc_send = new int[nsend]; int *length_send = new int[nsend]; int *num_send = new int[nsend]; int *index_send = new int[n]; int *offset_send = new int[n]; // list still stores size of message for procs I send to // proc_send = procs I send to // length_send = total size of message I send to each proc // to balance pattern of send messages: // each proc begins with iproc > me, continues until iproc = me // reset list to store which send message each proc corresponds to int iproc = me; int isend = 0; for (i = 0; i < nprocs; i++) { iproc++; if (iproc == nprocs) iproc = 0; if (list[iproc] > 0) { proc_send[isend] = iproc; length_send[isend] = list[iproc]; list[iproc] = isend; isend++; } } // num_send = # of datums I send to each proc for (i = 0; i < nsend; i++) num_send[i] = 0; for (i = 0; i < n; i++) { isend = list[proclist[i]]; num_send[isend]++; } // count = offsets into n-length index_send for each proc I send to // index_send = list of which datums to send to each proc // 1st N1 values are datum indices for 1st proc, // next N2 values are datum indices for 2nd proc, etc // offset_send = where each datum starts in send buffer count[0] = 0; for (i = 1; i < nsend; i++) count[i] = count[i-1] + num_send[i-1]; for (i = 0; i < n; i++) { isend = list[proclist[i]]; index_send[count[isend]++] = i; if (i) offset_send[i] = offset_send[i-1] + sizes[i-1]; else offset_send[i] = 0; } // tell receivers how much data I send // sendmax = largest # of datums I send in a single message int sendmax = 0; for (i = 0; i < nsend; i++) { MPI_Send(&length_send[i],1,MPI_INT,proc_send[i],0,world); sendmax = MAX(sendmax,length_send[i]); } // receive incoming messages // proc_recv = procs I recv from // length_recv = total size of message each proc sends me // nrecvsize = total size of data I recv *nrecvsize = 0; for (i = 0; i < nrecv; i++) { MPI_Recv(&length_recv[i],1,MPI_INT,MPI_ANY_SOURCE,0,world,status); proc_recv[i] = status->MPI_SOURCE; *nrecvsize += length_recv[i]; } // barrier to insure all MPI_ANY_SOURCE messages are received // else another proc could proceed to irregular_perform() and send to me MPI_Barrier(world); // free work vectors delete [] count; delete [] list; // initialize plan and return it plan->nsend = nsend; plan->nrecv = nrecv; plan->sendmax = sendmax; plan->proc_send = proc_send; plan->length_send = length_send; plan->num_send = num_send; plan->index_send = index_send; plan->offset_send = offset_send; plan->proc_recv = proc_recv; plan->length_recv = length_recv; plan->request = request; plan->status = status; return plan; } /* ---------------------------------------------------------------------- perform irregular communication plan = previouly computed plan via irregular_create() sendbuf = list of atoms to send sizes = # of doubles for each atom recvbuf = received atoms ------------------------------------------------------------------------- */ void Comm::irregular_perform(Plan *plan, double *sendbuf, int *sizes, double *recvbuf) { int i,m,n,offset; // post all receives offset = 0; for (int irecv = 0; irecv < plan->nrecv; irecv++) { MPI_Irecv(&recvbuf[offset],plan->length_recv[irecv],MPI_DOUBLE, plan->proc_recv[irecv],0,world,&plan->request[irecv]); offset += plan->length_recv[irecv]; } // allocate buf for largest send double *buf = (double *) memory->smalloc(plan->sendmax*sizeof(double), "comm::irregular"); // send each message // pack buf with list of datums (datum = one atom) // m = index of datum in sendbuf n = 0; for (int isend = 0; isend < plan->nsend; isend++) { offset = 0; for (i = 0; i < plan->num_send[isend]; i++) { m = plan->index_send[n++]; memcpy(&buf[offset],&sendbuf[plan->offset_send[m]], sizes[m]*sizeof(double)); offset += sizes[m]; } MPI_Send(buf,plan->length_send[isend],MPI_DOUBLE, plan->proc_send[isend],0,world); } // free temporary send buffer memory->sfree(buf); // wait on all incoming messages if (plan->nrecv) MPI_Waitall(plan->nrecv,plan->request,plan->status); } /* ---------------------------------------------------------------------- destroy an irregular communication plan ------------------------------------------------------------------------- */ void Comm::irregular_destroy(Plan *plan) { delete [] plan->proc_send; delete [] plan->length_send; delete [] plan->num_send; delete [] plan->index_send; delete [] plan->offset_send; delete [] plan->proc_recv; delete [] plan->length_recv; delete [] plan->request; delete [] plan->status; memory->sfree(plan); } /* ---------------------------------------------------------------------- determine which proc owns atom with x coord x will be in box (orthogonal) or lamda coords (triclinic) ------------------------------------------------------------------------- */ int Comm::irregular_lookup(double *x) { int loc[3]; if (triclinic == 0) { double *boxlo = domain->boxlo; double *boxhi = domain->boxhi; loc[0] = static_cast (procgrid[0] * (x[0]-boxlo[0]) / (boxhi[0]-boxlo[0])); loc[1] = static_cast (procgrid[1] * (x[1]-boxlo[1]) / (boxhi[1]-boxlo[1])); loc[2] = static_cast (procgrid[2] * (x[2]-boxlo[2]) / (boxhi[2]-boxlo[2])); } else { loc[0] = static_cast (procgrid[0] * x[0]); loc[1] = static_cast (procgrid[1] * x[1]); loc[2] = static_cast (procgrid[2] * x[2]); } if (loc[0] < 0) loc[0] = 0; if (loc[0] >= procgrid[0]) loc[0] = procgrid[0] - 1; if (loc[1] < 0) loc[1] = 0; if (loc[1] >= procgrid[1]) loc[1] = procgrid[1] - 1; if (loc[2] < 0) loc[2] = 0; if (loc[2] >= procgrid[2]) loc[2] = procgrid[2] - 1; return grid2proc[loc[0]][loc[1]][loc[2]]; } /* ---------------------------------------------------------------------- assign nprocs to 3d xprd,yprd,zprd box so as to minimize surface area area = surface area of each of 3 faces of simulation box for triclinic, area = cross product of 2 edge vectors stored in h matrix ------------------------------------------------------------------------- */ void Comm::procs2box() { double area[3]; if (domain->triclinic == 0) { area[0] = domain->xprd * domain->yprd; area[1] = domain->xprd * domain->zprd; area[2] = domain->yprd * domain->zprd; } else { double *h = domain->h; double x,y,z; cross(h[0],0.0,0.0,h[5],h[1],0.0,x,y,z); area[0] = sqrt(x*x + y*y + z*z); cross(h[0],0.0,0.0,h[4],h[3],h[2],x,y,z); area[1] = sqrt(x*x + y*y + z*z); cross(h[5],h[1],0.0,h[4],h[3],h[2],x,y,z); area[2] = sqrt(x*x + y*y + z*z); } double bestsurf = 2.0 * (area[0]+area[1]+area[2]); // loop thru all possible factorizations of nprocs // surf = surface area of a proc sub-domain // for 2d, insure ipz = 1 int ipx,ipy,ipz,nremain; double surf; ipx = 1; while (ipx <= nprocs) { if (nprocs % ipx == 0) { nremain = nprocs/ipx; ipy = 1; while (ipy <= nremain) { if (nremain % ipy == 0) { ipz = nremain/ipy; if (domain->dimension == 3 || ipz == 1) { surf = area[0]/ipx/ipy + area[1]/ipx/ipz + area[2]/ipy/ipz; if (surf < bestsurf) { bestsurf = surf; procgrid[0] = ipx; procgrid[1] = ipy; procgrid[2] = ipz; } } } ipy++; } } ipx++; } } /* ---------------------------------------------------------------------- vector cross product: c = a x b ------------------------------------------------------------------------- */ void Comm::cross(double ax, double ay, double az, double bx, double by, double bz, double &cx, double &cy, double &cz) { cx = ay*bz - az*by; cy = az*bx - ax*bz; cz = ax*by - ay*bx; } /* ---------------------------------------------------------------------- realloc the size of the send buffer as needed with BUFFACTOR & BUFEXTRA if flag = 1, realloc if flag = 0, don't need to realloc with copy, just free/malloc ------------------------------------------------------------------------- */ void Comm::grow_send(int n, int flag) { maxsend = static_cast (BUFFACTOR * n); if (flag) buf_send = (double *) memory->srealloc(buf_send,(maxsend+BUFEXTRA)*sizeof(double), "comm:buf_send"); else { memory->sfree(buf_send); buf_send = (double *) memory->smalloc((maxsend+BUFEXTRA)*sizeof(double), "comm:buf_send"); } } /* ---------------------------------------------------------------------- free/malloc the size of the recv buffer as needed with BUFFACTOR ------------------------------------------------------------------------- */ void Comm::grow_recv(int n) { maxrecv = static_cast (BUFFACTOR * n); memory->sfree(buf_recv); buf_recv = (double *) memory->smalloc(maxrecv*sizeof(double), "comm:buf_recv"); } /* ---------------------------------------------------------------------- realloc the size of the iswap sendlist as needed with BUFFACTOR ------------------------------------------------------------------------- */ void Comm::grow_list(int iswap, int n) { maxsendlist[iswap] = static_cast (BUFFACTOR * n); sendlist[iswap] = (int *) memory->srealloc(sendlist[iswap],maxsendlist[iswap]*sizeof(int), "comm:sendlist[iswap]"); } /* ---------------------------------------------------------------------- realloc the buffers needed for swaps ------------------------------------------------------------------------- */ void Comm::grow_swap(int n) { free_swap(); allocate_swap(n); if (style == MULTI) { free_multi(); allocate_multi(n); } sendlist = (int **) memory->srealloc(sendlist,n*sizeof(int *),"comm:sendlist"); maxsendlist = (int *) memory->srealloc(maxsendlist,n*sizeof(int),"comm:maxsendlist"); for (int i = maxswap; i < n; i++) { maxsendlist[i] = BUFMIN; sendlist[i] = (int *) memory->smalloc(BUFMIN*sizeof(int),"comm:sendlist[i]"); } maxswap = n; } /* ---------------------------------------------------------------------- allocation of swap info ------------------------------------------------------------------------- */ void Comm::allocate_swap(int n) { sendnum = (int *) memory->smalloc(n*sizeof(int),"comm:sendnum"); recvnum = (int *) memory->smalloc(n*sizeof(int),"comm:recvnum"); sendproc = (int *) memory->smalloc(n*sizeof(int),"comm:sendproc"); recvproc = (int *) memory->smalloc(n*sizeof(int),"comm:recvproc"); size_comm_recv = (int *) memory->smalloc(n*sizeof(int),"comm:size"); size_reverse_send = (int *) memory->smalloc(n*sizeof(int),"comm:size"); size_reverse_recv = (int *) memory->smalloc(n*sizeof(int),"comm:size"); slablo = (double *) memory->smalloc(n*sizeof(double),"comm:slablo"); slabhi = (double *) memory->smalloc(n*sizeof(double),"comm:slabhi"); firstrecv = (int *) memory->smalloc(n*sizeof(int),"comm:firstrecv"); pbc_flag = (int *) memory->smalloc(n*sizeof(int),"comm:pbc_flag"); pbc = (int **) memory->create_2d_int_array(n,6,"comm:pbc"); } /* ---------------------------------------------------------------------- allocation of multi-type swap info ------------------------------------------------------------------------- */ void Comm::allocate_multi(int n) { multilo = memory->create_2d_double_array(n,atom->ntypes+1,"comm:multilo"); multihi = memory->create_2d_double_array(n,atom->ntypes+1,"comm:multihi"); } /* ---------------------------------------------------------------------- free memory for swaps ------------------------------------------------------------------------- */ void Comm::free_swap() { memory->sfree(sendnum); memory->sfree(recvnum); memory->sfree(sendproc); memory->sfree(recvproc); memory->sfree(size_comm_recv); memory->sfree(size_reverse_send); memory->sfree(size_reverse_recv); memory->sfree(slablo); memory->sfree(slabhi); memory->sfree(firstrecv); memory->sfree(pbc_flag); memory->destroy_2d_int_array(pbc); } /* ---------------------------------------------------------------------- free memory for multi-type swaps ------------------------------------------------------------------------- */ void Comm::free_multi() { memory->destroy_2d_double_array(multilo); memory->destroy_2d_double_array(multihi); } /* ---------------------------------------------------------------------- set communication style ------------------------------------------------------------------------- */ void Comm::set(int narg, char **arg) { - if (narg != 1) error->all("Illegal communicate command"); + if (narg < 1) error->all("Illegal communicate command"); if (strcmp(arg[0],"single") == 0) style = SINGLE; else if (strcmp(arg[0],"multi") == 0) style = MULTI; else error->all("Illegal communicate command"); + + int iarg = 1; + while (iarg < narg) { + if (strcmp(arg[iarg],"group") == 0) { + if (iarg+2 > narg) error->all("Illegal communicate command"); + igroup = group->find(arg[iarg+1]); + if (igroup < 0) + error->all("Invalid group ID in communicate command"); + iarg += 2; + } else error->all("Illegal communicate command"); + } } /* ---------------------------------------------------------------------- return # of bytes of allocated memory ------------------------------------------------------------------------- */ double Comm::memory_usage() { double bytes = 0.0; for (int i = 0; i < nswap; i++) bytes += maxsendlist[i] * sizeof(int); bytes += maxsend * sizeof(double); bytes += maxrecv * sizeof(double); return bytes; } diff --git a/src/comm.h b/src/comm.h index 09cf8046e..b6b8ef057 100644 --- a/src/comm.h +++ b/src/comm.h @@ -1,120 +1,121 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 COMM_H #define COMM_H #include "pointers.h" namespace LAMMPS_NS { class Comm : protected Pointers { public: int me,nprocs; // proc info int style; // single vs multi-type comm int procgrid[3]; // assigned # of procs in each dim int user_procgrid[3]; // user request for procs in each dim int myloc[3]; // which proc I am in each dim int procneigh[3][2]; // my 6 neighboring procs int nswap; // # of swaps to perform int need[3]; // procs I need atoms from in each dim int maxforward_fix,maxreverse_fix; // comm sizes called from Fix,Pair int maxforward_pair,maxreverse_pair; double cutghost[3]; // cutoffs used for acquiring ghost atoms Comm(class LAMMPS *); ~Comm(); void init(); void set_procs(); // setup 3d grid of procs void setup(); // setup 3d communication pattern void communicate(); // communication of atom coords void reverse_communicate(); // reverse communication of forces void exchange(); // move atoms to new procs void borders(); // setup list of atoms to communicate void comm_pair(class Pair *); // forward comm from a Pair void reverse_comm_pair(class Pair *); // reverse comm from a Pair void comm_fix(class Fix *); // forward comm from a Fix void reverse_comm_fix(class Fix *); // reverse comm from a Fix void comm_compute(class Compute *); // forward comm from a Compute void reverse_comm_compute(class Compute *); // reverse comm from a Compute void irregular(); // irregular communication across all procs void set(int, char **); // set communication style double memory_usage(); private: int triclinic; // 0 if domain is orthog, 1 if triclinic int maxswap; // max # of swaps memory is allocated for int *sendnum,*recvnum; // # of atoms to send/recv in each swap int *sendproc,*recvproc; // proc to send/recv to/from at each swap int *size_comm_recv; // # of values to recv in each forward comm int *size_reverse_send; // # to send in each reverse comm int *size_reverse_recv; // # to recv in each reverse comm double *slablo,*slabhi; // bounds of slab to send at each swap double **multilo,**multihi; // bounds of slabs for multi-type swap double **cutghostmulti; // cutghost on a per-type basis int *pbc_flag; // general flag for sending atoms thru PBC int **pbc; // dimension flags for PBC adjustments int comm_x_only,comm_f_only; // 1 if only exchange x,f in for/rev comm int map_style; // non-0 if global->local mapping is done int ***grid2proc; // which proc owns i,j,k loc in 3d grid + int igroup,groupbit; // only communicate this group int *firstrecv; // where to put 1st recv atom in each swap int **sendlist; // list of atoms to send in each swap int *maxsendlist; // max size of send list for each swap double *buf_send; // send buffer for all comm double *buf_recv; // recv buffer for all comm int maxsend,maxrecv; // current size of send/recv buffer int maxforward,maxreverse; // max # of datums in forward/reverse comm struct Plan { // plan for irregular communication int nsend; // # of messages to send int nrecv; // # of messages to recv int sendmax; // # of doubles in largest send message int *proc_send; // procs to send to int *length_send; // # of doubles to send to each proc int *num_send; // # of datums to send to each proc int *index_send; // list of which datums to send to each proc int *offset_send; // where each datum starts in send buffer int *proc_recv; // procs to recv from int *length_recv; // # of doubles to recv from each proc MPI_Request *request; // MPI requests for posted recvs MPI_Status *status; // MPI statuses for WaitAll }; void procs2box(); // map procs to 3d box void cross(double, double, double, double, double, double, double &, double &, double &); // cross product void grow_send(int,int); // reallocate send buffer void grow_recv(int); // free/allocate recv buffer void grow_list(int, int); // reallocate one sendlist void grow_swap(int); // grow swap and multi arrays void allocate_swap(int); // allocate swap arrays void allocate_multi(int); // allocate multi arrays void free_swap(); // free swap arrays void free_multi(); // free multi arrays struct Plan *irregular_create(int, int *, int *, int *); void irregular_perform(Plan *, double *, int *, double *); void irregular_destroy(Plan *); int irregular_lookup(double *); }; } #endif diff --git a/src/neigh_derive.cpp b/src/neigh_derive.cpp index d2c11f129..cb1a2a670 100644 --- a/src/neigh_derive.cpp +++ b/src/neigh_derive.cpp @@ -1,498 +1,493 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- build half list from full list pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) works if full list is a skip list ------------------------------------------------------------------------- */ void Neighbor::half_from_full_no_newton(NeighList *list) { int i,j,ii,jj,n,jnum; int *neighptr,*jlist; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *ilist_full = list->listfull->ilist; int *numneigh_full = list->listfull->numneigh; int **firstneigh_full = list->listfull->firstneigh; int inum_full = list->listfull->inum; int inum = 0; int npage = 0; int npnt = 0; // loop over atoms in full list for (ii = 0; ii < inum_full; ii++) { if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; // loop over parent full list i = ilist_full[ii]; jlist = firstneigh_full[i]; jnum = numneigh_full[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; if (j > i) neighptr[n++] = j; } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- build half list from full list pair stored once if i,j are both owned and i < j if j is ghost, only store if j coords are "above and to the right" of i works if full list is a skip list ------------------------------------------------------------------------- */ void Neighbor::half_from_full_newton(NeighList *list) { int i,j,ii,jj,n,jnum,joriginal; int *neighptr,*jlist; double xtmp,ytmp,ztmp; double **x = atom->x; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *ilist_full = list->listfull->ilist; int *numneigh_full = list->listfull->numneigh; int **firstneigh_full = list->listfull->firstneigh; int inum_full = list->listfull->inum; int inum = 0; int npage = 0; int npnt = 0; // loop over parent full list for (ii = 0; ii < inum_full; ii++) { if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; i = ilist_full[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over full neighbor list jlist = firstneigh_full[i]; jnum = numneigh_full[i]; for (jj = 0; jj < jnum; jj++) { j = joriginal = jlist[jj]; if (j < nlocal) { if (i > j) continue; } else { if (j >= nall) j %= nall; if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } neighptr[n++] = joriginal; } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip this is for half and full lists ------------------------------------------------------------------------- */ void Neighbor::skip_from(NeighList *list) { int i,j,ii,jj,n,itype,jnum,joriginal; int *neighptr,*jlist; int *type = atom->type; int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *ilist_skip = list->listskip->ilist; int *numneigh_skip = list->listskip->numneigh; int **firstneigh_skip = list->listskip->firstneigh; int inum_skip = list->listskip->inum; int *iskip = list->iskip; int **ijskip = list->ijskip; int inum = 0; int npage = 0; int npnt = 0; // loop over atoms in other list // skip I atom entirely if iskip is set for type[I] // skip I,J pair if ijskip is set for type[I],type[J] for (ii = 0; ii < inum_skip; ii++) { i = ilist_skip[ii]; itype = type[i]; if (iskip[itype]) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; // loop over parent non-skip list jlist = firstneigh_skip[i]; jnum = numneigh_skip[i]; for (jj = 0; jj < jnum; jj++) { j = joriginal = jlist[jj]; if (j >= nall) j %= nall; if (ijskip[itype][type[j]]) continue; neighptr[n++] = joriginal; } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip this is for granular lists with history, copy the history values from parent ------------------------------------------------------------------------- */ void Neighbor::skip_from_granular(NeighList *list) { int i,j,ii,jj,n,nn,itype,jnum,joriginal; int *neighptr,*jlist,*touchptr,*touchptr_skip; double *shearptr,*shearptr_skip; int *type = atom->type; int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *ilist_skip = list->listskip->ilist; int *numneigh_skip = list->listskip->numneigh; int **firstneigh_skip = list->listskip->firstneigh; int **firsttouch_skip = list->listskip->listgranhistory->firstneigh; double **firstshear_skip = list->listskip->listgranhistory->firstdouble; int inum_skip = list->listskip->inum; int *iskip = list->iskip; int **ijskip = list->ijskip; NeighList *listgranhistory = list->listgranhistory; int **firsttouch = listgranhistory->firstneigh; double **firstshear = listgranhistory->firstdouble; int **pages_touch = listgranhistory->pages; double **pages_shear = listgranhistory->dpages; int inum = 0; int npage = 0; int npnt = 0; // loop over atoms in other list // skip I atom entirely if iskip is set for type[I] // skip I,J pair if ijskip is set for type[I],type[J] for (ii = 0; ii < inum_skip; ii++) { i = ilist_skip[ii]; itype = type[i]; if (iskip[itype]) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) { pages = list->add_pages(); pages_touch = listgranhistory->add_pages(); pages_shear = listgranhistory->dpages; } } n = 0; neighptr = &pages[npage][npnt]; nn = 0; touchptr = &pages_touch[npage][npnt]; shearptr = &pages_shear[npage][3*npnt]; // loop over parent non-skip granular list and its history info touchptr_skip = firsttouch_skip[i]; shearptr_skip = firstshear_skip[i]; jlist = firstneigh_skip[i]; jnum = numneigh_skip[i]; for (jj = 0; jj < jnum; jj++) { j = joriginal = jlist[jj]; if (j >= nall) j %= nall; if (ijskip[itype][type[j]]) continue; neighptr[n] = joriginal; touchptr[n++] = touchptr_skip[jj]; shearptr[nn++] = shearptr_skip[3*jj]; shearptr[nn++] = shearptr_skip[3*jj+1]; shearptr[nn++] = shearptr_skip[3*jj+2]; } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; firsttouch[i] = touchptr; firstshear[i] = shearptr; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip this is for respa lists, copy the inner/middle values from parent ------------------------------------------------------------------------- */ void Neighbor::skip_from_respa(NeighList *list) { int i,j,ii,jj,n,itype,jnum,joriginal,n_inner,n_middle; int *neighptr,*jlist,*neighptr_inner,*neighptr_middle; int *type = atom->type; int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *ilist_skip = list->listskip->ilist; int *numneigh_skip = list->listskip->numneigh; int **firstneigh_skip = list->listskip->firstneigh; int inum_skip = list->listskip->inum; int *iskip = list->iskip; int **ijskip = list->ijskip; NeighList *listinner = list->listinner; int *numneigh_inner = listinner->numneigh; int **firstneigh_inner = listinner->firstneigh; int **pages_inner = listinner->pages; int *numneigh_inner_skip = list->listskip->listinner->numneigh; int **firstneigh_inner_skip = list->listskip->listinner->firstneigh; NeighList *listmiddle; int *numneigh_middle,**firstneigh_middle,**pages_middle; int *numneigh_middle_skip,**firstneigh_middle_skip; int respamiddle = list->respamiddle; if (respamiddle) { listmiddle = list->listmiddle; numneigh_middle = listmiddle->numneigh; firstneigh_middle = listmiddle->firstneigh; pages_middle = listmiddle->pages; numneigh_middle_skip = list->listskip->listmiddle->numneigh; firstneigh_middle_skip = list->listskip->listmiddle->firstneigh; } int inum = 0; int npage = 0; int npnt = 0; int npage_inner = 0; int npnt_inner = 0; int npage_middle = 0; int npnt_middle = 0; // loop over atoms in other list // skip I atom entirely if iskip is set for type[I] // skip I,J pair if ijskip is set for type[I],type[J] for (ii = 0; ii < inum_skip; ii++) { i = ilist_skip[ii]; itype = type[i]; if (iskip[itype]) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; if (pgsize - npnt_inner < oneatom) { npnt_inner = 0; npage_inner++; if (npage_inner == listinner->maxpage) pages_inner = listinner->add_pages(); } neighptr_inner = &pages_inner[npage_inner][npnt_inner]; n_inner = 0; if (respamiddle) { if (pgsize - npnt_middle < oneatom) { npnt_middle = 0; npage_middle++; if (npage_middle == listmiddle->maxpage) pages_middle = listmiddle->add_pages(); } neighptr_middle = &pages_middle[npage_middle][npnt_middle]; n_middle = 0; } // loop over parent outer rRESPA list jlist = firstneigh_skip[i]; jnum = numneigh_skip[i]; for (jj = 0; jj < jnum; jj++) { j = joriginal = jlist[jj]; if (j >= nall) j %= nall; if (ijskip[itype][type[j]]) continue; neighptr[n++] = joriginal; } // loop over parent inner rRESPA list jlist = firstneigh_inner_skip[i]; jnum = numneigh_inner_skip[i]; for (jj = 0; jj < jnum; jj++) { j = joriginal = jlist[jj]; if (j >= nall) j %= nall; if (ijskip[itype][type[j]]) continue; neighptr_inner[n_inner++] = joriginal; } // loop over parent middle rRESPA list if (respamiddle) { jlist = firstneigh_middle_skip[i]; jnum = numneigh_middle_skip[i]; for (jj = 0; jj < jnum; jj++) { j = joriginal = jlist[jj]; if (j >= nall) j %= nall; if (ijskip[itype][type[j]]) continue; neighptr_middle[n_middle++] = joriginal; } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; if (npnt_inner >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); if (respamiddle) { firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; if (npnt_middle >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } } list->inum = inum; } /* ---------------------------------------------------------------------- create list which is simply a copy of parent list ------------------------------------------------------------------------- */ void Neighbor::copy_from(NeighList *list) { NeighList *listcopy = list->listcopy; list->inum = listcopy->inum; list->ilist = listcopy->ilist; list->numneigh = listcopy->numneigh; list->firstneigh = listcopy->firstneigh; list->firstdouble = listcopy->firstdouble; list->pages = listcopy->pages; list->dpages = listcopy->dpages; } diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp index 878ab0819..81c6c9586 100644 --- a/src/neigh_full.cpp +++ b/src/neigh_full.cpp @@ -1,284 +1,284 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- N^2 search for all neighbors every neighbor pair appears in list of both atoms i and j ------------------------------------------------------------------------- */ void Neighbor::full_nsq(NeighList *list) { int i,j,n,itype,jtype,which; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over all atoms, owned and ghost // skip i = j for (j = 0; j < nall; j++) { if (i == j) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- binned neighbor list construction for all neighbors every neighbor pair appears in list of both atoms i and j ------------------------------------------------------------------------- */ void Neighbor::full_bin(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &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 surrounding bins in stencil including self // skip i = j ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { if (i == j) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- binned neighbor list construction for all neighbors multi-type stencil is itype dependent and is distance checked every neighbor pair appears in list of both atoms i and j ------------------------------------------------------------------------- */ void Neighbor::full_multi(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *nstencil_multi = list->nstencil_multi; int **stencil_multi = list->stencil_multi; double **distsq_multi = list->distsq_multi; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &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 // skip if i,j neighbor cutoff is less than bin distance // skip i = j ibin = coord2bin(x[i]); s = stencil_multi[itype]; distsq = distsq_multi[itype]; cutsq = cutneighsq[itype]; ns = nstencil_multi[itype]; for (k = 0; k < ns; k++) { for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { jtype = type[j]; if (cutsq[jtype] < distsq[k]) continue; if (i == j) continue; 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } diff --git a/src/neigh_gran.cpp b/src/neigh_gran.cpp index 915caa8e4..8adbc4501 100644 --- a/src/neigh_gran.cpp +++ b/src/neigh_gran.cpp @@ -1,605 +1,603 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "fix_shear_history.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- granular particles N^2 / 2 search for neighbor pairs with partial Newton's 3rd law shear history must be accounted for when a neighbor pair is added pair added to list if atoms i and j are both owned and i < j pair added if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ void Neighbor::granular_nsq_no_newton(NeighList *list) { int i,j,m,n,nn; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr,*touchptr; double *shearptr; NeighList *listgranhistory; int *npartner,**partner; double ***shearpartner; int **firsttouch; double **firstshear; int **pages_touch; double **pages_shear; double **x = atom->x; double *radius = atom->radius; int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; FixShearHistory *fix_history = list->fix_history; if (fix_history) { npartner = fix_history->npartner; partner = fix_history->partner; shearpartner = fix_history->shearpartner; listgranhistory = list->listgranhistory; firsttouch = listgranhistory->firstneigh; firstshear = listgranhistory->firstdouble; pages_touch = listgranhistory->pages; pages_shear = listgranhistory->dpages; } int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) { pages = list->add_pages(); if (fix_history) { pages_touch = listgranhistory->add_pages(); pages_shear = listgranhistory->dpages; } } } n = 0; neighptr = &pages[npage][npnt]; if (fix_history) { nn = 0; touchptr = &pages_touch[npage][npnt]; shearptr = &pages_shear[npage][3*npnt]; } xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; // loop over remaining atoms, owned and ghost for (j = i+1; j < nall; j++) { if (exclude && exclusion(i,j,type[i],type[j],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; radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { neighptr[n] = j; if (fix_history) { if (rsq < radsum*radsum) { for (m = 0; m < npartner[i]; m++) if (partner[i][m] == tag[j]) break; if (m < npartner[i]) { touchptr[n] = 1; shearptr[nn++] = shearpartner[i][m][0]; shearptr[nn++] = shearpartner[i][m][1]; shearptr[nn++] = shearpartner[i][m][2]; } else { touchptr[n] = 0; shearptr[nn++] = 0.0; shearptr[nn++] = 0.0; shearptr[nn++] = 0.0; } } else { touchptr[n] = 0; shearptr[nn++] = 0.0; shearptr[nn++] = 0.0; shearptr[nn++] = 0.0; } } n++; } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; if (fix_history) { firsttouch[i] = touchptr; firstshear[i] = shearptr; } - - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- granular particles N^2 / 2 search for neighbor pairs with full Newton's 3rd law no shear history is allowed for this option pair added to list if atoms i and j are both owned and i < j if j is ghost only me or other proc adds pair decision based on itag,jtag tests ------------------------------------------------------------------------- */ void Neighbor::granular_nsq_newton(NeighList *list) { int i,j,n,itag,jtag; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; double **x = atom->x; double *radius = atom->radius; int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } n = 0; neighptr = &pages[npage][npnt]; itag = tag[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; // loop over remaining atoms, owned and ghost for (j = i+1; j < nall; j++) { if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; } else { if (x[j][2] < ztmp) continue; else if (x[j][2] == ztmp && x[j][1] < ytmp) continue; else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } } if (exclude && exclusion(i,j,type[i],type[j],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; radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) neighptr[n++] = j; } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- granular particles binned neighbor list construction with partial Newton's 3rd law shear history must be accounted for when a neighbor pair is added each owned atom i checks own bin and surrounding bins in non-Newton stencil pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ void Neighbor::granular_bin_no_newton(NeighList *list) { int i,j,k,m,n,nn,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr,*touchptr; double *shearptr; NeighList *listgranhistory; int *npartner,**partner; double ***shearpartner; int **firsttouch; double **firstshear; int **pages_touch; double **pages_shear; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; double *radius = atom->radius; int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; FixShearHistory *fix_history = list->fix_history; if (fix_history) { npartner = fix_history->npartner; partner = fix_history->partner; shearpartner = fix_history->shearpartner; listgranhistory = list->listgranhistory; firsttouch = listgranhistory->firstneigh; firstshear = listgranhistory->firstdouble; pages_touch = listgranhistory->pages; pages_shear = listgranhistory->dpages; } int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) { pages = list->add_pages(); if (fix_history) { pages_touch = listgranhistory->add_pages(); pages_shear = listgranhistory->dpages; } } } n = 0; neighptr = &pages[npage][npnt]; if (fix_history) { nn = 0; touchptr = &pages_touch[npage][npnt]; shearptr = &pages_shear[npage][3*npnt]; } xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; ibin = coord2bin(x[i]); // loop over all atoms in surrounding bins in stencil including self // only store pair if i < j // stores own/own pairs only once // stores own/ghost pairs on both procs for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { if (j <= i) continue; if (exclude && exclusion(i,j,type[i],type[j],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; radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { neighptr[n] = j; if (fix_history) { if (rsq < radsum*radsum) { for (m = 0; m < npartner[i]; m++) if (partner[i][m] == tag[j]) break; if (m < npartner[i]) { touchptr[n] = 1; shearptr[nn++] = shearpartner[i][m][0]; shearptr[nn++] = shearpartner[i][m][1]; shearptr[nn++] = shearpartner[i][m][2]; } else { touchptr[n] = 0; shearptr[nn++] = 0.0; shearptr[nn++] = 0.0; shearptr[nn++] = 0.0; } } else { touchptr[n] = 0; shearptr[nn++] = 0.0; shearptr[nn++] = 0.0; shearptr[nn++] = 0.0; } } n++; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; if (fix_history) { firsttouch[i] = touchptr; firstshear[i] = shearptr; } - - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- granular particles binned neighbor list construction with full Newton's 3rd law no shear history is allowed for this option each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::granular_bin_newton(NeighList *list) { int i,j,k,n,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } n = 0; neighptr = &pages[npage][npnt]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i for (j = bins[i]; j >= 0; j = bins[j]) { if (j >= nlocal) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; if (exclude && exclusion(i,j,type[i],type[j],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; radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) neighptr[n++] = j; } // loop over all atoms in other bins in stencil, store every pair ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { if (exclude && exclusion(i,j,type[i],type[j],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; radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) neighptr[n++] = j; } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- granular particles binned neighbor list construction with Newton's 3rd law for triclinic no shear history is allowed for this option each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::granular_bin_newton_tri(NeighList *list) { int i,j,k,n,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } n = 0; neighptr = &pages[npage][npnt]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; // loop over all atoms in bins in stencil // pairs for atoms j "below" i are excluded // below = lower z or (equal z and lower y) or (equal zy and <= x) // this excludes self-self interaction ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue; if (exclude && exclusion(i,j,type[i],type[j],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; radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) neighptr[n++] = j; } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } diff --git a/src/neigh_half_bin.cpp b/src/neigh_half_bin.cpp index c0e5ac721..3f19e4cd6 100644 --- a/src/neigh_half_bin.cpp +++ b/src/neigh_half_bin.cpp @@ -1,321 +1,321 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- binned neighbor list construction with partial Newton's 3rd law each owned atom i checks own bin and other bins in stencil pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ void Neighbor::half_bin_no_newton(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &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 // only store pair if i < j // stores own/own pairs only once // stores own/ghost pairs on both procs 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- binned neighbor list construction with full Newton's 3rd law each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::half_bin_newton(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i for (j = bins[i]; j >= 0; j = bins[j]) { if (j >= nlocal) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } // loop over all atoms in other bins in stencil, store every pair ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- binned neighbor list construction with Newton's 3rd law for triclinic each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::half_bin_newton_tri(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &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 bins in stencil // pairs for atoms j "below" i are excluded // below = lower z or (equal z and lower y) or (equal zy and <= x) // this excludes self-self interaction ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } diff --git a/src/neigh_half_multi.cpp b/src/neigh_half_multi.cpp index 8ecad243a..cda4965eb 100644 --- a/src/neigh_half_multi.cpp +++ b/src/neigh_half_multi.cpp @@ -1,347 +1,347 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- binned neighbor list construction with partial Newton's 3rd law each owned atom i checks own bin and other bins in stencil multi-type stencil is itype dependent and is distance checked pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ void Neighbor::half_multi_no_newton(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *nstencil_multi = list->nstencil_multi; int **stencil_multi = list->stencil_multi; double **distsq_multi = list->distsq_multi; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &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 // only store pair if i < j // skip if i,j neighbor cutoff is less than bin distance // stores own/own pairs only once // stores own/ghost pairs on both procs ibin = coord2bin(x[i]); s = stencil_multi[itype]; distsq = distsq_multi[itype]; cutsq = cutneighsq[itype]; ns = nstencil_multi[itype]; for (k = 0; k < ns; k++) { for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { if (j <= i) continue; jtype = type[j]; if (cutsq[jtype] < distsq[k]) continue; 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- binned neighbor list construction with full Newton's 3rd law each owned atom i checks its own bin and other bins in Newton stencil multi-type stencil is itype dependent and is distance checked every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::half_multi_newton(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *nstencil_multi = list->nstencil_multi; int **stencil_multi = list->stencil_multi; double **distsq_multi = list->distsq_multi; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i for (j = bins[i]; j >= 0; j = bins[j]) { if (j >= nlocal) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } // loop over all atoms in other bins in stencil, store every pair // skip if i,j neighbor cutoff is less than bin distance ibin = coord2bin(x[i]); s = stencil_multi[itype]; distsq = distsq_multi[itype]; cutsq = cutneighsq[itype]; ns = nstencil_multi[itype]; for (k = 0; k < ns; k++) { for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { jtype = type[j]; if (cutsq[jtype] < distsq[k]) continue; 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- binned neighbor list construction with Newton's 3rd law for triclinic each owned atom i checks its own bin and other bins in triclinic stencil multi-type stencil is itype dependent and is distance checked every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::half_multi_newton_tri(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int *nstencil_multi = list->nstencil_multi; int **stencil_multi = list->stencil_multi; double **distsq_multi = list->distsq_multi; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &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 bins, including self, in stencil // skip if i,j neighbor cutoff is less than bin distance // bins below self are excluded from stencil // pairs for atoms j below i are excluded ibin = coord2bin(x[i]); s = stencil_multi[itype]; distsq = distsq_multi[itype]; cutsq = cutneighsq[itype]; ns = nstencil_multi[itype]; for (k = 0; k < ns; k++) { for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { jtype = type[j]; if (cutsq[jtype] < distsq[k]) continue; if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue; 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp index ea7ae68e0..aaee116dc 100644 --- a/src/neigh_half_nsq.cpp +++ b/src/neigh_half_nsq.cpp @@ -1,188 +1,188 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- N^2 / 2 search for neighbor pairs with partial Newton's 3rd law pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ void Neighbor::half_nsq_no_newton(NeighList *list) { int i,j,n,itype,jtype,which; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &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 for (j = i+1; j < nall; j++) { 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } /* ---------------------------------------------------------------------- N^2 / 2 search for neighbor pairs with full Newton's 3rd law every pair stored exactly once by some processor decision on ghost atoms based on itag,jtag tests ------------------------------------------------------------------------- */ void Neighbor::half_nsq_newton(NeighList *list) { int i,j,n,itype,jtype,itag,jtag,which; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; double **x = atom->x; int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int inum = 0; int npage = 0; int npnt = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; itag = tag[i]; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over remaining atoms, owned and ghost // itag = jtag is possible for long cutoffs that include images of self for (j = i+1; j < nall; j++) { if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; } else { if (x[j][2] < ztmp) continue; else if (x[j][2] == ztmp && x[j][1] < ytmp) continue; else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } list->inum = inum; } diff --git a/src/neigh_respa.cpp b/src/neigh_respa.cpp index 38bb3232e..125bfecc9 100644 --- a/src/neigh_respa.cpp +++ b/src/neigh_respa.cpp @@ -1,819 +1,819 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- multiple respa lists N^2 / 2 search for neighbor pairs with partial Newton's 3rd law pair added to list if atoms i and j are both owned and i < j pair added if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ void Neighbor::respa_nsq_no_newton(NeighList *list) { int i,j,n,itype,jtype,which,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; NeighList *listinner = list->listinner; int *numneigh_inner = listinner->numneigh; int **firstneigh_inner = listinner->firstneigh; int **pages_inner = listinner->pages; NeighList *listmiddle; int *numneigh_middle,**firstneigh_middle,**pages_middle; int respamiddle = list->respamiddle; if (respamiddle) { listmiddle = list->listmiddle; numneigh_middle = listmiddle->numneigh; firstneigh_middle = listmiddle->firstneigh; pages_middle = listmiddle->pages; } int inum = 0; int npage = 0; int npnt = 0; int npage_inner = 0; int npnt_inner = 0; int npage_middle = 0; int npnt_middle = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; if (pgsize - npnt_inner < oneatom) { npnt_inner = 0; npage_inner++; if (npage_inner == listinner->maxpage) pages_inner = listinner->add_pages(); } neighptr_inner = &pages_inner[npage_inner][npnt_inner]; n_inner = 0; if (respamiddle) { if (pgsize - npnt_middle < oneatom) { npnt_middle = 0; npage_middle++; if (npage_middle == listmiddle->maxpage) pages_middle = listmiddle->add_pages(); } neighptr_middle = &pages_middle[npage_middle][npnt_middle]; n_middle = 0; } itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over remaining atoms, owned and ghost for (j = i+1; j < nall; j++) { 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; if (npnt_inner >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); if (respamiddle) { firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; if (npnt_middle >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } } list->inum = inum; } /* ---------------------------------------------------------------------- multiple respa lists N^2 / 2 search for neighbor pairs with full Newton's 3rd law pair added to list if atoms i and j are both owned and i < j if j is ghost only me or other proc adds pair decision based on itag,jtag tests ------------------------------------------------------------------------- */ void Neighbor::respa_nsq_newton(NeighList *list) { int i,j,n,itype,jtype,itag,jtag,which,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; double **x = atom->x; int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; NeighList *listinner = list->listinner; int *numneigh_inner = listinner->numneigh; int **firstneigh_inner = listinner->firstneigh; int **pages_inner = listinner->pages; NeighList *listmiddle; int *numneigh_middle,**firstneigh_middle,**pages_middle; int respamiddle = list->respamiddle; if (respamiddle) { listmiddle = list->listmiddle; numneigh_middle = listmiddle->numneigh; firstneigh_middle = listmiddle->firstneigh; pages_middle = listmiddle->pages; } int inum = 0; int npage = 0; int npnt = 0; int npage_inner = 0; int npnt_inner = 0; int npage_middle = 0; int npnt_middle = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; if (pgsize - npnt_inner < oneatom) { npnt_inner = 0; npage_inner++; if (npage_inner == listinner->maxpage) pages_inner = listinner->add_pages(); } neighptr_inner = &pages_inner[npage_inner][npnt_inner]; n_inner = 0; if (respamiddle) { if (pgsize - npnt_middle < oneatom) { npnt_middle = 0; npage_middle++; if (npage_middle == listmiddle->maxpage) pages_middle = listmiddle->add_pages(); } neighptr_middle = &pages_middle[npage_middle][npnt_middle]; n_middle = 0; } itag = tag[i]; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over remaining atoms, owned and ghost for (j = i+1; j < nall; j++) { if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; } else { if (x[j][2] < ztmp) continue; else if (x[j][2] == ztmp && x[j][1] < ytmp) continue; else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; if (npnt_inner >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); if (respamiddle) { firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; if (npnt_middle >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } } list->inum = inum; } /* ---------------------------------------------------------------------- multiple respa lists binned neighbor list construction with partial Newton's 3rd law each owned atom i checks own bin and surrounding bins in non-Newton stencil pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ void Neighbor::respa_bin_no_newton(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; NeighList *listinner = list->listinner; int *numneigh_inner = listinner->numneigh; int **firstneigh_inner = listinner->firstneigh; int **pages_inner = listinner->pages; NeighList *listmiddle; int *numneigh_middle,**firstneigh_middle,**pages_middle; int respamiddle = list->respamiddle; if (respamiddle) { listmiddle = list->listmiddle; numneigh_middle = listmiddle->numneigh; firstneigh_middle = listmiddle->firstneigh; pages_middle = listmiddle->pages; } int inum = 0; int npage = 0; int npnt = 0; int npage_inner = 0; int npnt_inner = 0; int npage_middle = 0; int npnt_middle = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; if (pgsize - npnt_inner < oneatom) { npnt_inner = 0; npage_inner++; if (npage_inner == listinner->maxpage) pages_inner = listinner->add_pages(); } neighptr_inner = &pages_inner[npage_inner][npnt_inner]; n_inner = 0; if (respamiddle) { if (pgsize - npnt_middle < oneatom) { npnt_middle = 0; npage_middle++; if (npage_middle == listmiddle->maxpage) pages_middle = listmiddle->add_pages(); } neighptr_middle = &pages_middle[npage_middle][npnt_middle]; n_middle = 0; } itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; ibin = coord2bin(x[i]); // loop over all atoms in surrounding bins in stencil including self // only store pair if i < j // stores own/own pairs only once // stores own/ghost pairs on both procs 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; if (npnt_inner >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); if (respamiddle) { firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; if (npnt_middle >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } } list->inum = inum; } /* ---------------------------------------------------------------------- multiple respa lists binned neighbor list construction with full Newton's 3rd law each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::respa_bin_newton(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; NeighList *listinner = list->listinner; int *numneigh_inner = listinner->numneigh; int **firstneigh_inner = listinner->firstneigh; int **pages_inner = listinner->pages; NeighList *listmiddle; int *numneigh_middle,**firstneigh_middle,**pages_middle; int respamiddle = list->respamiddle; if (respamiddle) { listmiddle = list->listmiddle; numneigh_middle = listmiddle->numneigh; firstneigh_middle = listmiddle->firstneigh; pages_middle = listmiddle->pages; } int inum = 0; int npage = 0; int npnt = 0; int npage_inner = 0; int npnt_inner = 0; int npage_middle = 0; int npnt_middle = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; if (pgsize - npnt_inner < oneatom) { npnt_inner = 0; npage_inner++; if (npage_inner == listinner->maxpage) pages_inner = listinner->add_pages(); } neighptr_inner = &pages_inner[npage_inner][npnt_inner]; n_inner = 0; if (respamiddle) { if (pgsize - npnt_middle < oneatom) { npnt_middle = 0; npage_middle++; if (npage_middle == listmiddle->maxpage) pages_middle = listmiddle->add_pages(); } neighptr_middle = &pages_middle[npage_middle][npnt_middle]; n_middle = 0; } itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i for (j = bins[i]; j >= 0; j = bins[j]) { if (j >= nlocal) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } } // loop over all atoms in other bins in stencil, store every pair ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; if (npnt_inner >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); if (respamiddle) { firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; if (npnt_middle >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } } list->inum = inum; } /* ---------------------------------------------------------------------- multiple respa lists binned neighbor list construction with Newton's 3rd law for triclinic each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::respa_bin_newton_tri(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; // bin local & ghost atoms bin_atoms(); // loop over each atom, storing neighbors double **x = atom->x; int *type = atom->type; int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; int **pages = list->pages; int nstencil = list->nstencil; int *stencil = list->stencil; NeighList *listinner = list->listinner; int *numneigh_inner = listinner->numneigh; int **firstneigh_inner = listinner->firstneigh; int **pages_inner = listinner->pages; NeighList *listmiddle; int *numneigh_middle,**firstneigh_middle,**pages_middle; int respamiddle = list->respamiddle; if (respamiddle) { listmiddle = list->listmiddle; numneigh_middle = listmiddle->numneigh; firstneigh_middle = listmiddle->firstneigh; pages_middle = listmiddle->pages; } int inum = 0; int npage = 0; int npnt = 0; int npage_inner = 0; int npnt_inner = 0; int npage_middle = 0; int npnt_middle = 0; for (i = 0; i < nlocal; i++) { + if (include_group && !(mask[i] & include_groupbit)) continue; if (pgsize - npnt < oneatom) { npnt = 0; npage++; if (npage == list->maxpage) pages = list->add_pages(); } neighptr = &pages[npage][npnt]; n = 0; if (pgsize - npnt_inner < oneatom) { npnt_inner = 0; npage_inner++; if (npage_inner == listinner->maxpage) pages_inner = listinner->add_pages(); } neighptr_inner = &pages_inner[npage_inner][npnt_inner]; n_inner = 0; if (respamiddle) { if (pgsize - npnt_middle < oneatom) { npnt_middle = 0; npage_middle++; if (npage_middle == listmiddle->maxpage) pages_middle = listmiddle->add_pages(); } neighptr_middle = &pages_middle[npage_middle][npnt_middle]; n_middle = 0; } itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; // loop over all atoms in bins in stencil // pairs for atoms j "below" i are excluded // below = lower z or (equal z and lower y) or (equal zy and <= x) // this excludes self-self interaction ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) 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(i,j); else which = 0; if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } } } - ilist[inum] = i; + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; - inum++; npnt += n; if (npnt >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; if (npnt_inner >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); if (respamiddle) { firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; if (npnt_middle >= pgsize) error->one("Neighbor list overflow, boost neigh_modify one or page"); } } list->inum = inum; } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index b9af0647b..0c0b6b2b8 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -1,1517 +1,1535 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author (triclinic and multi-neigh) : Pieter in 't Veld (SNL) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "stdlib.h" #include "string.h" #include "limits.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "atom.h" #include "atom_vec.h" #include "comm.h" #include "force.h" #include "pair.h" #include "domain.h" #include "group.h" #include "modify.h" #include "fix.h" #include "compute.h" #include "update.h" #include "respa.h" #include "output.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define RQDELTA 1 #define EXDELTA 1 #define LB_FACTOR 1.5 #define SMALL 1.0e-6 #define BIG 1.0e20 #define CUT2BIN_RATIO 100 #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) enum{NSQ,BIN,MULTI}; // also in neigh_list.cpp //#define NEIGH_LIST_DEBUG 1 /* ---------------------------------------------------------------------- */ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); style = BIN; every = 1; delay = 10; dist_check = 1; pgsize = 100000; oneatom = 2000; binsizeflag = 0; cutneighsq = NULL; cuttype = NULL; cuttypesq = NULL; fixchecklist = NULL; // coords at last neighboring maxhold = 0; xhold = NULL; // binning maxhead = 0; binhead = NULL; maxbin = 0; bins = NULL; // pair exclusion list info + include_group = 0; + nex_type = maxex_type = 0; ex1_type = ex2_type = NULL; ex_type = NULL; nex_group = maxex_group = 0; ex1_group = ex2_group = ex1_bit = ex2_bit = NULL; nex_mol = maxex_mol = 0; ex_mol_group = ex_mol_bit = NULL; // pair lists maxlocal = 0; nlist = 0; lists = NULL; pair_build = NULL; stencil_create = NULL; blist = glist = slist = NULL; nrequest = maxrequest = 0; requests = NULL; old_style = BIN; old_nrequest = 0; old_requests = NULL; // bond lists maxbond = 0; bondlist = NULL; maxangle = 0; anglelist = NULL; maxdihedral = 0; dihedrallist = NULL; maximproper = 0; improperlist = NULL; } /* ---------------------------------------------------------------------- */ Neighbor::~Neighbor() { memory->destroy_2d_double_array(cutneighsq); delete [] cuttype; delete [] cuttypesq; delete [] fixchecklist; memory->destroy_2d_double_array(xhold); memory->sfree(binhead); memory->sfree(bins); memory->sfree(ex1_type); memory->sfree(ex2_type); memory->destroy_2d_int_array(ex_type); memory->sfree(ex1_group); memory->sfree(ex2_group); delete [] ex1_bit; delete [] ex2_bit; memory->sfree(ex_mol_group); delete [] ex_mol_bit; for (int i = 0; i < nlist; i++) delete lists[i]; delete [] lists; delete [] pair_build; delete [] stencil_create; delete [] blist; delete [] glist; delete [] slist; for (int i = 0; i < nrequest; i++) delete requests[i]; memory->sfree(requests); for (int i = 0; i < old_nrequest; i++) delete old_requests[i]; memory->sfree(old_requests); memory->destroy_2d_int_array(bondlist); memory->destroy_2d_int_array(anglelist); memory->destroy_2d_int_array(dihedrallist); memory->destroy_2d_int_array(improperlist); } /* ---------------------------------------------------------------------- */ void Neighbor::init() { int i,j,m,n; ncalls = ndanger = 0; dimension = domain->dimension; triclinic = domain->triclinic; newton_pair = force->newton_pair; // error check if (delay > 0 && (delay % every) != 0) error->all("Neighbor delay must be 0 or multiple of every setting"); if (pgsize < 10*oneatom) error->all("Neighbor page size must be >= 10x the one atom setting"); // ------------------------------------------------------------------ // settings // set neighbor cutoffs (force cutoff + skin) // trigger determines when atoms migrate and neighbor lists are rebuilt // cutneigh and cutneighsq determine what pairs go into neighbor list // set to 0 if cutforce = 0 // cutneighmin/max used for neighbor bin sizes // cutghost determines comm distance = max of cutneigh & skin // skin is only > cutneighmax when cutneighmax = 0.0 due to pair = NULL // in this case, are running a bond-only simulation // still need ghost atom communication for bonds and atom exchange // triggered by atoms moving 1/2 skin distance, but no neigh lists formed // user should set skin big enough to find all bonds triggersq = 0.25*skin*skin; n = atom->ntypes; if (cutneighsq == NULL) { cutneighsq = memory->create_2d_double_array(n+1,n+1,"neigh:cutneighsq"); cuttype = new double[n+1]; cuttypesq = new double[n+1]; } double cutoff,delta,cut; cutneighmin = BIG; cutneighmax = 0.0; for (i = 1; i <= n; i++) { cuttype[i] = cuttypesq[i] = 0.0; for (j = 1; j <= n; j++) { if (force->pair) cutoff = sqrt(force->pair->cutsq[i][j]); else cutoff = 0.0; if (cutoff > 0.0) delta = skin; else delta = 0.0; cut = cutoff + delta; cutneighsq[i][j] = cut*cut; cuttype[i] = MAX(cuttype[i],cut); cuttypesq[i] = MAX(cuttypesq[i],cut*cut); cutneighmin = MIN(cutneighmin,cut); cutneighmax = MAX(cutneighmax,cut); } } cutghost = MAX(cutneighmax,skin); cutneighmaxsq = cutneighmax * cutneighmax; // check other classes that can induce reneighboring in decide() restart_check = 0; if (output->restart_every) restart_check = 1; delete [] fixchecklist; fixchecklist = NULL; fixchecklist = new int[modify->nfix]; fix_check = 0; for (i = 0; i < modify->nfix; i++) if (modify->fix[i]->force_reneighbor) fixchecklist[fix_check++] = i; must_check = 0; if (restart_check || fix_check) must_check = 1; // set special_flag for 1-2, 1-3, 1-4 neighbors // flag[0] is not used, flag[1] = 1-2, flag[2] = 1-3, flag[3] = 1-4 // flag = 0 if both LJ/Coulomb special values are 0.0 // flag = 1 if both LJ/Coulomb special values are 1.0 // flag = 2 otherwise or if KSpace solver is enabled // pairwise portion of KSpace solver uses all 1-2,1-3,1-4 neighbors if (force->special_lj[1] == 0.0 && force->special_coul[1] == 0.0) special_flag[1] = 0; else if (force->special_lj[1] == 1.0 && force->special_coul[1] == 1.0) special_flag[1] = 1; else special_flag[1] = 2; if (force->special_lj[2] == 0.0 && force->special_coul[2] == 0.0) special_flag[2] = 0; else if (force->special_lj[2] == 1.0 && force->special_coul[2] == 1.0) special_flag[2] = 1; else special_flag[2] = 2; if (force->special_lj[3] == 0.0 && force->special_coul[3] == 0.0) special_flag[3] = 0; else if (force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0) special_flag[3] = 1; else special_flag[3] = 2; if (force->kspace) special_flag[1] = special_flag[2] = special_flag[3] = 2; // rRESPA cutoffs int respa = 0; if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) { if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; } if (respa) { double *cut_respa = ((Respa *) update->integrate)->cutoff; cut_inner_sq = (cut_respa[1] + skin) * (cut_respa[1] + skin); cut_middle_sq = (cut_respa[3] + skin) * (cut_respa[3] + skin); cut_middle_inside_sq = (cut_respa[0] - skin) * (cut_respa[0] - skin); } // ------------------------------------------------------------------ // xhold, bins, exclusion lists // free xhold and bins if not needed for this run if (dist_check == 0) { memory->destroy_2d_double_array(xhold); maxhold = 0; xhold = NULL; } if (style == NSQ) { memory->sfree(bins); memory->sfree(binhead); maxbin = maxhead = 0; binhead = NULL; bins = NULL; } // 1st time allocation of xhold and bins if (dist_check) { if (maxhold == 0) { maxhold = atom->nmax; xhold = memory->create_2d_double_array(maxhold,3,"neigh:xhold"); } } if (style != NSQ) { if (maxbin == 0) { maxbin = atom->nmax; bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins"); } } // exclusion lists for type, group, molecule settings from neigh_modify n = atom->ntypes; + include_groupbit = group->bitmask[include_group]; + if (nex_type == 0 && nex_group == 0 && nex_mol == 0) exclude = 0; else exclude = 1; if (nex_type) { memory->destroy_2d_int_array(ex_type); ex_type = (int **) memory->create_2d_int_array(n+1,n+1,"neigh:ex_type"); for (i = 1; i <= n; i++) for (j = 1; j <= n; j++) ex_type[i][j] = 0; for (i = 0; i < nex_type; i++) { if (ex1_type[i] <= 0 || ex1_type[i] > n || ex2_type[i] <= 0 || ex2_type[i] > n) error->all("Invalid atom type in neighbor exclusion list"); ex_type[ex1_type[i]][ex2_type[i]] = 1; ex_type[ex2_type[i]][ex1_type[i]] = 1; } } if (nex_group) { delete [] ex1_bit; delete [] ex2_bit; ex1_bit = new int[nex_group]; ex2_bit = new int[nex_group]; for (i = 0; i < nex_group; i++) { ex1_bit[i] = group->bitmask[ex1_group[i]]; ex2_bit[i] = group->bitmask[ex2_group[i]]; } } if (nex_mol) { delete [] ex_mol_bit; ex_mol_bit = new int[nex_mol]; for (i = 0; i < nex_mol; i++) ex_mol_bit[i] = group->bitmask[ex_mol_group[i]]; } // ------------------------------------------------------------------ // pairwise lists // test if pairwise lists need to be re-created // no need to re-create if: // neigh style has not changed and current requests = old requests int same = 1; if (style != old_style) same = 0; if (nrequest != old_nrequest) same = 0; else for (i = 0; i < nrequest; i++) if (requests[i]->identical(old_requests[i]) == 0) same = 0; #ifdef NEIGH_LIST_DEBUG if (comm->me == 0) printf("SAME flag %d\n",same); #endif // if old and new are not the same, create new pairwise lists if (!same) { // delete old lists and create new ones for (i = 0; i < nlist; i++) delete lists[i]; delete [] lists; delete [] pair_build; delete [] stencil_create; nlist = nrequest; lists = new NeighList*[nlist]; pair_build = new PairPtr[nlist]; stencil_create = new StencilPtr[nlist]; // create individual lists, one per request // copy dnum setting from request to list // pass list ptr back to requestor (except for Command class) for (i = 0; i < nlist; i++) { lists[i] = new NeighList(lmp,pgsize); lists[i]->index = i; lists[i]->dnum = requests[i]->dnum; if (requests[i]->pair) { Pair *pair = (Pair *) requests[i]->requestor; pair->init_list(requests[i]->id,lists[i]); } else if (requests[i]->fix) { Fix *fix = (Fix *) requests[i]->requestor; fix->init_list(requests[i]->id,lists[i]); } else if (requests[i]->compute) { Compute *compute = (Compute *) requests[i]->requestor; compute->init_list(requests[i]->id,lists[i]); } } // detect lists that are connected to other lists // if-then-else sequence is important // since don't want to re-process skip or copy lists further down // copy: point this list at request->otherlist, could be a skip list // skip: point this list at request->otherlist, copy skip info from request // half_from_full: point this list at preceeding full list // granhistory: set preceeding list's listgranhistory to this list // also set precedding list's ptr to FixShearHistory // respaouter: point this list at preceeding 1/2 inner/middle lists // pair and half: if there is a full non-occasional non-skip list // change this list to half_from_full and point at the full list // parent could be copy list or pair or fix // fix/compute requests: // kind of request = half or full, occasional or not doesn't matter // if request = half and non-skip pair half/respaouter exists, // become copy of that list // if request = full and non-skip pair full exists, // become copy of that list // if request = half and non-skip pair full exists, // become half_from_full of that list // if no matches, do nothing, fix/compute list will be built directly // ok if parent is copy list for (i = 0; i < nlist; i++) { if (requests[i]->copy) lists[i]->listcopy = lists[requests[i]->otherlist]; else if (requests[i]->skip) { lists[i]->listskip = lists[requests[i]->otherlist]; lists[i]->copy_skip_info(requests[i]->iskip,requests[i]->ijskip); } else if (requests[i]->half_from_full) lists[i]->listfull = lists[i-1]; else if (requests[i]->granhistory) { lists[i-1]->listgranhistory = lists[i]; for (int ifix = 0; ifix < modify->nfix; ifix++) if (strcmp(modify->fix[ifix]->style,"SHEAR_HISTORY") == 0) lists[i-1]->fix_history = (FixShearHistory *) modify->fix[ifix]; } else if (requests[i]->respaouter) { if (requests[i-1]->respainner) { lists[i]->respamiddle = 0; lists[i]->listinner = lists[i-1]; } else { lists[i]->respamiddle = 1; lists[i]->listmiddle = lists[i-1]; lists[i]->listinner = lists[i-2]; } } else if (requests[i]->pair && requests[i]->half) { for (j = 0; j < nlist; j++) if (requests[j]->full && requests[j]->occasional == 0 && requests[j]->skip == 0) break; if (j < nlist) { requests[i]->half = 0; requests[i]->half_from_full = 1; lists[i]->listfull = lists[j]; } } else if (requests[i]->fix || requests[i]->compute) { for (j = 0; j < nlist; j++) { if (requests[i]->half && requests[j]->pair && requests[j]->skip == 0 && requests[j]->half) break; if (requests[i]->full && requests[j]->pair && requests[j]->skip == 0 && requests[j]->full) break; if (requests[i]->half && requests[j]->pair && requests[j]->skip == 0 && requests[j]->respaouter) break; } if (j < nlist) { requests[i]->copy = 1; lists[i]->listcopy = lists[j]; } else { for (j = 0; j < nlist; j++) { if (requests[i]->half && requests[j]->pair && requests[j]->skip == 0 && requests[j]->full) break; } if (j < nlist) { requests[i]->half = 0; requests[i]->half_from_full = 1; lists[i]->listfull = lists[j]; } } } } // set ptrs to pair_build and stencil_create functions for each list // ptrs set to NULL if not set explicitly for (i = 0; i < nlist; i++) { choose_build(i,requests[i]); if (style != NSQ) choose_stencil(i,requests[i]); else stencil_create[i] = NULL; } // set each list's build/grow/stencil flags based on neigh request // buildflag = 1 if its pair_build() invoked every reneighbor // growflag = 1 if it stores atom-based arrays and pages // stencilflag = 1 if it stores stencil arrays for (i = 0; i < nlist; i++) { lists[i]->buildflag = 1; if (pair_build[i] == NULL) lists[i]->buildflag = 0; if (requests[i]->occasional) lists[i]->buildflag = 0; lists[i]->growflag = 1; if (requests[i]->copy) lists[i]->growflag = 0; lists[i]->stencilflag = 1; if (style == NSQ) lists[i]->stencilflag = 0; if (stencil_create[i] == NULL) lists[i]->stencilflag = 0; } #ifdef NEIGH_LIST_DEBUG for (i = 0; i < nlist; i++) lists[i]->print_attributes(); #endif // allocate atom arrays and 1st pages of lists that store them maxlocal = atom->nmax; for (i = 0; i < nlist; i++) if (lists[i]->growflag) { lists[i]->grow(maxlocal); lists[i]->add_pages(); } // setup 3 vectors of pairwise neighbor lists // blist = lists whose pair_build() is invoked every reneighbor // glist = lists who store atom arrays which are used every reneighbor // slist = lists who store stencil arrays which are used every reneighbor // blist and glist vectors are used by neighbor::build() // slist vector is used by neighbor::setup_bins() nblist = nglist = nslist = 0; delete [] blist; delete [] glist; delete [] slist; blist = new int[nlist]; glist = new int[nlist]; slist = new int[nlist]; for (i = 0; i < nlist; i++) { if (lists[i]->buildflag) blist[nblist++] = i; if (lists[i]->growflag && requests[i]->occasional == 0) glist[nglist++] = i; if (lists[i]->stencilflag && requests[i]->occasional == 0) slist[nslist++] = i; } #ifdef NEIGH_LIST_DEBUG print_lists_of_lists(); #endif // reorder build vector if necessary // relevant for lists that copy/skip/half-full from parent // the derived list must appear in blist after the parent list // no occasional lists are in build vector // swap two lists within blist when dependency is mis-ordered // done when entire pass thru blist results in no swaps int done = 0; while (!done) { done = 1; for (i = 0; i < nblist; i++) { NeighList *ptr = NULL; if (lists[blist[i]]->listfull) ptr = lists[blist[i]]->listfull; if (lists[blist[i]]->listcopy) ptr = lists[blist[i]]->listcopy; if (lists[blist[i]]->listskip) ptr = lists[blist[i]]->listskip; if (ptr == NULL) continue; for (m = 0; m < nlist; m++) if (ptr == lists[m]) break; for (j = 0; j < nblist; j++) if (m == blist[j]) break; if (j < i) continue; int tmp = blist[i]; blist[i] = blist[j]; blist[j] = tmp; done = 0; break; } } #ifdef NEIGH_LIST_DEBUG print_lists_of_lists(); #endif } // delete old requests // copy current requests and style to old for next run for (i = 0; i < old_nrequest; i++) delete old_requests[i]; memory->sfree(old_requests); old_nrequest = nrequest; old_requests = requests; nrequest = maxrequest = 0; requests = NULL; old_style = style; // ------------------------------------------------------------------ // topology lists // 1st time allocation of topology lists if (atom->molecular && atom->nbonds && maxbond == 0) { if (nprocs == 1) maxbond = atom->nbonds; else maxbond = static_cast (LB_FACTOR * atom->nbonds / nprocs); bondlist = memory->create_2d_int_array(maxbond,3,"neigh:bondlist"); } if (atom->molecular && atom->nangles && maxangle == 0) { if (nprocs == 1) maxangle = atom->nangles; else maxangle = static_cast (LB_FACTOR * atom->nangles / nprocs); anglelist = memory->create_2d_int_array(maxangle,4,"neigh:anglelist"); } if (atom->molecular && atom->ndihedrals && maxdihedral == 0) { if (nprocs == 1) maxdihedral = atom->ndihedrals; else maxdihedral = static_cast (LB_FACTOR * atom->ndihedrals / nprocs); dihedrallist = memory->create_2d_int_array(maxdihedral,5,"neigh:dihedrallist"); } if (atom->molecular && atom->nimpropers && maximproper == 0) { if (nprocs == 1) maximproper = atom->nimpropers; else maximproper = static_cast (LB_FACTOR * atom->nimpropers / nprocs); improperlist = memory->create_2d_int_array(maximproper,5,"neigh:improperlist"); } // set flags that determine which topology neighboring routines to use // SHAKE sets bonds and angles negative // bond_quartic sets bonds to 0 // delete_bonds sets all interactions negative int bond_off = 0; int angle_off = 0; for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"shake") == 0) bond_off = angle_off = 1; if (force->bond && force->bond_match("quartic")) bond_off = 1; if (atom->avec->bonds_allow) { for (i = 0; i < atom->nlocal; i++) { if (bond_off) break; for (m = 0; m < atom->num_bond[i]; m++) if (atom->bond_type[i][m] <= 0) bond_off = 1; } } if (atom->avec->angles_allow) { for (i = 0; i < atom->nlocal; i++) { if (angle_off) break; for (m = 0; m < atom->num_angle[i]; m++) if (atom->angle_type[i][m] <= 0) angle_off = 1; } } int dihedral_off = 0; if (atom->avec->dihedrals_allow) { for (i = 0; i < atom->nlocal; i++) { if (dihedral_off) break; for (m = 0; m < atom->num_dihedral[i]; m++) if (atom->dihedral_type[i][m] <= 0) dihedral_off = 1; } } int improper_off = 0; if (atom->avec->impropers_allow) { for (i = 0; i < atom->nlocal; i++) { if (improper_off) break; for (m = 0; m < atom->num_improper[i]; m++) if (atom->improper_type[i][m] <= 0) improper_off = 1; } } // set ptrs to topology build functions if (bond_off) bond_build = &Neighbor::bond_partial; else bond_build = &Neighbor::bond_all; if (angle_off) angle_build = &Neighbor::angle_partial; else angle_build = &Neighbor::angle_all; if (dihedral_off) dihedral_build = &Neighbor::dihedral_partial; else dihedral_build = &Neighbor::dihedral_all; if (improper_off) improper_build = &Neighbor::improper_partial; else improper_build = &Neighbor::improper_all; // set topology neighbor list counts to 0 // in case all are turned off but potential is still defined nbondlist = nanglelist = ndihedrallist = nimproperlist = 0; } /* ---------------------------------------------------------------------- */ int Neighbor::request(void *requestor) { if (nrequest == maxrequest) { maxrequest += RQDELTA; requests = (NeighRequest **) memory->srealloc(requests,maxrequest*sizeof(NeighRequest *), "neighbor:requests"); } requests[nrequest] = new NeighRequest(lmp); requests[nrequest]->requestor = requestor; nrequest++; return nrequest-1; } /* ---------------------------------------------------------------------- determine which pair_build function each neigh list needs based on settings of neigh request copy -> copy_from function skip -> granular function if gran with granhistory respa function if respaouter skip_from function for everything else half_from_full, half, full, gran, respaouter -> choose by newton and tri style NSQ options = newton off, newton on style BIN options = newton off, newton on and not tri, newton on and tri stlye MULTI options = same options as BIN if none of these, ptr = NULL since pair_build is not invoked for this list use "else if" b/c skip,copy can be set in addition to half,full,etc ------------------------------------------------------------------------- */ void Neighbor::choose_build(int index, NeighRequest *rq) { PairPtr pb = NULL; if (rq->copy) pb = &Neighbor::copy_from; else if (rq->skip) { if (rq->gran && lists[index]->listgranhistory) pb = &Neighbor::skip_from_granular; else if (rq->respaouter) pb = &Neighbor::skip_from_respa; else pb = &Neighbor::skip_from; } else if (rq->half_from_full) { if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton; else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton; } else if (rq->half) { if (style == NSQ) { if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton; else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton; } else if (style == BIN) { if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton; else if (triclinic == 0) pb = &Neighbor::half_bin_newton; else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri; } else if (style == MULTI) { if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton; else if (triclinic == 0) pb = &Neighbor::half_multi_newton; else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri; } } else if (rq->full) { if (style == NSQ) pb = &Neighbor::full_nsq; else if (style == BIN) pb = &Neighbor::full_bin; else if (style == MULTI) pb = &Neighbor::full_multi; } else if (rq->gran) { if (style == NSQ) { if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton; else if (newton_pair == 1) pb = &Neighbor::granular_nsq_newton; } else if (style == BIN) { if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton; else if (triclinic == 0) pb = &Neighbor::granular_bin_newton; else if (triclinic == 1) pb = &Neighbor::granular_bin_newton_tri; } else if (style == MULTI) error->all("Neighbor multi not yet enabled for granular"); } else if (rq->respaouter) { if (style == NSQ) { if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton; else if (newton_pair == 1) pb = &Neighbor::respa_nsq_newton; } else if (style == BIN) { if (newton_pair == 0) pb = &Neighbor::respa_bin_no_newton; else if (triclinic == 0) pb = &Neighbor::respa_bin_newton; else if (triclinic == 1) pb = &Neighbor::respa_bin_newton_tri; } else if (style == MULTI) error->all("Neighbor multi not yet enabled for rRESPA"); } pair_build[index] = pb; } /* ---------------------------------------------------------------------- determine which stencil_create function each neigh list needs based on settings of neigh request, only called if style != NSQ skip or copy or half_from_full -> no stencil half, gran, respaouter, full -> choose by newton and tri and dimension if none of these, ptr = NULL since this list needs no stencils use "else if" b/c skip,copy can be set in addition to half,full,etc ------------------------------------------------------------------------- */ void Neighbor::choose_stencil(int index, NeighRequest *rq) { StencilPtr sc = NULL; if (rq->skip || rq->copy || rq->half_from_full) sc = NULL; else if (rq->half || rq->gran || rq->respaouter) { if (style == BIN) { if (newton_pair == 0) { if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_no_newton; else if (dimension == 3) sc = &Neighbor::stencil_half_bin_3d_no_newton; } else if (triclinic == 0) { if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_newton; else if (dimension == 3) sc = &Neighbor::stencil_half_bin_3d_newton; } else if (triclinic == 1) { if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_newton_tri; else if (dimension == 3) sc = &Neighbor::stencil_half_bin_3d_newton_tri; } } else if (style == MULTI) { if (newton_pair == 0) { if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_no_newton; else if (dimension == 3) sc = &Neighbor::stencil_half_multi_3d_no_newton; } else if (triclinic == 0) { if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_newton; else if (dimension == 3) sc = &Neighbor::stencil_half_multi_3d_newton; } else if (triclinic == 1) { if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_newton_tri; else if (dimension == 3) sc = &Neighbor::stencil_half_multi_3d_newton_tri; } } } else if (rq->full) { if (style == BIN) { if (dimension == 2) sc = &Neighbor::stencil_full_bin_2d; else if (dimension == 3) sc = &Neighbor::stencil_full_bin_3d; } else if (style == MULTI) { if (dimension == 2) sc = &Neighbor::stencil_full_multi_2d; else if (dimension == 3) sc = &Neighbor::stencil_full_multi_3d; } } stencil_create[index] = sc; } /* ---------------------------------------------------------------------- */ void Neighbor::print_lists_of_lists() { if (comm->me == 0) { printf("Build lists = %d: ",nblist); for (int i = 0; i < nblist; i++) printf("%d ",blist[i]); printf("\n"); printf("Grow lists = %d: ",nglist); for (int i = 0; i < nglist; i++) printf("%d ",glist[i]); printf("\n"); printf("Stencil lists = %d: ",nslist); for (int i = 0; i < nslist; i++) printf("%d ",slist[i]); printf("\n"); } } /* ---------------------------------------------------------------------- */ int Neighbor::decide() { if (must_check) { int n = update->ntimestep; if (restart_check && n == output->next_restart) return 1; for (int i = 0; i < fix_check; i++) if (n == modify->fix[fixchecklist[i]]->next_reneighbor) return 1; } ago++; if (ago >= delay && ago % every == 0) { if (dist_check == 0) return 1; else return check_distance(); } else return 0; } /* ---------------------------------------------------------------------- */ int Neighbor::check_distance() { double delx,dely,delz,rsq; - int nlocal = atom->nlocal; double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; int flag = 0; - for (int i = 0; i < nlocal; i++) { - delx = x[i][0] - xhold[i][0]; - dely = x[i][1] - xhold[i][1]; - delz = x[i][2] - xhold[i][2]; - rsq = delx*delx + dely*dely + delz*delz; - if (rsq > triggersq) flag = 1; + + if (include_group) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & include_groupbit) { + delx = x[i][0] - xhold[i][0]; + dely = x[i][1] - xhold[i][1]; + delz = x[i][2] - xhold[i][2]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq > triggersq) flag = 1; + } + } else { + for (int i = 0; i < nlocal; i++) { + delx = x[i][0] - xhold[i][0]; + dely = x[i][1] - xhold[i][1]; + delz = x[i][2] - xhold[i][2]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq > triggersq) flag = 1; + } } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world); if (flagall && ago == MAX(every,delay)) ndanger++; return flagall; } /* ---------------------------------------------------------------------- build all perpetual neighbor lists every few timesteps pairwise & topology lists are created as needed ------------------------------------------------------------------------- */ void Neighbor::build() { int i; ago = 0; ncalls++; // store current atom positions if needed if (dist_check) { double **x = atom->x; int nlocal = atom->nlocal; if (nlocal > maxhold) { maxhold = atom->nmax; memory->destroy_2d_double_array(xhold); xhold = memory->create_2d_double_array(maxhold,3,"neigh:xhold"); } for (i = 0; i < nlocal; i++) { xhold[i][0] = x[i][0]; xhold[i][1] = x[i][1]; xhold[i][2] = x[i][2]; } } // if necessary, extend atom arrays in pairwise lists // only for lists with growflag set and which are used every reneighbor if (atom->nlocal > maxlocal) { maxlocal = atom->nmax; for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxlocal); } // extend atom bin list if necessary if (style != NSQ && atom->nmax > maxbin) { maxbin = atom->nmax; memory->sfree(bins); bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins"); } // invoke building of pair and molecular neighbor lists // only for pairwise lists with buildflag set for (i = 0; i < nblist; i++) (this->*pair_build[blist[i]])(lists[blist[i]]); if (atom->molecular) { if (atom->nbonds) (this->*bond_build)(); if (atom->nangles) (this->*angle_build)(); if (atom->ndihedrals) (this->*dihedral_build)(); if (atom->nimpropers) (this->*improper_build)(); } } /* ---------------------------------------------------------------------- build a single occasional pairwise neighbor list indexed by I called by other classes ------------------------------------------------------------------------- */ void Neighbor::build_one(int i) { // update stencils and grow atom arrays and bins as needed // only for relevant settings of stencilflag and growflag // do not reset maxlocal to atom->nmax, since not all lists are being grown if (lists[i]->stencilflag) { lists[i]->stencil_allocate(smax,style); (this->*stencil_create[i])(lists[i],sx,sy,sz); } if (lists[i]->growflag) lists[i]->grow(maxlocal); if (style != NSQ && atom->nmax > maxbin) { maxbin = atom->nmax; memory->sfree(bins); bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins"); } (this->*pair_build[i])(lists[i]); } /* ---------------------------------------------------------------------- setup neighbor binning parameters bin numbering in each dimension is global: 0 = 0.0 to binsize, 1 = binsize to 2*binsize, etc nbin-1,nbin,etc = bbox-binsize to binsize, bbox to bbox+binsize, etc -1,-2,etc = -binsize to 0.0, -2*size to -size, etc code will work for any binsize since next(xyz) and stencil extend as far as necessary binsize = 1/2 of cutoff is roughly optimal for orthogonal boxes: a dim must be filled exactly by integer # of bins in periodic, procs on both sides of PBC must see same bin boundary in non-periodic, coord2bin() still assumes this by use of nbin xyz for triclinic boxes: tilted simulation box cannot contain integer # of bins stencil & neigh list built differently to account for this mbinlo = lowest global bin any of my ghost atoms could fall into mbinhi = highest global bin any of my ghost atoms could fall into mbin = number of bins I need in a dimension ------------------------------------------------------------------------- */ void Neighbor::setup_bins() { // bbox lo/hi = bounding box of entire domain // bbox = size of bbox of entire domain // bsubbox lo/hi = bounding box of my subdomain extended by ghost atoms // for triclinic: // bbox bounds all 8 corners of tilted box // subdomain is in lamda coords // include dimension-dependent extension via comm->cutghost // domain->bbox() converts lamda extent to box coords and computes bbox double bbox[3],bsubboxlo[3],bsubboxhi[3]; if (triclinic == 0) { bboxlo = domain->boxlo; bboxhi = domain->boxhi; bsubboxlo[0] = domain->sublo[0] - cutghost; bsubboxlo[1] = domain->sublo[1] - cutghost; bsubboxlo[2] = domain->sublo[2] - cutghost; bsubboxhi[0] = domain->subhi[0] + cutghost; bsubboxhi[1] = domain->subhi[1] + cutghost; bsubboxhi[2] = domain->subhi[2] + cutghost; } else { bboxlo = domain->boxlo_bound; bboxhi = domain->boxhi_bound; double lo[3],hi[3]; lo[0] = domain->sublo_lamda[0] - comm->cutghost[0]; lo[1] = domain->sublo_lamda[1] - comm->cutghost[1]; lo[2] = domain->sublo_lamda[2] - comm->cutghost[2]; hi[0] = domain->subhi_lamda[0] + comm->cutghost[0]; hi[1] = domain->subhi_lamda[1] + comm->cutghost[1]; hi[2] = domain->subhi_lamda[2] + comm->cutghost[2]; domain->bbox(lo,hi,bsubboxlo,bsubboxhi); } bbox[0] = bboxhi[0] - bboxlo[0]; bbox[1] = bboxhi[1] - bboxlo[1]; bbox[2] = bboxhi[2] - bboxlo[2]; // optimal bin size is roughly 1/2 the cutoff // for BIN style, binsize = 1/2 of max neighbor cutoff // for MULTI style, binsize = 1/2 of min neighbor cutoff // special case of all cutoffs = 0.0, binsize = box size double binsize_optimal; if (binsizeflag) binsize_optimal = binsize_user; else if (style == BIN) binsize_optimal = 0.5*cutneighmax; else binsize_optimal = 0.5*cutneighmin; if (binsize_optimal == 0.0) binsize_optimal = bbox[0]; double binsizeinv = 1.0/binsize_optimal; // test for too many global bins in any dimension due to huge global domain if (bbox[0]*binsizeinv > INT_MAX || bbox[1]*binsizeinv > INT_MAX || bbox[2]*binsizeinv > INT_MAX) error->all("Domain too large for neighbor bins"); // create actual bins // always have one bin even if cutoff > bbox // for 2d, nbinz = 1 nbinx = static_cast (bbox[0]*binsizeinv); nbiny = static_cast (bbox[1]*binsizeinv); if (dimension == 3) nbinz = static_cast (bbox[2]*binsizeinv); else nbinz = 1; if (nbinx == 0) nbinx = 1; if (nbiny == 0) nbiny = 1; if (nbinz == 0) nbinz = 1; // compute actual bin size for nbins to fit into box exactly // error if actual bin size << cutoff, since will create a zillion bins // this happens when nbin = 1 and box size << cutoff // typically due to non-periodic, flat system in a particular dim // in that extreme case, should use NSQ not BIN neighbor style binsizex = bbox[0]/nbinx; binsizey = bbox[1]/nbiny; binsizez = bbox[2]/nbinz; bininvx = 1.0 / binsizex; bininvy = 1.0 / binsizey; bininvz = 1.0 / binsizez; if (binsize_optimal*bininvx > CUT2BIN_RATIO || binsize_optimal*bininvy > CUT2BIN_RATIO || binsize_optimal*bininvz > CUT2BIN_RATIO) error->all("Cannot use neighbor bins - box size << cutoff"); // mbinlo/hi = lowest and highest global bins my ghost atoms could be in // coord = lowest and highest values of coords for my ghost atoms // static_cast(-1.5) = -1, so subract additional -1 // add in SMALL for round-off safety int mbinxhi,mbinyhi,mbinzhi; double coord; coord = bsubboxlo[0] - SMALL*bbox[0]; mbinxlo = static_cast ((coord-bboxlo[0])*bininvx); if (coord < bboxlo[0]) mbinxlo = mbinxlo - 1; coord = bsubboxhi[0] + SMALL*bbox[0]; mbinxhi = static_cast ((coord-bboxlo[0])*bininvx); coord = bsubboxlo[1] - SMALL*bbox[1]; mbinylo = static_cast ((coord-bboxlo[1])*bininvy); if (coord < bboxlo[1]) mbinylo = mbinylo - 1; coord = bsubboxhi[1] + SMALL*bbox[1]; mbinyhi = static_cast ((coord-bboxlo[1])*bininvy); if (dimension == 3) { coord = bsubboxlo[2] - SMALL*bbox[2]; mbinzlo = static_cast ((coord-bboxlo[2])*bininvz); if (coord < bboxlo[2]) mbinzlo = mbinzlo - 1; coord = bsubboxhi[2] + SMALL*bbox[2]; mbinzhi = static_cast ((coord-bboxlo[2])*bininvz); } // extend bins by 1 to insure stencil extent is included // if 2d, only 1 bin in z mbinxlo = mbinxlo - 1; mbinxhi = mbinxhi + 1; mbinx = mbinxhi - mbinxlo + 1; mbinylo = mbinylo - 1; mbinyhi = mbinyhi + 1; mbiny = mbinyhi - mbinylo + 1; if (dimension == 3) { mbinzlo = mbinzlo - 1; mbinzhi = mbinzhi + 1; } else mbinzlo = mbinzhi = 0; mbinz = mbinzhi - mbinzlo + 1; // memory for bin ptrs mbins = mbinx*mbiny*mbinz; if (mbins > maxhead) { maxhead = mbins; memory->sfree(binhead); binhead = (int *) memory->smalloc(maxhead*sizeof(int),"neigh:binhead"); } // create stencil of bins to search over in neighbor list construction // sx,sy,sz = max range of stencil extent // smax = // stencil is empty if cutneighmax = 0.0 sx = static_cast (cutneighmax*bininvx); if (sx*binsizex < cutneighmax) sx++; sy = static_cast (cutneighmax*bininvy); if (sy*binsizey < cutneighmax) sy++; sz = static_cast (cutneighmax*bininvz); if (sz*binsizez < cutneighmax) sz++; if (dimension == 2) sz = 0; smax = (2*sx+1) * (2*sy+1) * (2*sz+1); // create stencils for pairwise neighbor lists // only done for lists with stencilflag and buildflag set for (int i = 0; i < nslist; i++) { lists[slist[i]]->stencil_allocate(smax,style); (this->*stencil_create[slist[i]])(lists[slist[i]],sx,sy,sz); } } /* ---------------------------------------------------------------------- compute closest distance between central bin (0,0,0) and bin (i,j,k) ------------------------------------------------------------------------- */ double Neighbor::bin_distance(int i, int j, int k) { double delx,dely,delz; if (i > 0) delx = (i-1)*binsizex; else if (i == 0) delx = 0.0; else delx = (i+1)*binsizex; if (j > 0) dely = (j-1)*binsizey; else if (j == 0) dely = 0.0; else dely = (j+1)*binsizey; if (k > 0) delz = (k-1)*binsizez; else if (k == 0) delz = 0.0; else delz = (k+1)*binsizez; return (delx*delx + dely*dely + delz*delz); } /* ---------------------------------------------------------------------- set neighbor style and skin distance ------------------------------------------------------------------------- */ void Neighbor::set(int narg, char **arg) { if (narg != 2) error->all("Illegal neighbor command"); skin = atof(arg[0]); if (skin < 0.0) error->all("Illegal neighbor command"); if (strcmp(arg[1],"nsq") == 0) style = NSQ; else if (strcmp(arg[1],"bin") == 0) style = BIN; else if (strcmp(arg[1],"multi") == 0) style = MULTI; else error->all("Illegal neighbor command"); } /* ---------------------------------------------------------------------- modify parameters of the pair-wise neighbor build ------------------------------------------------------------------------- */ void Neighbor::modify_params(int narg, char **arg) { int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"every") == 0) { if (iarg+2 > narg) error->all("Illegal neigh_modify command"); every = atoi(arg[iarg+1]); if (every <= 0) error->all("Illegal neigh_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"delay") == 0) { if (iarg+2 > narg) error->all("Illegal neigh_modify command"); delay = atoi(arg[iarg+1]); if (delay < 0) error->all("Illegal neigh_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"check") == 0) { if (iarg+2 > narg) error->all("Illegal neigh_modify command"); if (strcmp(arg[iarg+1],"yes") == 0) dist_check = 1; else if (strcmp(arg[iarg+1],"no") == 0) dist_check = 0; else error->all("Illegal neigh_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"page") == 0) { if (iarg+2 > narg) error->all("Illegal neigh_modify command"); pgsize = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"one") == 0) { if (iarg+2 > narg) error->all("Illegal neigh_modify command"); oneatom = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"binsize") == 0) { if (iarg+2 > narg) error->all("Illegal neigh_modify command"); binsize_user = atof(arg[iarg+1]); if (binsize_user <= 0.0) binsizeflag = 0; else binsizeflag = 1; iarg += 2; + } else if (strcmp(arg[iarg],"include") == 0) { + if (iarg+2 > narg) error->all("Illegal neigh_modify command"); + include_group = group->find(arg[iarg+1]); + if (include_group == -1) + error->all("Invalid group ID in neigh_modify command"); + iarg += 2; } else if (strcmp(arg[iarg],"exclude") == 0) { if (iarg+2 > narg) error->all("Illegal neigh_modify command"); if (strcmp(arg[iarg+1],"type") == 0) { if (iarg+4 > narg) error->all("Illegal neigh_modify command"); if (nex_type == maxex_type) { maxex_type += EXDELTA; ex1_type = (int *) memory->srealloc(ex1_type,maxex_type*sizeof(int), "neigh:ex1_type"); ex2_type = (int *) memory->srealloc(ex2_type,maxex_type*sizeof(int), "neigh:ex2_type"); } ex1_type[nex_type] = atoi(arg[iarg+2]); ex2_type[nex_type] = atoi(arg[iarg+3]); nex_type++; iarg += 4; } else if (strcmp(arg[iarg+1],"group") == 0) { if (iarg+4 > narg) error->all("Illegal neigh_modify command"); if (nex_group == maxex_group) { maxex_group += EXDELTA; ex1_group = (int *) memory->srealloc(ex1_group,maxex_group*sizeof(int), "neigh:ex1_group"); ex2_group = (int *) memory->srealloc(ex2_group,maxex_group*sizeof(int), "neigh:ex2_group"); } ex1_group[nex_group] = group->find(arg[iarg+2]); ex2_group[nex_group] = group->find(arg[iarg+3]); if (ex1_group[nex_group] == -1 || ex2_group[nex_group] == -1) error->all("Invalid group ID in neigh_modify command"); nex_group++; iarg += 4; } else if (strcmp(arg[iarg+1],"molecule") == 0) { if (iarg+3 > narg) error->all("Illegal neigh_modify command"); if (atom->molecular == 0) { char *str = (char *) "Must use molecular atom style with neigh_modify exclude molecule"; error->all(str); } if (nex_mol == maxex_mol) { maxex_mol += EXDELTA; ex_mol_group = (int *) memory->srealloc(ex_mol_group,maxex_mol*sizeof(int), "neigh:ex_mol_group"); } ex_mol_group[nex_mol] = group->find(arg[iarg+2]); if (ex_mol_group[nex_mol] == -1) error->all("Invalid group ID in neigh_modify command"); nex_mol++; iarg += 3; } else if (strcmp(arg[iarg+1],"none") == 0) { nex_type = nex_group = nex_mol = 0; iarg += 2; } else error->all("Illegal neigh_modify command"); } else error->all("Illegal neigh_modify command"); } } /* ---------------------------------------------------------------------- determine if atom j is in special list of atom i if it is not, return 0 if it is and special flag is 0 (both coeffs are 0.0), return -1 if it is and special flag is 1 (both coeffs are 1.0), return 0 if it is and special flag is 2 (otherwise), return 1,2,3 for which neighbor it is (and which coeff it maps to) ------------------------------------------------------------------------- */ int Neighbor::find_special(int i, int j) { int *list = atom->special[i]; int n1 = atom->nspecial[i][0]; int n2 = atom->nspecial[i][1]; int n3 = atom->nspecial[i][2]; int tag = atom->tag[j]; for (int i = 0; i < n3; i++) { if (list[i] == tag) { if (i < n1) { if (special_flag[1] == 0) return -1; else if (special_flag[1] == 1) return 0; else return 1; } else if (i < n2) { if (special_flag[2] == 0) return -1; else if (special_flag[2] == 1) return 0; else return 2; } else { if (special_flag[3] == 0) return -1; else if (special_flag[3] == 1) return 0; else return 3; } } } return 0; } /* ---------------------------------------------------------------------- bin owned and ghost atoms ------------------------------------------------------------------------- */ void Neighbor::bin_atoms() { int i,ibin,nall; - double **x; - - nall = atom->nlocal + atom->nghost; - x = atom->x; for (i = 0; i < mbins; i++) binhead[i] = -1; // bin in reverse order so linked list will be in forward order // also puts ghost atoms at end of list, which is necessary - for (i = nall-1; i >= 0; i--) { - ibin = coord2bin(x[i]); - bins[i] = binhead[ibin]; - binhead[ibin] = i; - } - - /* - for (i = nlocal; i < nall; i++) { - ibin = coord2bin(x[i]); - bins[i] = binhead[ibin]; - binhead[ibin] = i; - } + double **x = atom->x; + int *mask = atom->mask; + nall = atom->nlocal + atom->nghost; - for (i = 0; i < nlocal; i++) { - ibin = coord2bin(x[i]); - bins[i] = binhead[ibin]; - binhead[ibin] = i; + if (include_group) { + for (i = nall-1; i >= 0; i--) + if (mask[i] & include_groupbit) { + ibin = coord2bin(x[i]); + bins[i] = binhead[ibin]; + binhead[ibin] = i; + } + } else { + for (i = nall-1; i >= 0; i--) { + ibin = coord2bin(x[i]); + bins[i] = binhead[ibin]; + binhead[ibin] = i; + } } - */ } /* ---------------------------------------------------------------------- convert atom coords into local bin # for orthogonal, only ghost atoms will have coord >= bboxhi or coord < bboxlo take special care to insure ghosts are put in correct bins this is necessary so that both procs on either side of PBC treat a pair of atoms straddling the PBC in a consistent way for triclinic, doesn't matter since stencil & neigh list built differently ------------------------------------------------------------------------- */ int Neighbor::coord2bin(double *x) { int ix,iy,iz; if (x[0] >= bboxhi[0]) ix = static_cast ((x[0]-bboxhi[0])*bininvx) + nbinx - mbinxlo; else if (x[0] >= bboxlo[0]) ix = static_cast ((x[0]-bboxlo[0])*bininvx) - mbinxlo; else ix = static_cast ((x[0]-bboxlo[0])*bininvx) - mbinxlo - 1; if (x[1] >= bboxhi[1]) iy = static_cast ((x[1]-bboxhi[1])*bininvy) + nbiny - mbinylo; else if (x[1] >= bboxlo[1]) iy = static_cast ((x[1]-bboxlo[1])*bininvy) - mbinylo; else iy = static_cast ((x[1]-bboxlo[1])*bininvy) - mbinylo - 1; if (x[2] >= bboxhi[2]) iz = static_cast ((x[2]-bboxhi[2])*bininvz) + nbinz - mbinzlo; else if (x[2] >= bboxlo[2]) iz = static_cast ((x[2]-bboxlo[2])*bininvz) - mbinzlo; else iz = static_cast ((x[2]-bboxlo[2])*bininvz) - mbinzlo - 1; return (iz*mbiny*mbinx + iy*mbinx + ix + 1); } /* ---------------------------------------------------------------------- test if atom pair i,j is excluded from neighbor list due to type, group, molecule settings from neigh_modify command return 1 if should be excluded, 0 if included ------------------------------------------------------------------------- */ int Neighbor::exclusion(int i, int j, int itype, int jtype, int *mask, int *molecule) { int m; if (nex_type && ex_type[itype][jtype]) return 1; if (nex_group) { for (m = 0; m < nex_group; m++) { if (mask[i] & ex1_bit[m] && mask[j] & ex2_bit[m]) return 1; if (mask[i] & ex2_bit[m] && mask[j] & ex1_bit[m]) return 1; } } if (nex_mol) { for (m = 0; m < nex_mol; m++) if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] && molecule[i] == molecule[j]) return 1; } return 0; } /* ---------------------------------------------------------------------- return # of bytes of allocated memory ------------------------------------------------------------------------- */ double Neighbor::memory_usage() { double bytes = 0.0; bytes += maxhold*3 * sizeof(double); if (style != NSQ) { bytes += maxbin * sizeof(int); bytes += maxhead * sizeof(int); } for (int i = 0; i < nlist; i++) bytes += lists[i]->memory_usage(); bytes += maxbond*3 * sizeof(int); bytes += maxangle*4 * sizeof(int); bytes += maxdihedral*5 * sizeof(int); bytes += maximproper*5 * sizeof(int); return bytes; } diff --git a/src/neighbor.h b/src/neighbor.h index 39f578898..3fc62c442 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -1,243 +1,246 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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 NEIGHBOR_H #define NEIGHBOR_H #include "pointers.h" namespace LAMMPS_NS { class Neighbor : protected Pointers { public: int style; // 0,1,2 = nsq, bin, multi int every; // build every this many steps int delay; // delay build for this many steps int dist_check; // 0 = always build, 1 = only if 1/2 dist int ago; // how many steps ago neighboring occurred int pgsize; // size of neighbor page int oneatom; // max # of neighbors for one atom double skin; // skin distance double cutghost; // distance for acquiring ghost atoms double *cuttype; // for each type, max neigh cut w/ others int ncalls; // # of times build has been called int ndanger; // # of dangerous builds int nrequest; // requests for pairwise neighbor lists class NeighRequest **requests; // from Pair, Fix, Compute, Command classes int maxrequest; int old_style; // previous run info to avoid int old_nrequest; // re-creation of pairwise neighbor lists class NeighRequest **old_requests; int nlist; // pairwise neighbor lists class NeighList **lists; int nbondlist; // list of bonds to compute int **bondlist; int nanglelist; // list of angles to compute int **anglelist; int ndihedrallist; // list of dihedrals to compute int **dihedrallist; int nimproperlist; // list of impropers to compute int **improperlist; Neighbor(class LAMMPS *); ~Neighbor(); void init(); int request(void *); // another class requests a neighbor list void print_lists_of_lists(); // debug print out int decide(); // decide whether to build or not int check_distance(); // check max distance moved since last build void setup_bins(); // setup bins based on box and cutoff void build(); // create all neighbor lists (pair,bond) void build_one(int); // create a single neighbor list void set(int, char **); // set neighbor style and skin distance void modify_params(int, char**); // modify parameters that control builds double memory_usage(); private: int me,nprocs; int maxlocal; // size of atom-based NeighList arrays int maxbond,maxangle,maxdihedral,maximproper; // size of bond lists int must_check; // 1 if must check other classes to reneigh int restart_check; // 1 if restart enabled, 0 if no int fix_check; // # of fixes that induce reneigh int *fixchecklist; // which fixes to check double **cutneighsq; // neighbor cutneigh sq for each type pair double cutneighmin; // min neighbor cutoff (for cutforce > 0) double cutneighmax; // max neighbor cutoff for all type pairs double cutneighmaxsq; // cutneighmax squared double *cuttypesq; // cuttype squared double triggersq; // trigger = build when atom moves this dist double **xhold; // atom coords at last neighbor build int maxhold; // size of xhold array int nbinx,nbiny,nbinz; // # of global bins int *bins; // ptr to next atom in each bin int maxbin; // size of bins array int *binhead; // ptr to 1st atom in each bin int maxhead; // size of binhead array int mbins; // # of local bins and offset int mbinx,mbiny,mbinz; int mbinxlo,mbinylo,mbinzlo; int binsizeflag; // user-chosen bin size double binsize_user; double binsizex,binsizey,binsizez; // actual bin sizes and inverse sizes double bininvx,bininvy,bininvz; int sx,sy,sz,smax; // bin stencil extents int dimension; // 2/3 for 2d/3d int triclinic; // 0 if domain is orthog, 1 if triclinic int newton_pair; // 0 if newton off, 1 if on for pairwise double *bboxlo,*bboxhi; // copy of full domain bounding box double inner[2],middle[2]; // rRESPA cutoffs for extra lists double cut_inner_sq; // outer cutoff for inner neighbor list double cut_middle_sq; // outer cutoff for middle neighbor list double cut_middle_inside_sq; // inner cutoff for middle neighbor list int special_flag[4]; // flags for 1-2, 1-3, 1-4 neighbors + int include_group; // only form neighbor lists on this group + int include_groupbit; + int exclude; // 0 if no type/group exclusions, 1 if yes int nex_type; // # of entries in type exclusion list int maxex_type; // max # in type list int *ex1_type,*ex2_type; // pairs of types to exclude int **ex_type; // 2d array of excluded type pairs int nex_group; // # of entries in group exclusion list int maxex_group; // max # in group list int *ex1_group,*ex2_group; // pairs of group #'s to exclude int *ex1_bit,*ex2_bit; // pairs of group bits to exclude int nex_mol; // # of entries in molecule exclusion list int maxex_mol; // max # in molecule list int *ex_mol_group; // molecule group #'s to exclude int *ex_mol_bit; // molecule group bits to exclude int nblist,nglist,nslist; // # of pairwise neigh lists of various kinds int *blist; // lists to build every reneighboring int *glist; // lists to grow atom arrays every reneigh int *slist; // lists to grow stencil arrays every reneigh void bin_atoms(); // bin all atoms double bin_distance(int, int, int); // distance between binx int coord2bin(double *); // mapping atom coord to a bin int find_special(int, int); // look for special neighbor int exclusion(int, int, int, int, int *, int *); // test for pair exclusion void choose_build(int, class NeighRequest *); void choose_stencil(int, class NeighRequest *); // pairwise build functions typedef void (Neighbor::*PairPtr)(class NeighList *); PairPtr *pair_build; void half_nsq_no_newton(class NeighList *); void half_nsq_newton(class NeighList *); void half_bin_no_newton(class NeighList *); void half_bin_newton(class NeighList *); void half_bin_newton_tri(class NeighList *); void half_multi_no_newton(class NeighList *); void half_multi_newton(class NeighList *); void half_multi_newton_tri(class NeighList *); void full_nsq(class NeighList *); void full_bin(class NeighList *); void full_multi(class NeighList *); void half_from_full_no_newton(class NeighList *); void half_from_full_newton(class NeighList *); void skip_from(class NeighList *); void skip_from_granular(class NeighList *); void skip_from_respa(class NeighList *); void copy_from(class NeighList *); void granular_nsq_no_newton(class NeighList *); void granular_nsq_newton(class NeighList *); void granular_bin_no_newton(class NeighList *); void granular_bin_newton(class NeighList *); void granular_bin_newton_tri(class NeighList *); void respa_nsq_no_newton(class NeighList *); void respa_nsq_newton(class NeighList *); void respa_bin_no_newton(class NeighList *); void respa_bin_newton(class NeighList *); void respa_bin_newton_tri(class NeighList *); // pairwise stencil creation functions typedef void (Neighbor::*StencilPtr)(class NeighList *, int, int, int); StencilPtr *stencil_create; void stencil_half_bin_2d_no_newton(class NeighList *, int, int, int); void stencil_half_bin_3d_no_newton(class NeighList *, int, int, int); void stencil_half_bin_2d_newton(class NeighList *, int, int, int); void stencil_half_bin_3d_newton(class NeighList *, int, int, int); void stencil_half_bin_2d_newton_tri(class NeighList *, int, int, int); void stencil_half_bin_3d_newton_tri(class NeighList *, int, int, int); void stencil_half_multi_2d_no_newton(class NeighList *, int, int, int); void stencil_half_multi_3d_no_newton(class NeighList *, int, int, int); void stencil_half_multi_2d_newton(class NeighList *, int, int, int); void stencil_half_multi_3d_newton(class NeighList *, int, int, int); void stencil_half_multi_2d_newton_tri(class NeighList *, int, int, int); void stencil_half_multi_3d_newton_tri(class NeighList *, int, int, int); void stencil_full_bin_2d(class NeighList *, int, int, int); void stencil_full_bin_3d(class NeighList *, int, int, int); void stencil_full_multi_2d(class NeighList *, int, int, int); void stencil_full_multi_3d(class NeighList *, int, int, int); // topology build functions typedef void (Neighbor::*BondPtr)(); // ptrs to topology build functions BondPtr bond_build; // ptr to bond list functions void bond_all(); // bond list with all bonds void bond_partial(); // exclude certain bonds BondPtr angle_build; // ptr to angle list functions void angle_all(); // angle list with all angles void angle_partial(); // exclude certain angles BondPtr dihedral_build; // ptr to dihedral list functions void dihedral_all(); // dihedral list with all dihedrals void dihedral_partial(); // exclude certain dihedrals BondPtr improper_build; // ptr to improper list functions void improper_all(); // improper list with all impropers void improper_partial(); // exclude certain impropers }; } #endif diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index 6ab2b0bbe..24985e302 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -1,722 +1,722 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov 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. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Paul Crozier (SNL) ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_lj_cut.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "update.h" #include "integrate.h" #include "respa.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) /* ---------------------------------------------------------------------- */ PairLJCut::PairLJCut(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; } /* ---------------------------------------------------------------------- */ PairLJCut::~PairLJCut() { if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); memory->destroy_2d_double_array(cut); memory->destroy_2d_double_array(epsilon); memory->destroy_2d_double_array(sigma); memory->destroy_2d_double_array(lj1); memory->destroy_2d_double_array(lj2); memory->destroy_2d_double_array(lj3); memory->destroy_2d_double_array(lj4); memory->destroy_2d_double_array(offset); } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj; int *ilist,*jlist,*numneigh,**firstneigh; evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } if (eflag) { evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,0.0,fpair,delx,dely,delz); } } } if (vflag_fdotr) virial_compute(); } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_inner() { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; int *ilist,*jlist,*numneigh,**firstneigh; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = listinner->inum; ilist = listinner->ilist; numneigh = listinner->numneigh; firstneigh = listinner->firstneigh; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cut_out_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; jtype = type[j]; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; if (rsq > cut_out_on_sq) { rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw); } f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } } } } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_middle() { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; int *ilist,*jlist,*numneigh,**firstneigh; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = listmiddle->inum; ilist = listmiddle->ilist; numneigh = listmiddle->numneigh; firstneigh = listmiddle->firstneigh; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; double cut_out_on = cut_respa[2]; double cut_out_off = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; jtype = type[j]; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; if (rsq < cut_in_on_sq) { rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; fpair *= rsw*rsw*(3.0 - 2.0*rsw); } if (rsq > cut_out_on_sq) { rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); } f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } } } } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_outer(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; int *ilist,*jlist,*numneigh,**firstneigh; evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = 0; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = listouter->inum; ilist = listouter->ilist; numneigh = listouter->numneigh; firstneigh = listouter->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { if (rsq > cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; if (rsq < cut_in_on_sq) { rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; fpair *= rsw*rsw*(3.0 - 2.0*rsw); } f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } } if (eflag) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } if (vflag) { if (rsq <= cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; } else if (rsq < cut_in_on_sq) fpair = factor_lj*forcelj*r2inv; } if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,0.0,fpair,delx,dely,delz); } } } } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairLJCut::allocate() { allocated = 1; int n = atom->ntypes; setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); epsilon = memory->create_2d_double_array(n+1,n+1,"pair:epsilon"); sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); lj1 = memory->create_2d_double_array(n+1,n+1,"pair:lj1"); lj2 = memory->create_2d_double_array(n+1,n+1,"pair:lj2"); lj3 = memory->create_2d_double_array(n+1,n+1,"pair:lj3"); lj4 = memory->create_2d_double_array(n+1,n+1,"pair:lj4"); offset = memory->create_2d_double_array(n+1,n+1,"pair:offset"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairLJCut::settings(int narg, char **arg) { if (narg != 1) error->all("Illegal pair_style command"); cut_global = atof(arg[0]); // reset cutoffs that have been explicitly set if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; } } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairLJCut::coeff(int narg, char **arg) { if (narg < 4 || narg > 5) error->all("Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); double epsilon_one = atof(arg[2]); double sigma_one = atof(arg[3]); double cut_one = cut_global; if (narg == 5) cut_one = atof(arg[4]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut[i][j] = cut_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all("Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairLJCut::init_style() { // request regular or rRESPA neighbor lists int irequest; if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) { int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; if (respa == 0) irequest = neighbor->request(this); else if (respa == 1) { irequest = neighbor->request(this); neighbor->requests[irequest]->id = 1; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respainner = 1; irequest = neighbor->request(this); neighbor->requests[irequest]->id = 3; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respaouter = 1; } else { irequest = neighbor->request(this); neighbor->requests[irequest]->id = 1; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respainner = 1; irequest = neighbor->request(this); neighbor->requests[irequest]->id = 2; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respamiddle = 1; irequest = neighbor->request(this); neighbor->requests[irequest]->id = 3; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respaouter = 1; } } else irequest = neighbor->request(this); } /* ---------------------------------------------------------------------- neighbor callback to inform pair style of neighbor list to use regular or rRESPA ------------------------------------------------------------------------- */ void PairLJCut::init_list(int id, NeighList *ptr) { if (id == 0) list = ptr; else if (id == 1) listinner = ptr; else if (id == 2) listmiddle = ptr; else if (id == 3) listouter = ptr; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairLJCut::init_one(int i, int j) { if (setflag[i][j] == 0) { epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], sigma[i][i],sigma[j][j]); sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); cut[i][j] = mix_distance(cut[i][i],cut[j][j]); } lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); if (offset_flag) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; lj3[j][i] = lj3[i][j]; lj4[j][i] = lj4[i][j]; offset[j][i] = offset[i][j]; // set & error check interior rRESPA cutoff if (strcmp(update->integrate_style,"respa") == 0) { if (((Respa *) update->integrate)->level_inner >= 0) { cut_respa = ((Respa *) update->integrate)->cutoff; if (cut[i][j] < cut_respa[3]) error->all("Pair cutoff < Respa interior cutoff"); } } else cut_respa = NULL; // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce if (tail_flag) { int *type = atom->type; int nlocal = atom->nlocal; double count[2],all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); double PI = 4.0*atan(1.0); double sig2 = sigma[i][j]*sigma[i][j]; double sig6 = sig2*sig2*sig2; double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; double rc6 = rc3*rc3; double rc9 = rc3*rc6; etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); ptail_ij = 16.0*PI*all[0]*all[1]*epsilon[i][j] * sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); } return cut[i][j]; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJCut::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&epsilon[i][j],sizeof(double),1,fp); fwrite(&sigma[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJCut::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) { fread(&epsilon[i][j],sizeof(double),1,fp); fread(&sigma[i][j],sizeof(double),1,fp); fread(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJCut::write_restart_settings(FILE *fp) { fwrite(&cut_global,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJCut::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { fread(&cut_global,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- */ double PairLJCut::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double r2inv,r6inv,forcelj,philj; r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fforce = factor_lj*forcelj*r2inv; philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; return factor_lj*philj; } /* ---------------------------------------------------------------------- */ void *PairLJCut::extract(char *str) { - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"diameter") == 0) return (void *) sigma; return NULL; }