diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp index ac877f0bcc..589fd83af8 100644 --- a/src/fix_shake.cpp +++ b/src/fix_shake.cpp @@ -38,6 +38,10 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; +// allocate space for static class variable + +FixShake *FixShake::fsptr; + #define BIG 1.0e20 #define MASSDELTA 0.1 @@ -521,7 +525,7 @@ void FixShake::post_force(int vflag) int m; for (int i = 0; i < nlist; i++) { m = list[i]; - if (shake_flag[m] == 2) shake2(m); + if (shake_flag[m] == 2) shake(m); else if (shake_flag[m] == 3) shake3(m); else if (shake_flag[m] == 4) shake4(m); else shake3angle(m); @@ -560,7 +564,7 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop) int m; for (int i = 0; i < nlist; i++) { m = list[i]; - if (shake_flag[m] == 2) shake2(m); + if (shake_flag[m] == 2) shake(m); else if (shake_flag[m] == 3) shake3(m); else if (shake_flag[m] == 4) shake4(m); else shake3angle(m); @@ -656,7 +660,7 @@ void FixShake::find_clusters() int max = 0; for (i = 0; i < nlocal; i++) max = MAX(max,nspecial[i][0]); - int *npartner,*nshake; + int *npartner; memory->create(npartner,nlocal,"shake:npartner"); memory->create(nshake,nlocal,"shake:nshake"); @@ -719,10 +723,7 @@ void FixShake::find_clusters() } } - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"shake:buf"); // fill buffer with info @@ -745,39 +746,9 @@ void FixShake::find_clusters() } // cycle buffer around ring of procs back to self - // when receive buffer, scan bond partner IDs for atoms I own - // if I own partner: - // fill in mask and type and massflag - // search for bond with 1st atom and fill in bondtype - messtag = 1; - for (loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - m = atom->map(buf[i+1]); - if (m >= 0 && m < nlocal) { - buf[i+2] = mask[m]; - buf[i+3] = type[m]; - if (nmass) { - if (rmass) massone = rmass[m]; - else massone = mass[type[m]]; - buf[i+4] = masscheck(massone); - } - if (buf[i+5] == 0) { - n = bondfind(m,buf[i],buf[i+1]); - if (n >= 0) buf[i+5] = bond_type[m][n]; - } - } - i += nper; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + fsptr = this; + comm->ring(size,sizeof(int),buf,1,ring_bonds,buf); // store partner info returned to me @@ -793,8 +764,7 @@ void FixShake::find_clusters() m += nper; } - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); // error check for unfilled partner info // if partner_type not set, is an error @@ -885,10 +855,7 @@ void FixShake::find_clusters() } } - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"shake:buf"); // fill buffer with info @@ -905,28 +872,12 @@ void FixShake::find_clusters() } // cycle buffer around ring of procs back to self - // when receive buffer, scan bond partner IDs for atoms I own - // if I own partner, fill in nshake value - messtag = 2; - for (loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - m = atom->map(buf[i+1]); - if (m >= 0 && m < nlocal) buf[i+2] = nshake[m]; - i += 3; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + fsptr = this; + comm->ring(size,sizeof(int),buf,2,ring_nshake,buf); // store partner info returned to me - + m = 0; while (m < size) { i = atom->map(buf[m]); @@ -936,8 +887,7 @@ void FixShake::find_clusters() m += 3; } - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); // ----------------------------------------------------- // error checks @@ -1050,10 +1000,7 @@ void FixShake::find_clusters() } } - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"shake:buf"); // fill buffer with info @@ -1079,37 +1026,11 @@ void FixShake::find_clusters() } // cycle buffer around ring of procs back to self - // when receive buffer, scan for ID that I own - // if I own ID, fill in shake array values - messtag = 3; - for (loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - m = atom->map(buf[i]); - if (m >= 0 && m < nlocal) { - shake_flag[m] = buf[i+1]; - shake_atom[m][0] = buf[i+2]; - shake_atom[m][1] = buf[i+3]; - shake_atom[m][2] = buf[i+4]; - shake_atom[m][3] = buf[i+5]; - shake_type[m][0] = buf[i+6]; - shake_type[m][1] = buf[i+7]; - shake_type[m][2] = buf[i+8]; - } - i += 9; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + fsptr = this; + comm->ring(size,sizeof(int),buf,3,ring_shake,NULL); - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); // ----------------------------------------------------- // free local memory @@ -1196,6 +1117,99 @@ void FixShake::find_clusters() } } +/* ---------------------------------------------------------------------- + when receive buffer, scan bond partner IDs for atoms I own + if I own partner: + fill in mask and type and massflag + search for bond with 1st atom and fill in bondtype +------------------------------------------------------------------------- */ + +void FixShake::ring_bonds(int ndatum, char *cbuf) +{ + Atom *atom = fsptr->atom; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *mask = atom->mask; + int **bond_type = atom->bond_type; + int *type = atom->type; + int nlocal = atom->nlocal; + int nmass = fsptr->nmass; + + int *buf = (int *) cbuf; + int m,n; + double massone; + + for (int i = 0; i < ndatum; i += 6) { + m = atom->map(buf[i+1]); + if (m >= 0 && m < nlocal) { + buf[i+2] = mask[m]; + buf[i+3] = type[m]; + if (nmass) { + if (rmass) massone = rmass[m]; + else massone = mass[type[m]]; + buf[i+4] = fsptr->masscheck(massone); + } + if (buf[i+5] == 0) { + n = fsptr->bondfind(m,buf[i],buf[i+1]); + if (n >= 0) buf[i+5] = bond_type[m][n]; + } + } + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan bond partner IDs for atoms I own + if I own partner, fill in nshake value +------------------------------------------------------------------------- */ + +void FixShake::ring_nshake(int ndatum, char *cbuf) +{ + Atom *atom = fsptr->atom; + int nlocal = atom->nlocal; + + int *nshake = fsptr->nshake; + + int *buf = (int *) cbuf; + int m; + + for (int i = 0; i < ndatum; i += 3) { + m = atom->map(buf[i+1]); + if (m >= 0 && m < nlocal) buf[i+2] = nshake[m]; + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan bond partner IDs for atoms I own + if I own partner, fill in nshake value +------------------------------------------------------------------------- */ + +void FixShake::ring_shake(int ndatum, char *cbuf) +{ + Atom *atom = fsptr->atom; + int nlocal = atom->nlocal; + + int *shake_flag = fsptr->shake_flag; + int **shake_atom = fsptr->shake_atom; + int **shake_type = fsptr->shake_type; + + int *buf = (int *) cbuf; + int m; + + for (int i = 0; i < ndatum; i += 9) { + m = atom->map(buf[i]); + if (m >= 0 && m < nlocal) { + shake_flag[m] = buf[i+1]; + shake_atom[m][0] = buf[i+2]; + shake_atom[m][1] = buf[i+3]; + shake_atom[m][2] = buf[i+4]; + shake_atom[m][3] = buf[i+5]; + shake_type[m][0] = buf[i+6]; + shake_type[m][1] = buf[i+7]; + shake_type[m][2] = buf[i+8]; + } + } +} + /* ---------------------------------------------------------------------- check if massone is within MASSDELTA of any mass in mass_list return 1 if yes, 0 if not @@ -1252,7 +1266,7 @@ void FixShake::unconstrained_update_respa(int ilevel) // for levels > 0 this includes more than one velocity update // xshake = predicted position from call to this routine at level N = // x + dt0 (v + dtN/m fN + 1/2 dt(N-1)/m f(N-1) + ... + 1/2 dt0/m f0) - // also set dtfsq = dt0*dtN so that shake2,shake3,etc can use it + // also set dtfsq = dt0*dtN so that shake,shake3,etc can use it double ***f_level = ((FixRespa *) modify->fix[ifix_respa])->f_level; dtfsq = dtf_inner * step_respa[ilevel]; @@ -1298,7 +1312,7 @@ void FixShake::unconstrained_update_respa(int ilevel) /* ---------------------------------------------------------------------- */ -void FixShake::shake2(int m) +void FixShake::shake(int m) { int nlist,list[2]; double v[6]; diff --git a/src/fix_shake.h b/src/fix_shake.h index c10d8c127f..fcd2a93f38 100644 --- a/src/fix_shake.h +++ b/src/fix_shake.h @@ -82,6 +82,7 @@ class FixShake : public Fix { // for angle cluster, 3rd value // is angletype double **xshake; // unconstrained atom coords + int *nshake; // count int vflag; // virial flag double dtv,dtfsq; // timesteps for trial move @@ -102,13 +103,21 @@ class FixShake : public Fix { int masscheck(double); void unconstrained_update(); void unconstrained_update_respa(int); - void shake2(int); + void shake(int); void shake3(int); void shake4(int); void shake3angle(int); void stats(); int bondfind(int, int, int); int anglefind(int, int, int); + + // static variable for ring communication callback to access class data + // callback functions for ring communication + + static FixShake *fsptr; + static void ring_bonds(int, char *); + static void ring_nshake(int, char *); + static void ring_shake(int, char *); }; } diff --git a/src/read_data.cpp b/src/read_data.cpp index 0065472021..f03805cbca 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -34,9 +34,9 @@ #include "angle.h" #include "dihedral.h" #include "improper.h" +#include "special.h" #include "error.h" #include "memory.h" -#include "special.h" using namespace LAMMPS_NS; diff --git a/src/special.cpp b/src/special.cpp index eb7d8f826f..5bcb4265b5 100644 --- a/src/special.cpp +++ b/src/special.cpp @@ -17,11 +17,16 @@ #include "atom.h" #include "atom_vec.h" #include "force.h" +#include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; +// allocate space for static class variable + +Special *Special::sptr; + /* ---------------------------------------------------------------------- */ Special::Special(LAMMPS *lmp) : Pointers(lmp) @@ -49,10 +54,9 @@ Special::~Special() void Special::build() { - int i,j,k,m,n,loop,size,original; - int num12,num13,num14; - int max,maxall,messtag,nbuf,nbufmax; - int *buf,*bufcopy,*count; + int i,j,k,m,n,size; + int max,maxall,nbuf; + int *buf; MPI_Request request; MPI_Status status; @@ -100,10 +104,7 @@ void Special::build() nbuf = 0; for (i = 0; i < nlocal; i++) nbuf += num_bond[i]; - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"special:buf"); // fill buffer with global tags of bond partners of my atoms @@ -116,23 +117,10 @@ void Special::build() // when receive buffer, scan tags for atoms I own // when find one, increment nspecial count for that atom - messtag = 1; - for (loop = 0; loop < nprocs; loop++) { - for (i = 0; i < size; i++) { - m = atom->map(buf[i]); - if (m >= 0 && m < nlocal) nspecial[m][0]++; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + sptr = this; + comm->ring(size,sizeof(int),buf,1,ring_one,NULL); - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); } // ---------------------------------------------------- @@ -153,7 +141,7 @@ void Special::build() // count = accumulating counter - count = new int[nlocal]; + memory->create(count,nlocal,"special:count"); for (i = 0; i < nlocal; i++) count[i] = 0; // add bond partners stored by atom to onetwo list @@ -172,10 +160,7 @@ void Special::build() nbuf = 0; for (i = 0; i < nlocal; i++) nbuf += 2*num_bond[i]; - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"special:buf"); // fill buffer with global tags of both atoms in bond @@ -190,26 +175,13 @@ void Special::build() // when receive buffer, scan 2nd-atom tags for atoms I own // when find one, add 1st-atom tag to onetwo list for 2nd atom - messtag = 2; - for (loop = 0; loop < nprocs; loop++) { - for (i = 1; i < size; i += 2) { - m = atom->map(buf[i]); - if (m >= 0 && m < nlocal) onetwo[m][count[m]++] = buf[i-1]; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + sptr = this; + comm->ring(size,sizeof(int),buf,2,ring_two,NULL); - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); } - delete [] count; + memory->destroy(count); // ----------------------------------------------------- // done if special_bonds for 1-3, 1-4 are set to 1.0 @@ -230,10 +202,7 @@ void Special::build() nbuf = 0; for (i = 0; i < nlocal; i++) nbuf += 2 + nspecial[i][0]; - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"special:buf"); // fill buffer with: // (1) = counter for 1-3 neighbors, initialized to 0 @@ -252,27 +221,8 @@ void Special::build() // when find one, increment 1-3 count by # of 1-2 neighbors of my atom, // subtracting one since my list will contain original atom - messtag = 3; - for (loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - n = buf[i]; - num12 = buf[i+1]; - for (j = 0; j < num12; j++) { - m = atom->map(buf[i+2+j]); - if (m >= 0 && m < nlocal) n += nspecial[m][0] - 1; - } - buf[i] = n; - i += 2 + num12; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + sptr = this; + comm->ring(size,sizeof(int),buf,3,ring_three,buf); // extract count from buffer that has cycled back to me // nspecial[i][1] = # of 1-3 neighbors of atom i @@ -283,8 +233,7 @@ void Special::build() j += 2 + nspecial[i][0]; } - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); // ---------------------------------------------------- // create onethree[i] = list of 1-3 neighbors for atom i @@ -307,10 +256,7 @@ void Special::build() nbuf = 0; for (i = 0; i < nlocal; i++) nbuf += 4 + nspecial[i][0] + nspecial[i][1]; - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"special:buf"); // fill buffer with: // (1) = global tag of original atom @@ -337,32 +283,8 @@ void Special::build() // exclude the atom whose tag = original // this process may include duplicates but they will be culled later - messtag = 4; - for (loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - original = buf[i]; - num12 = buf[i+1]; - num13 = buf[i+2]; - n = buf[i+3]; - for (j = 0; j < num12; j++) { - m = atom->map(buf[i+4+j]); - if (m >= 0 && m < nlocal) - for (k = 0; k < nspecial[m][0]; k++) - if (onetwo[m][k] != original) - buf[i+4+num12+(n++)] = onetwo[m][k]; - } - buf[i+3] = n; - i += 4 + num12 + num13; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + sptr = this; + comm->ring(size,sizeof(int),buf,4,ring_four,buf); // fill onethree with buffer values that have been returned to me // sanity check: accumulated buf[i+3] count should equal @@ -377,8 +299,7 @@ void Special::build() onethree[i][k] = buf[j++]; } - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); // done if special_bonds for 1-4 are set to 1.0 @@ -397,10 +318,7 @@ void Special::build() nbuf = 0; for (i = 0; i < nlocal; i++) nbuf += 2 + nspecial[i][1]; - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"special:buf"); // fill buffer with: // (1) = counter for 1-4 neighbors, initialized to 0 @@ -419,27 +337,8 @@ void Special::build() // when find one, increment 1-4 count by # of 1-2 neighbors of my atom // may include duplicates and original atom but they will be culled later - messtag = 5; - for (loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - n = buf[i]; - num13 = buf[i+1]; - for (j = 0; j < num13; j++) { - m = atom->map(buf[i+2+j]); - if (m >= 0 && m < nlocal) n += nspecial[m][0]; - } - buf[i] = n; - i += 2 + num13; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + sptr = this; + comm->ring(size,sizeof(int),buf,5,ring_five,buf); // extract count from buffer that has cycled back to me // nspecial[i][2] = # of 1-4 neighbors of atom i @@ -450,8 +349,7 @@ void Special::build() j += 2 + nspecial[i][1]; } - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); // ---------------------------------------------------- // create onefour[i] = list of 1-4 neighbors for atom i @@ -475,10 +373,7 @@ void Special::build() nbuf = 0; for (i = 0; i < nlocal; i++) nbuf += 3 + nspecial[i][1] + nspecial[i][2]; - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - buf = new int[nbufmax]; - bufcopy = new int[nbufmax]; + memory->create(buf,nbuf,"special:buf"); // fill buffer with: // (1) = # of 1-3 neighbors @@ -502,30 +397,8 @@ void Special::build() // incrementing the count in buf(i+4) // this process may include duplicates but they will be culled later - messtag = 6; - for (loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - num13 = buf[i]; - num14 = buf[i+1]; - n = buf[i+2]; - for (j = 0; j < num13; j++) { - m = atom->map(buf[i+3+j]); - if (m >= 0 && m < nlocal) - for (k = 0; k < nspecial[m][0]; k++) - buf[i+3+num13+(n++)] = onetwo[m][k]; - } - buf[i+2] = n; - i += 3 + num13 + num14; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + sptr = this; + comm->ring(size,sizeof(int),buf,6,ring_six,buf); // fill onefour with buffer values that have been returned to me // sanity check: accumulated buf[i+2] count should equal @@ -540,8 +413,7 @@ void Special::build() onefour[i][k] = buf[j++]; } - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); combine(); if (force->special_angle) angle_trim(); @@ -699,7 +571,7 @@ void Special::combine() void Special::angle_trim() { - int i,j,m,n,iglobal,jglobal,ilocal,jlocal; + int i,j,m,n; MPI_Request request; MPI_Status status; @@ -743,7 +615,6 @@ void Special::angle_trim() int maxcount = 0; for (i = 0; i < nlocal; i++) maxcount = MAX(maxcount,nspecial[i][1]-nspecial[i][0]); - int **dflag; memory->create(dflag,nlocal,maxcount,"special::dflag"); for (i = 0; i < nlocal; i++) { @@ -757,12 +628,8 @@ void Special::angle_trim() int nbuf = 0; for (i = 0; i < nlocal; i++) nbuf += 2*num_angle[i] + 2*2*num_dihedral[i]; - - int nbufmax; - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - int *buf = new int[nbufmax]; - int *bufcopy = new int[nbufmax]; + int *buf; + memory->create(buf,nbuf,"special:buf"); // fill buffer with list of 1,3 atoms in each angle // and with list of 1,3 and 2,4 atoms in each dihedral @@ -785,41 +652,8 @@ void Special::angle_trim() // when receive buffer, scan list of 1,3 atoms looking for atoms I own // when find one, scan its 1-3 neigh list and mark I,J as in an angle - int next = me + 1; - int prev = me -1; - if (next == nprocs) next = 0; - if (prev < 0) prev = nprocs - 1; - - int messtag = 7; - for (int loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - iglobal = buf[i]; - jglobal = buf[i+1]; - ilocal = atom->map(iglobal); - jlocal = atom->map(jglobal); - if (ilocal >= 0 && ilocal < nlocal) - for (m = nspecial[ilocal][0]; m < nspecial[ilocal][1]; m++) - if (jglobal == special[ilocal][m]) { - dflag[ilocal][m-nspecial[ilocal][0]] = 1; - break; - } - if (jlocal >= 0 && jlocal < nlocal) - for (m = nspecial[jlocal][0]; m < nspecial[jlocal][1]; m++) - if (iglobal == special[jlocal][m]) { - dflag[jlocal][m-nspecial[jlocal][0]] = 1; - break; - } - i += 2; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + sptr = this; + comm->ring(size,sizeof(int),buf,7,ring_seven,NULL); // delete 1-3 neighbors if they are not flagged in dflag // preserve 1-4 neighbors @@ -839,8 +673,7 @@ void Special::angle_trim() // clean up memory->destroy(dflag); - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); // if no angles or dihedrals are defined, // delete all 1-3 neighs, preserving 1-4 neighs @@ -880,7 +713,7 @@ void Special::angle_trim() void Special::dihedral_trim() { - int i,j,m,n,iglobal,jglobal,ilocal,jlocal; + int i,j,m,n; MPI_Request request; MPI_Status status; @@ -918,7 +751,6 @@ void Special::dihedral_trim() int maxcount = 0; for (i = 0; i < nlocal; i++) maxcount = MAX(maxcount,nspecial[i][2]-nspecial[i][1]); - int **dflag; memory->create(dflag,nlocal,maxcount,"special::dflag"); for (i = 0; i < nlocal; i++) { @@ -931,12 +763,8 @@ void Special::dihedral_trim() int nbuf = 0; for (i = 0; i < nlocal; i++) nbuf += 2*num_dihedral[i]; - - int nbufmax; - MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); - - int *buf = new int[nbufmax]; - int *bufcopy = new int[nbufmax]; + int *buf; + memory->create(buf,nbuf,"special:buf"); // fill buffer with list of 1,4 atoms in each dihedral @@ -951,41 +779,8 @@ void Special::dihedral_trim() // when receive buffer, scan list of 1,4 atoms looking for atoms I own // when find one, scan its 1-4 neigh list and mark I,J as in a dihedral - int next = me + 1; - int prev = me -1; - if (next == nprocs) next = 0; - if (prev < 0) prev = nprocs - 1; - - int messtag = 8; - for (int loop = 0; loop < nprocs; loop++) { - i = 0; - while (i < size) { - iglobal = buf[i]; - jglobal = buf[i+1]; - ilocal = atom->map(iglobal); - jlocal = atom->map(jglobal); - if (ilocal >= 0 && ilocal < nlocal) - for (m = nspecial[ilocal][1]; m < nspecial[ilocal][2]; m++) - if (jglobal == special[ilocal][m]) { - dflag[ilocal][m-nspecial[ilocal][1]] = 1; - break; - } - if (jlocal >= 0 && jlocal < nlocal) - for (m = nspecial[jlocal][1]; m < nspecial[jlocal][2]; m++) - if (iglobal == special[jlocal][m]) { - dflag[jlocal][m-nspecial[jlocal][1]] = 1; - break; - } - i += 2; - } - if (me != next) { - MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); - MPI_Send(buf,size,MPI_INT,next,messtag,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_INT,&size); - for (j = 0; j < size; j++) buf[j] = bufcopy[j]; - } - } + sptr = this; + comm->ring(size,sizeof(int),buf,8,ring_eight,NULL); // delete 1-4 neighbors if they are not flagged in dflag @@ -1000,8 +795,7 @@ void Special::dihedral_trim() // clean up memory->destroy(dflag); - delete [] buf; - delete [] bufcopy; + memory->destroy(buf); // if no dihedrals are defined, delete all 1-4 neighs @@ -1024,3 +818,254 @@ void Special::dihedral_trim() " %g = # of 1-4 neighbors after dihedral trim\n",allcount); } } + +/* ---------------------------------------------------------------------- + when receive buffer, scan tags for atoms I own + when find one, increment nspecial count for that atom +------------------------------------------------------------------------- */ + +void Special::ring_one(int ndatum, char *cbuf) +{ + Atom *atom = sptr->atom; + int **nspecial = atom->nspecial; + int nlocal = atom->nlocal; + + int *buf = (int *) cbuf; + int m; + + for (int i = 0; i < ndatum; i++) { + m = atom->map(buf[i]); + if (m >= 0 && m < nlocal) nspecial[m][0]++; + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan 2nd-atom tags for atoms I own + when find one, add 1st-atom tag to onetwo list for 2nd atom +------------------------------------------------------------------------- */ + +void Special::ring_two(int ndatum, char *cbuf) +{ + Atom *atom = sptr->atom; + int nlocal = atom->nlocal; + + int **onetwo = sptr->onetwo; + int *count = sptr->count; + + int *buf = (int *) cbuf; + int m; + + for (int i = 1; i < ndatum; i += 2) { + m = atom->map(buf[i]); + if (m >= 0 && m < nlocal) onetwo[m][count[m]++] = buf[i-1]; + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan list of 1-2 neighbors for atoms I own + when find one, increment 1-3 count by # of 1-2 neighbors of my atom, + subtracting one since my list will contain original atom +------------------------------------------------------------------------- */ + +void Special::ring_three(int ndatum, char *cbuf) +{ + Atom *atom = sptr->atom; + int **nspecial = atom->nspecial; + int nlocal = atom->nlocal; + + int *buf = (int *) cbuf; + int i,j,m,n,num12; + + i = 0; + while (i < ndatum) { + n = buf[i]; + num12 = buf[i+1]; + for (j = 0; j < num12; j++) { + m = atom->map(buf[i+2+j]); + if (m >= 0 && m < nlocal) n += nspecial[m][0] - 1; + } + buf[i] = n; + i += 2 + num12; + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan list of 1-2 neighbors for atoms I own + when find one, add its neighbors to 1-3 list + increment the count in buf(i+4) + exclude the atom whose tag = original + this process may include duplicates but they will be culled later +------------------------------------------------------------------------- */ + +void Special::ring_four(int ndatum, char *cbuf) +{ + Atom *atom = sptr->atom; + int **nspecial = atom->nspecial; + int nlocal = atom->nlocal; + + int **onetwo = sptr->onetwo; + + int *buf = (int *) cbuf; + int i,j,k,m,n,original,num12,num13; + + i = 0; + while (i < ndatum) { + original = buf[i]; + num12 = buf[i+1]; + num13 = buf[i+2]; + n = buf[i+3]; + for (j = 0; j < num12; j++) { + m = atom->map(buf[i+4+j]); + if (m >= 0 && m < nlocal) + for (k = 0; k < nspecial[m][0]; k++) + if (onetwo[m][k] != original) + buf[i+4+num12+(n++)] = onetwo[m][k]; + } + buf[i+3] = n; + i += 4 + num12 + num13; + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan list of 1-3 neighbors for atoms I own + when find one, increment 1-4 count by # of 1-2 neighbors of my atom + may include duplicates and original atom but they will be culled later +------------------------------------------------------------------------- */ + +void Special::ring_five(int ndatum, char *cbuf) +{ + Atom *atom = sptr->atom; + int **nspecial = atom->nspecial; + int nlocal = atom->nlocal; + + int *buf = (int *) cbuf; + int i,j,m,n,num13; + + i = 0; + while (i < ndatum) { + n = buf[i]; + num13 = buf[i+1]; + for (j = 0; j < num13; j++) { + m = atom->map(buf[i+2+j]); + if (m >= 0 && m < nlocal) n += nspecial[m][0]; + } + buf[i] = n; + i += 2 + num13; + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan list of 1-3 neighbors for atoms I own + when find one, add its neighbors to 1-4 list + incrementing the count in buf(i+4) + this process may include duplicates but they will be culled later +------------------------------------------------------------------------- */ + +void Special::ring_six(int ndatum, char *cbuf) +{ + Atom *atom = sptr->atom; + int **nspecial = atom->nspecial; + int nlocal = atom->nlocal; + + int **onetwo = sptr->onetwo; + + int *buf = (int *) cbuf; + int i,j,k,m,n,num13,num14; + + i = 0; + while (i < ndatum) { + num13 = buf[i]; + num14 = buf[i+1]; + n = buf[i+2]; + for (j = 0; j < num13; j++) { + m = atom->map(buf[i+3+j]); + if (m >= 0 && m < nlocal) + for (k = 0; k < nspecial[m][0]; k++) + buf[i+3+num13+(n++)] = onetwo[m][k]; + } + buf[i+2] = n; + i += 3 + num13 + num14; + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan list of 1,3 atoms looking for atoms I own + when find one, scan its 1-3 neigh list and mark I,J as in an angle +------------------------------------------------------------------------- */ + +void Special::ring_seven(int ndatum, char *cbuf) +{ + Atom *atom = sptr->atom; + int **nspecial = atom->nspecial; + int **special = atom->special; + int nlocal = atom->nlocal; + + int **dflag = sptr->dflag; + + int *buf = (int *) cbuf; + int i,m,iglobal,jglobal,ilocal,jlocal; + + i = 0; + while (i < ndatum) { + iglobal = buf[i]; + jglobal = buf[i+1]; + ilocal = atom->map(iglobal); + jlocal = atom->map(jglobal); + if (ilocal >= 0 && ilocal < nlocal) + for (m = nspecial[ilocal][0]; m < nspecial[ilocal][1]; m++) + if (jglobal == special[ilocal][m]) { + dflag[ilocal][m-nspecial[ilocal][0]] = 1; + break; + } + if (jlocal >= 0 && jlocal < nlocal) + for (m = nspecial[jlocal][0]; m < nspecial[jlocal][1]; m++) + if (iglobal == special[jlocal][m]) { + dflag[jlocal][m-nspecial[jlocal][0]] = 1; + break; + } + i += 2; + } +} + +/* ---------------------------------------------------------------------- + when receive buffer, scan list of 1,4 atoms looking for atoms I own + when find one, scan its 1-4 neigh list and mark I,J as in a dihedral +------------------------------------------------------------------------- */ + +void Special::ring_eight(int ndatum, char *cbuf) +{ + Atom *atom = sptr->atom; + int **nspecial = atom->nspecial; + int **special = atom->special; + int nlocal = atom->nlocal; + + int **dflag = sptr->dflag; + + int *buf = (int *) cbuf; + int i,m,iglobal,jglobal,ilocal,jlocal; + + i = 0; + while (i < ndatum) { + iglobal = buf[i]; + jglobal = buf[i+1]; + ilocal = atom->map(iglobal); + jlocal = atom->map(jglobal); + if (ilocal >= 0 && ilocal < nlocal) + for (m = nspecial[ilocal][1]; m < nspecial[ilocal][2]; m++) + if (jglobal == special[ilocal][m]) { + dflag[ilocal][m-nspecial[ilocal][1]] = 1; + break; + } + if (jlocal >= 0 && jlocal < nlocal) + for (m = nspecial[jlocal][1]; m < nspecial[jlocal][2]; m++) + if (iglobal == special[jlocal][m]) { + dflag[jlocal][m-nspecial[jlocal][1]] = 1; + break; + } + i += 2; + } +} + + + + diff --git a/src/special.h b/src/special.h index b550f1ac47..6d07cf00dd 100644 --- a/src/special.h +++ b/src/special.h @@ -29,9 +29,27 @@ class Special : protected Pointers { int **onetwo,**onethree,**onefour; int dihedral_flag; + // data used by ring callback methods + + int *count; + int **dflag; + void combine(); void angle_trim(); void dihedral_trim(); + + // static variable for ring communication callback to access class data + // callback functions for ring communication + + static Special *sptr; + static void ring_one(int, char *); + static void ring_two(int, char *); + static void ring_three(int, char *); + static void ring_four(int, char *); + static void ring_five(int, char *); + static void ring_six(int, char *); + static void ring_seven(int, char *); + static void ring_eight(int, char *); }; }