diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index f9d181ef96..4de368ca85 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -414,9 +414,9 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : if (comm->me == 0) { if (screen) - fprintf(screen," create_bodies CPU = %g secs\n",time2-time1); + fprintf(screen," create bodies CPU = %g secs\n",time2-time1); if (logfile) - fprintf(logfile," create_bodies CPU = %g secs\n",time2-time1); + fprintf(logfile," create bodies CPU = %g secs\n",time2-time1); } // set nlocal_body and allocate bodies I own @@ -1749,6 +1749,8 @@ int FixRigidSmall::rendezvous_body(int n, char *inbuf, memory->destroy(iclose); memory->destroy(rsqclose); + // flag = 2: new outbuf + rflag = 2; return nout; } diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index e0d1bf132b..66c92d42c5 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -219,8 +219,19 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : // identify all SHAKE clusters + double time1 = MPI_Wtime(); + find_clusters(); + double time2 = MPI_Wtime(); + + if (comm->me == 0) { + if (screen) + fprintf(screen," find clusters CPU = %g secs\n",time2-time1); + if (logfile) + fprintf(logfile," find clusters CPU = %g secs\n",time2-time1); + } + // initialize list of SHAKE clusters to constrain maxlist = 0; @@ -707,13 +718,6 @@ void FixShake::find_clusters() int nlocal = atom->nlocal; int angles_allow = atom->avec->angles_allow; - // setup ring of procs - - int next = me + 1; - int prev = me -1; - if (next == nprocs) next = 0; - if (prev < 0) prev = nprocs - 1; - // ----------------------------------------------------- // allocate arrays for self (1d) and bond partners (2d) // max = max # of bond partners for owned atoms = 2nd dim of partner arrays @@ -755,6 +759,10 @@ void FixShake::find_clusters() memory->create(partner_shake,nlocal,max,"shake:partner_shake"); memory->create(partner_nshake,nlocal,max,"shake:partner_nshake"); + // setup atomIDs and procowner vectors in rendezvous decomposition + + atom_owners(); + // ----------------------------------------------------- // set npartner and partner_tag from special arrays // ----------------------------------------------------- @@ -778,86 +786,13 @@ void FixShake::find_clusters() } // ----------------------------------------------------- - // set partner_mask, partner_type, partner_massflag, partner_bondtype - // for bonded partners - // requires communication for off-proc partners + // set partner_mask, partner_type, partner_massflag, + // partner_bondtype for all my bonded partners + // requires rendezvous communication for off-proc partners // ----------------------------------------------------- - // fill in mask, type, massflag, bondtype if own bond partner - // info to store in buf for each off-proc bond = nper = 6 - // 2 atoms IDs in bond, space for mask, type, massflag, bondtype - // nbufmax = largest buffer needed to hold info from any proc - - int nper = 6; - - nbuf = 0; - for (i = 0; i < nlocal; i++) { - for (j = 0; j < npartner[i]; j++) { - partner_mask[i][j] = 0; - partner_type[i][j] = 0; - partner_massflag[i][j] = 0; - partner_bondtype[i][j] = 0; - - m = atom->map(partner_tag[i][j]); - if (m >= 0 && m < nlocal) { - partner_mask[i][j] = mask[m]; - partner_type[i][j] = type[m]; - if (nmass) { - if (rmass) massone = rmass[m]; - else massone = mass[type[m]]; - partner_massflag[i][j] = masscheck(massone); - } - n = bondtype_findset(i,tag[i],partner_tag[i][j],0); - if (n) partner_bondtype[i][j] = n; - else { - n = bondtype_findset(m,tag[i],partner_tag[i][j],0); - if (n) partner_bondtype[i][j] = n; - } - } else nbuf += nper; - } - } - - memory->create(buf,nbuf,"shake:buf"); - - // fill buffer with info - - size = 0; - for (i = 0; i < nlocal; i++) { - for (j = 0; j < npartner[i]; j++) { - m = atom->map(partner_tag[i][j]); - if (m < 0 || m >= nlocal) { - buf[size] = tag[i]; - buf[size+1] = partner_tag[i][j]; - buf[size+2] = 0; - buf[size+3] = 0; - buf[size+4] = 0; - n = bondtype_findset(i,tag[i],partner_tag[i][j],0); - if (n) buf[size+5] = n; - else buf[size+5] = 0; - size += nper; - } - } - } - - // cycle buffer around ring of procs back to self - - comm->ring(size,sizeof(tagint),buf,1,ring_bonds,buf,(void *)this); - - // store partner info returned to me - - m = 0; - while (m < size) { - i = atom->map(buf[m]); - for (j = 0; j < npartner[i]; j++) - if (buf[m+1] == partner_tag[i][j]) break; - partner_mask[i][j] = buf[m+2]; - partner_type[i][j] = buf[m+3]; - partner_massflag[i][j] = buf[m+4]; - partner_bondtype[i][j] = buf[m+5]; - m += nper; - } - - memory->destroy(buf); + partner_info(npartner,partner_tag,partner_mask,partner_type, + partner_massflag,partner_bondtype); // error check for unfilled partner info // if partner_type not set, is an error @@ -868,17 +803,18 @@ void FixShake::find_clusters() // else it's an error flag = 0; + int flag2 = 0; for (i = 0; i < nlocal; i++) for (j = 0; j < npartner[i]; j++) { - if (partner_type[i][j] == 0) flag = 1; + if (partner_type[i][j] == 0) flag++; if (!(mask[i] & groupbit)) continue; if (!(partner_mask[i][j] & groupbit)) continue; - if (partner_bondtype[i][j] == 0) flag = 1; + if (partner_bondtype[i][j] == 0) flag2++; } MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all(FLERR,"Did not find fix shake partner info"); - + // ----------------------------------------------------- // identify SHAKEable bonds // set nshake[i] = # of SHAKE bonds attached to atom i @@ -931,56 +867,11 @@ void FixShake::find_clusters() // ----------------------------------------------------- // set partner_nshake for bonded partners - // requires communication for off-proc partners + // requires rendezvous communication for off-proc partners // ----------------------------------------------------- - // fill in partner_nshake if own bond partner - // info to store in buf for each off-proc bond = - // 2 atoms IDs in bond, space for nshake value - // nbufmax = largest buffer needed to hold info from any proc - - nbuf = 0; - for (i = 0; i < nlocal; i++) { - for (j = 0; j < npartner[i]; j++) { - m = atom->map(partner_tag[i][j]); - if (m >= 0 && m < nlocal) partner_nshake[i][j] = nshake[m]; - else nbuf += 3; - } - } - - memory->create(buf,nbuf,"shake:buf"); - - // fill buffer with info - - size = 0; - for (i = 0; i < nlocal; i++) { - for (j = 0; j < npartner[i]; j++) { - m = atom->map(partner_tag[i][j]); - if (m < 0 || m >= nlocal) { - buf[size] = tag[i]; - buf[size+1] = partner_tag[i][j]; - size += 3; - } - } - } - - // cycle buffer around ring of procs back to self - - comm->ring(size,sizeof(tagint),buf,2,ring_nshake,buf,(void *)this); - - // store partner info returned to me - - m = 0; - while (m < size) { - i = atom->map(buf[m]); - for (j = 0; j < npartner[i]; j++) - if (buf[m+1] == partner_tag[i][j]) break; - partner_nshake[i][j] = buf[m+2]; - m += 3; - } - - memory->destroy(buf); - + nshake_info(npartner,partner_tag,partner_nshake); + // ----------------------------------------------------- // error checks // no atom with nshake > 3 @@ -988,7 +879,7 @@ void FixShake::find_clusters() // ----------------------------------------------------- flag = 0; - for (i = 0; i < nlocal; i++) if (nshake[i] > 3) flag = 1; + for (i = 0; i < nlocal; i++) if (nshake[i] > 3) flag++; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all(FLERR,"Shake cluster of more than 4 atoms"); @@ -996,7 +887,7 @@ void FixShake::find_clusters() for (i = 0; i < nlocal; i++) { if (nshake[i] <= 1) continue; for (j = 0; j < npartner[i]; j++) - if (partner_shake[i][j] && partner_nshake[i][j] > 1) flag = 1; + if (partner_shake[i][j] && partner_nshake[i][j] > 1) flag++; } MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all(FLERR,"Shake clusters are connected"); @@ -1067,60 +958,7 @@ void FixShake::find_clusters() // requires communication for off-proc atoms // ----------------------------------------------------- - // fill in shake arrays for each bond partner I own - // info to store in buf for each off-proc bond = - // all values from shake_flag, shake_atom, shake_type - // nbufmax = largest buffer needed to hold info from any proc - - nbuf = 0; - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; - for (j = 0; j < npartner[i]; j++) { - if (partner_shake[i][j] == 0) continue; - m = atom->map(partner_tag[i][j]); - if (m >= 0 && m < nlocal) { - shake_flag[m] = shake_flag[i]; - shake_atom[m][0] = shake_atom[i][0]; - shake_atom[m][1] = shake_atom[i][1]; - shake_atom[m][2] = shake_atom[i][2]; - shake_atom[m][3] = shake_atom[i][3]; - shake_type[m][0] = shake_type[i][0]; - shake_type[m][1] = shake_type[i][1]; - shake_type[m][2] = shake_type[i][2]; - } else nbuf += 9; - } - } - - memory->create(buf,nbuf,"shake:buf"); - - // fill buffer with info - - size = 0; - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; - for (j = 0; j < npartner[i]; j++) { - if (partner_shake[i][j] == 0) continue; - m = atom->map(partner_tag[i][j]); - if (m < 0 || m >= nlocal) { - buf[size] = partner_tag[i][j]; - buf[size+1] = shake_flag[i]; - buf[size+2] = shake_atom[i][0]; - buf[size+3] = shake_atom[i][1]; - buf[size+4] = shake_atom[i][2]; - buf[size+5] = shake_atom[i][3]; - buf[size+6] = shake_type[i][0]; - buf[size+7] = shake_type[i][1]; - buf[size+8] = shake_type[i][2]; - size += 9; - } - } - } - - // cycle buffer around ring of procs back to self - - comm->ring(size,sizeof(tagint),buf,3,ring_shake,NULL,(void *)this); - - memory->destroy(buf); + shake_info(npartner,partner_tag,partner_shake); // ----------------------------------------------------- // free local memory @@ -1199,98 +1037,549 @@ 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 + setup atomIDs and procowner ------------------------------------------------------------------------- */ -void FixShake::ring_bonds(int ndatum, char *cbuf, void *ptr) +void FixShake::atom_owners() { - FixShake *fsptr = (FixShake *)ptr; - Atom *atom = fsptr->atom; + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + + int *proclist; + memory->create(proclist,nlocal,"shake:proclist"); + IDRvous *idbuf = (IDRvous *) + memory->smalloc((bigint) nlocal*sizeof(IDRvous),"shake:idbuf"); + + // setup input buf to rendezvous comm + // input datums = pairs of bonded atoms + // owning proc for each datum = random hash of atomID + // one datum for each owned atom: datum = owning proc, atomID + + for (int i = 0; i < nlocal; i++) { + proclist[i] = tag[i] % nprocs; + idbuf[i].me = me; + idbuf[i].atomID = tag[i]; + } + + // perform rendezvous operation + // each proc assigned every 1/Pth atom + + char *buf; + comm->rendezvous(nlocal,proclist, + (char *) idbuf,sizeof(IDRvous), + rendezvous_ids,buf,0,(void *) this); + + memory->destroy(proclist); + memory->sfree(idbuf); +} + +/* ---------------------------------------------------------------------- + setup partner_mask, partner_type, partner_massflag, partner_bondtype +------------------------------------------------------------------------- */ + +void FixShake::partner_info(int *npartner, tagint **partner_tag, + int **partner_mask, int **partner_type, + int **partner_massflag, int **partner_bondtype) +{ + int i,j,m,n; + int nlocal = atom->nlocal; + + // nsend = # of my datums to send + // one datum for every off-processor partner + + int nsend = 0; + for (i = 0; i < nlocal; i++) { + for (j = 0; j < npartner[i]; j++) { + m = atom->map(partner_tag[i][j]); + if (m < 0 || m >= nlocal) nsend++; + } + } + + int *proclist; + memory->create(proclist,nsend,"special:proclist"); + PartnerInfo *inbuf = (PartnerInfo *) + memory->smalloc((bigint) nsend*sizeof(PartnerInfo),"special:inbuf"); + + // set values in 4 partner arrays for all partner atoms I own + // also setup input buf to rendezvous comm + // input datums = pair of bonded atoms where I do not own partner + // owning proc for each datum = partner_tag % nprocs + // datum: atomID = partner_tag (off-proc), partnerID = tag (on-proc) + // 4 values for my owned atom + double *rmass = atom->rmass; double *mass = atom->mass; - int *mask = atom->mask; int *type = atom->type; + int *mask = atom->mask; + tagint *tag = atom->tag; + + double massone; + + nsend = 0; + for (i = 0; i < nlocal; i++) { + for (j = 0; j < npartner[i]; j++) { + partner_mask[i][j] = 0; + partner_type[i][j] = 0; + partner_massflag[i][j] = 0; + partner_bondtype[i][j] = 0; + + m = atom->map(partner_tag[i][j]); + + if (m >= 0 && m < nlocal) { + partner_mask[i][j] = mask[m]; + partner_type[i][j] = type[m]; + if (nmass) { + if (rmass) massone = rmass[m]; + else massone = mass[type[m]]; + partner_massflag[i][j] = masscheck(massone); + } + n = bondtype_findset(i,tag[i],partner_tag[i][j],0); + if (n) partner_bondtype[i][j] = n; + else { + n = bondtype_findset(m,tag[i],partner_tag[i][j],0); + if (n) partner_bondtype[i][j] = n; + } + + } else { + proclist[nsend] = partner_tag[i][j] % nprocs; + inbuf[nsend].atomID = partner_tag[i][j]; + inbuf[nsend].partnerID = tag[i]; + inbuf[nsend].mask = mask[i]; + inbuf[nsend].type = type[i]; + if (nmass) { + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + inbuf[nsend].massflag = masscheck(massone); + } else inbuf[nsend].massflag = 0; + + // my atom may own bond, in which case set partner_bondtype + // else receiver of this datum will own the bond and return the value + + n = bondtype_findset(i,tag[i],partner_tag[i][j],0); + if (n) { + partner_bondtype[i][j] = n; + inbuf[nsend].bondtype = n; + } else inbuf[nsend].bondtype = 0; + + nsend++; + } + } + } + + // perform rendezvous operation + // each proc owns random subset of atoms + // receives all data needed to populate un-owned partner 4 values + + char *buf; + int nreturn = comm->rendezvous(nsend,proclist, + (char *) inbuf,sizeof(PartnerInfo), + rendezvous_partners_info,buf,sizeof(PartnerInfo), + (void *) this); + PartnerInfo *outbuf = (PartnerInfo *) buf; + + memory->destroy(proclist); + memory->sfree(inbuf); + + // set partner 4 values for un-onwed partners based on output info + // outbuf.atomID = my owned atom, outbuf.partnerID = partner the info is for + + for (m = 0; m < nreturn; m++) { + i = atom->map(outbuf[m].atomID); + for (j = 0; j < npartner[i]; j++) + if (partner_tag[i][j] == outbuf[m].partnerID) break; + partner_mask[i][j] = outbuf[m].mask; + partner_type[i][j] = outbuf[m].type; + partner_massflag[i][j] = outbuf[m].massflag; + + // only set partner_bondtype if my atom did not set it when setting up rendezvous + // if this proc set it, then sender of this datum set outbuf.bondtype = 0 + + if (partner_bondtype[i][j] == 0) + partner_bondtype[i][j] = outbuf[m].bondtype; + } + + memory->sfree(outbuf); +} + +/* ---------------------------------------------------------------------- + setup partner_nshake +------------------------------------------------------------------------- */ + +void FixShake::nshake_info(int *npartner, tagint **partner_tag, + int **partner_nshake) +{ + int i,j,m,n; int nlocal = atom->nlocal; + + // nsend = # of my datums to send + // one datum for every off-processor partner + + int nsend = 0; + for (i = 0; i < nlocal; i++) { + for (j = 0; j < npartner[i]; j++) { + m = atom->map(partner_tag[i][j]); + if (m < 0 || m >= nlocal) nsend++; + } + } + + int *proclist; + memory->create(proclist,nsend,"special:proclist"); + NShakeInfo *inbuf = (NShakeInfo *) + memory->smalloc((bigint) nsend*sizeof(NShakeInfo),"special:inbuf"); + + // set partner_nshake for all partner atoms I own + // also setup input buf to rendezvous comm + // input datums = pair of bonded atoms where I do not own partner + // owning proc for each datum = partner_tag % nprocs + // datum: atomID = partner_tag (off-proc), partnerID = tag (on-proc) + // nshake value for my owned atom + + tagint *tag = atom->tag; + + nsend = 0; + for (i = 0; i < nlocal; i++) { + for (j = 0; j < npartner[i]; j++) { + partner_nshake[i][j] = 0; + m = atom->map(partner_tag[i][j]); + if (m >= 0 && m < nlocal) { + partner_nshake[i][j] = nshake[m]; + } else { + proclist[nsend] = partner_tag[i][j] % nprocs; + inbuf[nsend].atomID = partner_tag[i][j]; + inbuf[nsend].partnerID = tag[i]; + inbuf[nsend].nshake = nshake[i]; + nsend++; + } + } + } + + // perform rendezvous operation + // each proc owns random subset of atoms + // receives all data needed to populate un-owned partner nshake + + char *buf; + int nreturn = comm->rendezvous(nsend,proclist, + (char *) inbuf,sizeof(NShakeInfo), + rendezvous_nshake,buf,sizeof(NShakeInfo), + (void *) this); + NShakeInfo *outbuf = (NShakeInfo *) buf; + + memory->destroy(proclist); + memory->sfree(inbuf); + + // set partner nshake for un-onwed partners based on output info + // outbuf.atomID = my owned atom, outbuf.partnerID = partner the info is for + + for (m = 0; m < nreturn; m++) { + i = atom->map(outbuf[m].atomID); + for (j = 0; j < npartner[i]; j++) + if (partner_tag[i][j] == outbuf[m].partnerID) break; + partner_nshake[i][j] = outbuf[m].nshake; + } + + memory->sfree(outbuf); +} + +/* ---------------------------------------------------------------------- + setup shake_flag, shake_atom, shake_type +------------------------------------------------------------------------- */ + +void FixShake::shake_info(int *npartner, tagint **partner_tag, + int **partner_shake) +{ + int i,j,m,n; + int nlocal = atom->nlocal; + + // nsend = # of my datums to send + // one datum for every off-processor partner + + int nsend = 0; + for (i = 0; i < nlocal; i++) { + for (j = 0; j < npartner[i]; j++) { + m = atom->map(partner_tag[i][j]); + if (m < 0 || m >= nlocal) nsend++; + } + } + + int *proclist; + memory->create(proclist,nsend,"special:proclist"); + ShakeInfo *inbuf = (ShakeInfo *) + memory->smalloc((bigint) nsend*sizeof(ShakeInfo),"special:inbuf"); + + // set 3 shake arrays for all partner atoms I own + // also setup input buf to rendezvous comm + // input datums = partner atom where I do not own partner + // owning proc for each datum = partner_tag % nprocs + // datum: atomID = partner_tag (off-proc) + // values in 3 shake arrays + + nsend = 0; + for (i = 0; i < nlocal; i++) { + if (shake_flag[i] == 0) continue; + for (j = 0; j < npartner[i]; j++) { + if (partner_shake[i][j] == 0) continue; + m = atom->map(partner_tag[i][j]); + + if (m >= 0 && m < nlocal) { + shake_flag[m] = shake_flag[i]; + shake_atom[m][0] = shake_atom[i][0]; + shake_atom[m][1] = shake_atom[i][1]; + shake_atom[m][2] = shake_atom[i][2]; + shake_atom[m][3] = shake_atom[i][3]; + shake_type[m][0] = shake_type[i][0]; + shake_type[m][1] = shake_type[i][1]; + shake_type[m][2] = shake_type[i][2]; + + } else { + proclist[nsend] = partner_tag[i][j] % nprocs; + inbuf[nsend].atomID = partner_tag[i][j]; + inbuf[nsend].shake_flag = shake_flag[i]; + inbuf[nsend].shake_atom[0] = shake_atom[i][0]; + inbuf[nsend].shake_atom[1] = shake_atom[i][1]; + inbuf[nsend].shake_atom[2] = shake_atom[i][2]; + inbuf[nsend].shake_atom[3] = shake_atom[i][3]; + inbuf[nsend].shake_type[0] = shake_type[i][0]; + inbuf[nsend].shake_type[1] = shake_type[i][1]; + inbuf[nsend].shake_type[2] = shake_type[i][2]; + nsend++; + } + } + } + + // perform rendezvous operation + // each proc owns random subset of atoms + // receives all data needed to populate un-owned shake info + + char *buf; + int nreturn = comm->rendezvous(nsend,proclist, + (char *) inbuf,sizeof(ShakeInfo), + rendezvous_shake,buf,sizeof(ShakeInfo), + (void *) this); + ShakeInfo *outbuf = (ShakeInfo *) buf; + + memory->destroy(proclist); + memory->sfree(inbuf); + + // set shake info for un-onwed partners based on output info + + for (m = 0; m < nreturn; m++) { + i = atom->map(outbuf[m].atomID); + shake_flag[i] = outbuf[m].shake_flag; + shake_atom[i][0] = outbuf[m].shake_atom[0]; + shake_atom[i][1] = outbuf[m].shake_atom[1]; + shake_atom[i][2] = outbuf[m].shake_atom[2]; + shake_atom[i][3] = outbuf[m].shake_atom[3]; + shake_type[i][0] = outbuf[m].shake_type[0]; + shake_type[i][1] = outbuf[m].shake_type[1]; + shake_type[i][2] = outbuf[m].shake_type[2]; + } + + memory->sfree(outbuf); +} + +/* ---------------------------------------------------------------------- + process data for atoms assigned to me in rendezvous decomposition + inbuf = list of N IDRvous datums + no outbuf +------------------------------------------------------------------------- */ + +int FixShake::rendezvous_ids(int n, char *inbuf, + int &flag, int *&proclist, char *&outbuf, + void *ptr) +{ + FixShake *fsptr = (FixShake *) ptr; + Memory *memory = fsptr->memory; + + int *procowner; + tagint *atomIDs; + + memory->create(procowner,n,"special:procowner"); + memory->create(atomIDs,n,"special:atomIDs"); + + IDRvous *in = (IDRvous *) inbuf; + + for (int i = 0; i < n; i++) { + procowner[i] = in[i].me; + atomIDs[i] = in[i].atomID; + } + + // store rendezvous data in FixShake class + + fsptr->nrvous = n; + fsptr->procowner = procowner; + fsptr->atomIDs = atomIDs; + + // flag = 0: no 2nd irregular comm needed in comm->rendezvous + + flag = 0; + return 0; +} + +/* ---------------------------------------------------------------------- + process data for atoms assigned to me in rendezvous decomposition + inbuf = list of N PairRvous datums + outbuf = same list of N PairRvous datums, routed to different procs +------------------------------------------------------------------------- */ + +int FixShake::rendezvous_partners_info(int n, char *inbuf, + int &flag, int *&proclist, char *&outbuf, + void *ptr) +{ + int i,m; + + FixShake *fsptr = (FixShake *) ptr; + Atom *atom = fsptr->atom; + Memory *memory = fsptr->memory; + + // clear atom map so it can be here as a hash table + // faster than an STL map for large atom counts + + atom->map_clear(); + + // hash atom IDs stored in rendezvous decomposition + + int nrvous = fsptr->nrvous; + tagint *atomIDs = fsptr->atomIDs; + + for (i = 0; i < nrvous; i++) + atom->map_one(atomIDs[i],i); + + // proclist = owner of atomID in caller decomposition + // outbuf = info about owned atomID = 4 values + + PartnerInfo *in = (PartnerInfo *) inbuf; + int *procowner = fsptr->procowner; + memory->create(proclist,n,"shake:proclist"); + + double massone; int nmass = fsptr->nmass; - tagint *buf = (tagint *) 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->bondtype_findset(m,buf[i],buf[i+1],0); - if (n) buf[i+5] = n; - } - } + for (i = 0; i < n; i++) { + m = atom->map(in[i].atomID); + proclist[i] = procowner[m]; } + + outbuf = inbuf; + + // re-create atom map + + atom->map_init(0); + atom->nghost = 0; + atom->map_set(); + + // flag = 1: outbuf = inbuf + + flag = 1; + return n; } /* ---------------------------------------------------------------------- - when receive buffer, scan bond partner IDs for atoms I own - if I own partner, fill in nshake value + process data for atoms assigned to me in rendezvous decomposition + inbuf = list of N NShakeInfo datums + outbuf = same list of N NShakeInfo datums, routed to different procs ------------------------------------------------------------------------- */ -void FixShake::ring_nshake(int ndatum, char *cbuf, void *ptr) +int FixShake::rendezvous_nshake(int n, char *inbuf, + int &flag, int *&proclist, char *&outbuf, + void *ptr) { - FixShake *fsptr = (FixShake *)ptr; + int i,j,m; + + FixShake *fsptr = (FixShake *) ptr; Atom *atom = fsptr->atom; - int nlocal = atom->nlocal; + Memory *memory = fsptr->memory; - int *nshake = fsptr->nshake; + // clear atom map so it can be here as a hash table + // faster than an STL map for large atom counts - tagint *buf = (tagint *) cbuf; - int m; + atom->map_clear(); - for (int i = 0; i < ndatum; i += 3) { - m = atom->map(buf[i+1]); - if (m >= 0 && m < nlocal) buf[i+2] = nshake[m]; + // hash atom IDs stored in rendezvous decomposition + + int nrvous = fsptr->nrvous; + tagint *atomIDs = fsptr->atomIDs; + + for (i = 0; i < nrvous; i++) + atom->map_one(atomIDs[i],i); + + // proclist = owner of atomID in caller decomposition + // outbuf = info about owned atomID + + NShakeInfo *in = (NShakeInfo *) inbuf; + int *procowner = fsptr->procowner; + memory->create(proclist,n,"shake:proclist"); + + for (i = 0; i < n; i++) { + m = atom->map(in[i].atomID); + proclist[i] = procowner[m]; } + + outbuf = inbuf; + + // re-create atom map + + atom->map_init(0); + atom->nghost = 0; + atom->map_set(); + + // flag = 1: outbuf = inbuf + + flag = 1; + return n; } - /* ---------------------------------------------------------------------- - when receive buffer, scan bond partner IDs for atoms I own - if I own partner, fill in nshake value + process data for atoms assigned to me in rendezvous decomposition + inbuf = list of N PairRvous datums + outbuf = same list of N PairRvous datums, routed to different procs ------------------------------------------------------------------------- */ -void FixShake::ring_shake(int ndatum, char *cbuf, void *ptr) +int FixShake::rendezvous_shake(int n, char *inbuf, + int &flag, int *&proclist, char *&outbuf, + void *ptr) { - FixShake *fsptr = (FixShake *)ptr; + int i,j,m; + + FixShake *fsptr = (FixShake *) ptr; Atom *atom = fsptr->atom; - int nlocal = atom->nlocal; + Memory *memory = fsptr->memory; - int *shake_flag = fsptr->shake_flag; - tagint **shake_atom = fsptr->shake_atom; - int **shake_type = fsptr->shake_type; + // clear atom map so it can be here as a hash table + // faster than an STL map for large atom counts - tagint *buf = (tagint *) cbuf; - int m; + atom->map_clear(); - 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]; - } + // hash atom IDs stored in rendezvous decomposition + + int nrvous = fsptr->nrvous; + tagint *atomIDs = fsptr->atomIDs; + + for (i = 0; i < nrvous; i++) + atom->map_one(atomIDs[i],i); + + // proclist = owner of atomID in caller decomposition + // outbuf = info about owned atomID + + ShakeInfo *in = (ShakeInfo *) inbuf; + int *procowner = fsptr->procowner; + memory->create(proclist,n,"shake:proclist"); + + for (i = 0; i < n; i++) { + m = atom->map(in[i].atomID); + proclist[i] = procowner[m]; } + + outbuf = inbuf; + + // re-create atom map + + atom->map_init(0); + atom->nghost = 0; + atom->map_set(); + + // flag = 1: outbuf = inbuf; + + flag = 1; + return n; } /* ---------------------------------------------------------------------- diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index d4e7b85ec4..2baea90a4a 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -120,6 +120,11 @@ class FixShake : public Fix { int nmol; void find_clusters(); + void atom_owners(); + void partner_info(int *, tagint **, int **, int **, int **, int **); + void nshake_info(int *, tagint **, int **); + void shake_info(int *, tagint **, int **); + int masscheck(double); void unconstrained_update(); void unconstrained_update_respa(int); @@ -131,12 +136,40 @@ class FixShake : public Fix { int bondtype_findset(int, tagint, tagint, int); int angletype_findset(int, tagint, tagint, int); - // static variable for ring communication callback to access class data - // callback functions for ring communication + // data used by rendezvous callback methods - static void ring_bonds(int, char *, void *); - static void ring_nshake(int, char *, void *); - static void ring_shake(int, char *, void *); + int nrvous; + tagint *atomIDs; + int *procowner; + + struct IDRvous { + int me; + tagint atomID; + }; + + struct PartnerInfo { + tagint atomID,partnerID; + int mask,type,massflag,bondtype; + }; + + struct NShakeInfo { + tagint atomID,partnerID; + int nshake; + }; + + struct ShakeInfo { + tagint atomID; + int shake_flag; + int shake_atom[4]; + int shake_type[3]; + }; + + // callback functions for rendezvous communication + + static int rendezvous_ids(int, char *, int &, int *&, char *&, void *); + static int rendezvous_partners_info(int, char *, int &, int *&, char *&, void *); + static int rendezvous_nshake(int, char *, int &, int *&, char *&, void *); + static int rendezvous_shake(int, char *, int &, int *&, char *&, void *); }; } diff --git a/src/special.cpp b/src/special.cpp index 79d2f77e46..b0d5bc7dca 100644 --- a/src/special.cpp +++ b/src/special.cpp @@ -21,7 +21,6 @@ #include "modify.h" #include "fix.h" #include "accelerator_kokkos.h" -#include "hashlittle.h" #include "atom_masks.h" #include "memory.h" #include "error.h" @@ -177,25 +176,20 @@ void Special::atom_owners() // input datums = pairs of bonded atoms // owning proc for each datum = random hash of atomID // one datum for each owned atom: datum = owning proc, atomID - // one datum for each bond partner: datum = atomID, bond partner ID - // add inverted datum when netwon_bond on for (int i = 0; i < nlocal; i++) { - //proc = hashlittle(&tag[i],sizeof(tagint),0) % nprocs; proclist[i] = tag[i] % nprocs; idbuf[i].me = me; idbuf[i].atomID = tag[i]; } // perform rendezvous operation - // each proc owns random subset of atoms - // receives all info to form and return their onetwo lists + // each proc assigned every 1/Pth atom char *buf; comm->rendezvous(nlocal,proclist, - (char *) idbuf,sizeof(PairRvous), - rendezvous_ids,buf,sizeof(PairRvous), - (void *) this); + (char *) idbuf,sizeof(IDRvous), + rendezvous_ids,buf,0,(void *) this); memory->destroy(proclist); memory->sfree(idbuf); @@ -215,49 +209,45 @@ void Special::onetwo_build_newton() int **nspecial = atom->nspecial; int nlocal = atom->nlocal; - // ncount = # of my datums to send - // include nlocal datums with owner of each atom + // nsend = # of my datums to send - int ncount = 0; + int nsend = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < num_bond[i]; j++) { m = atom->map(bond_atom[i][j]); - if (m < 0 || m >= nlocal) ncount++; + if (m < 0 || m >= nlocal) nsend++; } } int *proclist; - memory->create(proclist,ncount,"special:proclist"); + memory->create(proclist,nsend,"special:proclist"); PairRvous *inbuf = (PairRvous *) - memory->smalloc((bigint) ncount*sizeof(PairRvous),"special:inbuf"); + memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf"); // setup input buf to rendezvous comm // input datums = pairs of bonded atoms - // owning proc for each datum = random hash of atomID - // one datum for each owned atom: datum = owning proc, atomID - // one datum for each bond partner: datum = atomID, bond partner ID - // add inverted datum when netwon_bond on + // owning proc for each datum = atomID % nprocs + // one datum for each bond partner: bond partner ID, atomID - ncount = 0; + nsend = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < num_bond[i]; j++) { m = atom->map(bond_atom[i][j]); if (m >= 0 && m < nlocal) continue; - proclist[ncount] = bond_atom[i][j] % nprocs; - inbuf[ncount].atomID = bond_atom[i][j]; - inbuf[ncount].partnerID = tag[i]; - ncount++; + proclist[nsend] = bond_atom[i][j] % nprocs; + inbuf[nsend].atomID = bond_atom[i][j]; + inbuf[nsend].partnerID = tag[i]; + nsend++; } } // perform rendezvous operation // each proc owns random subset of atoms - // receives all info to form and return their onetwo lists char *buf; - int nreturn = comm->rendezvous(ncount,proclist, + int nreturn = comm->rendezvous(nsend,proclist, (char *) inbuf,sizeof(PairRvous), - rendezvous_1234,buf,sizeof(PairRvous), + rendezvous_pairs,buf,sizeof(PairRvous), (void *) this); PairRvous *outbuf = (PairRvous *) buf; @@ -312,6 +302,28 @@ void Special::onetwo_build_newton() void Special::onetwo_build_newton_off() { + int i,j; + + int *num_bond = atom->num_bond; + tagint **bond_atom = atom->bond_atom; + int **nspecial = atom->nspecial; + int nlocal = atom->nlocal; + + int max = 0; + for (i = 0; i < nlocal; i++) + max = MAX(max,num_bond[i]); + + MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world); + memory->create(onetwo,nlocal,maxall,"special:onetwo"); + + // nsend = # of my datums to send + // include nlocal datums with owner of each atom + + for (i = 0; i < nlocal; i++) { + nspecial[i][0] = num_bond[i]; + for (j = 0; j < num_bond[i]; j++) + onetwo[i][j] = bond_atom[i][j]; + } } /* ---------------------------------------------------------------------- @@ -327,21 +339,20 @@ void Special::onethree_build() int **nspecial = atom->nspecial; int nlocal = atom->nlocal; - // ncount = # of my datums to send - // include nlocal datums with owner of each atom + // nsend = # of my datums to send - int ncount = 0; + int nsend = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < nspecial[i][0]; j++) { m = atom->map(onetwo[i][j]); - if (m < 0 || m >= nlocal) ncount += nspecial[i][0]-1; + if (m < 0 || m >= nlocal) nsend += nspecial[i][0]-1; } } int *proclist; - memory->create(proclist,ncount,"special:proclist"); + memory->create(proclist,nsend,"special:proclist"); PairRvous *inbuf = (PairRvous *) - memory->smalloc((bigint) ncount*sizeof(PairRvous),"special:inbuf"); + memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf"); // setup input buf to rendezvous comm // input datums = all pairs of onetwo atoms (they are 1-3 neighbors) @@ -349,7 +360,7 @@ void Special::onethree_build() // one datum for each owned atom: datum = owning proc, atomID // one datum for each onetwo pair: datum = atomID1, atomID2 - ncount = 0; + nsend = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < nspecial[i][0]; j++) { m = atom->map(onetwo[i][j]); @@ -357,10 +368,10 @@ void Special::onethree_build() proc = onetwo[i][j] % nprocs; for (k = 0; k < nspecial[i][0]; k++) { if (j == k) continue; - proclist[ncount] = proc; - inbuf[ncount].atomID = onetwo[i][j]; - inbuf[ncount].partnerID = onetwo[i][k]; - ncount++; + proclist[nsend] = proc; + inbuf[nsend].atomID = onetwo[i][j]; + inbuf[nsend].partnerID = onetwo[i][k]; + nsend++; } } } @@ -370,9 +381,9 @@ void Special::onethree_build() // receives all info to form and return their onethree lists char *buf; - int nreturn = comm->rendezvous(ncount,proclist, + int nreturn = comm->rendezvous(nsend,proclist, (char *) inbuf,sizeof(PairRvous), - rendezvous_1234,buf,sizeof(PairRvous), + rendezvous_pairs,buf,sizeof(PairRvous), (void *) this); PairRvous *outbuf = (PairRvous *) buf; @@ -434,21 +445,21 @@ void Special::onefour_build() int **nspecial = atom->nspecial; int nlocal = atom->nlocal; - // ncount = # of my datums to send + // nsend = # of my datums to send // include nlocal datums with owner of each atom - int ncount = 0; + int nsend = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < nspecial[i][1]; j++) { m = atom->map(onethree[i][j]); - if (m < 0 || m >= nlocal) ncount += nspecial[i][0]; + if (m < 0 || m >= nlocal) nsend += nspecial[i][0]; } } int *proclist; - memory->create(proclist,ncount,"special:proclist"); + memory->create(proclist,nsend,"special:proclist"); PairRvous *inbuf = (PairRvous *) - memory->smalloc((bigint) ncount*sizeof(PairRvous),"special:inbuf"); + memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf"); // setup input buf to rendezvous comm // input datums = all pairs of onethree and onetwo atoms (they're 1-4 neighbors) @@ -456,17 +467,17 @@ void Special::onefour_build() // one datum for each owned atom: datum = owning proc, atomID // one datum for each onethree/onetwo pair: datum = atomID1, atomID2 - ncount = 0; + nsend = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < nspecial[i][1]; j++) { m = atom->map(onethree[i][j]); if (m >= 0 && m < nlocal) continue; proc = onethree[i][j] % nprocs; for (k = 0; k < nspecial[i][0]; k++) { - proclist[ncount] = proc; - inbuf[ncount].atomID = onethree[i][j]; - inbuf[ncount].partnerID = onetwo[i][k]; - ncount++; + proclist[nsend] = proc; + inbuf[nsend].atomID = onethree[i][j]; + inbuf[nsend].partnerID = onetwo[i][k]; + nsend++; } } } @@ -476,9 +487,9 @@ void Special::onefour_build() // receives all info to form and return their onefour lists char *buf; - int nreturn = comm->rendezvous(ncount,proclist, + int nreturn = comm->rendezvous(nsend,proclist, (char *) inbuf,sizeof(PairRvous), - rendezvous_1234,buf,sizeof(PairRvous), + rendezvous_pairs,buf,sizeof(PairRvous), (void *) this); PairRvous *outbuf = (PairRvous *) buf; @@ -773,7 +784,7 @@ void Special::combine() void Special::angle_trim() { - int i,j,m,n,proc,index; + int i,j,k,m;; int *num_angle = atom->num_angle; int *num_dihedral = atom->num_dihedral; @@ -804,96 +815,89 @@ void Special::angle_trim() " %g = # of 1-3 neighbors before angle trim\n",allcount); } - // if angles or dihedrals are defined, - // flag each 1-3 neigh if it appears in an angle or dihedral + // if angles or dihedrals are defined + // rendezvous angle 1-3 and dihedral 1-3,2-4 pairs if ((num_angle && atom->nangles) || (num_dihedral && atom->ndihedrals)) { - // ncount = # of my datums to send in 3 parts for each owned atom - // proc owner, onethree list, angle end points - // angle end points are from angle list and 1-3 and 2-4 pairs in dihedrals - // latter is only for angles or dihedrlas where I own atom2 + // nsend = # of my datums to send + // latter is only for angles or dihedrlas where I own atom2 (newton bond off) - int ncount = nlocal; - for (i = 0; i < nlocal; i++) ncount += nspecial[i][1]; + int nsend = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < num_angle[i]; j++) { - index = atom->map(angle_atom2[i][j]); - if (index >= 0 && index < nlocal) ncount += 2; + if (tag[i] != angle_atom2[i][j]) continue; + m = atom->map(angle_atom1[i][j]); + if (m < 0 || m >= nlocal) nsend++; + m = atom->map(angle_atom3[i][j]); + if (m < 0 || m >= nlocal) nsend++; } for (j = 0; j < num_dihedral[i]; j++) { - index = atom->map(dihedral_atom2[i][j]); - if (index >= 0 && index < nlocal) ncount += 4; + if (tag[i] != dihedral_atom2[i][j]) continue; + m = atom->map(dihedral_atom1[i][j]); + if (m < 0 || m >= nlocal) nsend++; + m = atom->map(dihedral_atom3[i][j]); + if (m < 0 || m >= nlocal) nsend++; + m = atom->map(dihedral_atom4[i][j]); + if (m < 0 || m >= nlocal) nsend++; } } int *proclist; - memory->create(proclist,ncount,"special:proclist"); + memory->create(proclist,nsend,"special:proclist"); PairRvous *inbuf = (PairRvous *) - memory->smalloc((bigint) ncount*sizeof(PairRvous),"special:inbuf"); + memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf"); // setup input buf to rendezvous comm - // one datum for each owned atom: datum = proc, atomID - // sent to owner of atomID - // one datum for each 1-4 partner: datum = atomID, ID - // sent to owner of atomID - // two datums for each dihedral 1-4 endatoms : datum = atomID, ID - // sent to owner of atomID - m = 0; + nsend = 0; for (i = 0; i < nlocal; i++) { - for (j = 0; j < nspecial[i][1]; j++) { - proclist[m] = proc; - //inbuf[m].me = -1; - inbuf[m].atomID = tag[i]; - inbuf[m].partnerID = onethree[i][j]; - m++; - } - for (j = 0; j < num_angle[i]; j++) { - index = atom->map(angle_atom2[i][j]); - if (index < 0 || index >= nlocal) continue; + if (tag[i] != angle_atom2[i][j]) continue; - proclist[m] = hashlittle(&angle_atom1[i][j],sizeof(tagint),0) % nprocs; - //inbuf[m].me = -2; - inbuf[m].atomID = angle_atom1[i][j]; - inbuf[m].partnerID = angle_atom3[i][j]; - m++; + m = atom->map(angle_atom1[i][j]); + if (m < 0 || m >= nlocal) { + proclist[nsend] = angle_atom1[i][j] % nprocs; + inbuf[nsend].atomID = angle_atom1[i][j]; + inbuf[nsend].partnerID = angle_atom3[i][j]; + nsend++; + } - proclist[m] = hashlittle(&angle_atom3[i][j],sizeof(tagint),0) % nprocs; - //inbuf[m].me = -2; - inbuf[m].atomID = angle_atom3[i][j]; - inbuf[m].partnerID = angle_atom1[i][j]; - m++; + m = atom->map(angle_atom3[i][j]); + if (m < 0 || m >= nlocal) { + proclist[nsend] = angle_atom3[i][j] % nprocs; + inbuf[nsend].atomID = angle_atom3[i][j]; + inbuf[nsend].partnerID = angle_atom1[i][j]; + nsend++; + } } for (j = 0; j < num_dihedral[i]; j++) { - index = atom->map(dihedral_atom2[i][j]); - if (index < 0 || index >= nlocal) continue; + if (tag[i] != dihedral_atom2[i][j]) continue; - proclist[m] = hashlittle(&dihedral_atom1[i][j],sizeof(tagint),0) % nprocs; - //inbuf[m].me = -2; - inbuf[m].atomID = dihedral_atom1[i][j]; - inbuf[m].partnerID = dihedral_atom3[i][j]; - m++; + m = atom->map(dihedral_atom1[i][j]); + if (m < 0 || m >= nlocal) { + proclist[nsend] = dihedral_atom1[i][j] % nprocs; + inbuf[nsend].atomID = dihedral_atom1[i][j]; + inbuf[nsend].partnerID = dihedral_atom3[i][j]; + nsend++; + } - proclist[m] = hashlittle(&dihedral_atom2[i][j],sizeof(tagint),0) % nprocs; - //inbuf[m].me = -2; - inbuf[m].atomID = dihedral_atom2[i][j]; - inbuf[m].partnerID = dihedral_atom4[i][j]; - m++; + m = atom->map(dihedral_atom3[i][j]); + if (m < 0 || m >= nlocal) { + proclist[nsend] = dihedral_atom3[i][j] % nprocs; + inbuf[nsend].atomID = dihedral_atom3[i][j]; + inbuf[nsend].partnerID = dihedral_atom1[i][j]; + nsend++; + } - proclist[m] = hashlittle(&dihedral_atom3[i][j],sizeof(tagint),0) % nprocs; - //inbuf[m].me = -2; - inbuf[m].atomID = dihedral_atom3[i][j]; - inbuf[m].partnerID = dihedral_atom1[i][j]; - m++; - - proclist[m] = hashlittle(&dihedral_atom4[i][j],sizeof(tagint),0) % nprocs; - //inbuf[m].me = -2; - inbuf[m].atomID = dihedral_atom4[i][j]; - inbuf[m].partnerID = dihedral_atom2[i][j]; - m++; + m = atom->map(dihedral_atom4[i][j]); + if (m < 0 || m >= nlocal) { + proclist[nsend] = dihedral_atom4[i][j] % nprocs; + inbuf[nsend].atomID = dihedral_atom4[i][j]; + inbuf[nsend].partnerID = dihedral_atom2[i][j]; + nsend++; + } } } @@ -903,26 +907,112 @@ void Special::angle_trim() // when done: each atom has atom ID of owning atom of its body char *buf; - int nreturn = comm->rendezvous(ncount,proclist, + int nreturn = comm->rendezvous(nsend,proclist, (char *) inbuf,sizeof(PairRvous), - rendezvous_trim,buf,sizeof(PairRvous), + rendezvous_pairs,buf,sizeof(PairRvous), (void *) this); PairRvous *outbuf = (PairRvous *) buf; memory->destroy(proclist); memory->sfree(inbuf); - // reset nspecial[1] and onethree for all owned atoms based on output info + // flag all onethree atoms to keep - for (i = 0; i < nlocal; i++) nspecial[i][1] = 0; + int max = 0; + for (i = 0; i < nlocal; i++) + max = MAX(max,nspecial[i][1]); + MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world); + + int **flag; + memory->create(flag,nlocal,maxall,"special:flag"); + + for (i = 0; i < nlocal; i++) + for (j = 0; j < nspecial[i][1]; j++) + flag[i][j] = 0; + + // reset nspecial[1] and onethree for all owned atoms based on output info + + for (i = 0; i < nlocal; i++) { + for (j = 0; j < num_angle[i]; j++) { + if (tag[i] != angle_atom2[i][j]) continue; + + m = atom->map(angle_atom1[i][j]); + if (m >= 0 && m < nlocal) { + for (k = 0; k < nspecial[m][1]; k++) + if (onethree[m][k] == angle_atom3[i][j]) { + flag[m][k] = 1; + break; + } + } + + m = atom->map(angle_atom3[i][j]); + if (m >= 0 && m < nlocal) { + for (k = 0; k < nspecial[m][1]; k++) + if (onethree[m][k] == angle_atom1[i][j]) { + flag[m][k] = 1; + break; + } + } + } + + for (j = 0; j < num_dihedral[i]; j++) { + if (tag[i] != dihedral_atom2[i][j]) continue; + + m = atom->map(dihedral_atom1[i][j]); + if (m >= 0 && m < nlocal) { + for (k = 0; k < nspecial[m][1]; k++) + if (onethree[m][k] == dihedral_atom3[i][j]) { + flag[m][k] = 1; + break; + } + } + + m = atom->map(dihedral_atom3[i][j]); + if (m >= 0 && m < nlocal) { + for (k = 0; k < nspecial[m][1]; k++) + if (onethree[m][k] == dihedral_atom1[i][j]) { + flag[m][k] = 1; + break; + } + } + + m = atom->map(dihedral_atom4[i][j]); + if (m >= 0 && m < nlocal) { + for (k = 0; k < nspecial[m][1]; k++) + if (onethree[m][k] == dihedral_atom2[i][j]) { + flag[m][k] = 1; + break; + } + } + } + } for (m = 0; m < nreturn; m++) { i = atom->map(outbuf[m].atomID); - onethree[i][nspecial[i][1]++] = outbuf[m].partnerID; + for (k = 0; k < nspecial[i][1]; k++) + if (onethree[i][k] == outbuf[m].partnerID) { + flag[i][k] = 1; + break; + } } - + memory->destroy(outbuf); + // use flag values to compress onefour list for each atom + + for (i = 0; i < nlocal; i++) { + j = 0; + while (j < nspecial[i][1]) { + if (flag[i][j] == 0) { + onethree[i][j] = onethree[i][nspecial[i][1]-1]; + flag[i][j] = flag[i][nspecial[i][1]-1]; + nspecial[i][1]--; + } else j++; + } + } + + memory->destroy(flag); + // if no angles or dihedrals are defined, delete all 1-3 neighs } else { @@ -952,7 +1042,7 @@ void Special::angle_trim() void Special::dihedral_trim() { - int i,j,m,n,proc,index; + int i,j,k,m; int *num_dihedral = atom->num_dihedral; tagint **dihedral_atom1 = atom->dihedral_atom1; @@ -978,68 +1068,51 @@ void Special::dihedral_trim() " %g = # of 1-4 neighbors before dihedral trim\n",allcount); } - // if dihedrals are defined, rendezvous onefour list with dihedral 1-4 pairs + // if dihedrals are defined, rendezvous dihedral 1-4 pairs if (num_dihedral && atom->ndihedrals) { - // ncount = # of my datums to send in 3 parts for each owned atom - // onefour list, proc owner, dihedral end points - // latter is only for dihedrals where I own atom2 + // nsend = # of my datums to send + // latter is only for dihedrals where I own atom2 (newton bond off) - int ncount = nlocal; - for (i = 0; i < nlocal; i++) ncount += nspecial[i][2]; + int nsend = 0; for (i = 0; i < nlocal; i++) { for (j = 0; j < num_dihedral[i]; j++) { - index = atom->map(dihedral_atom2[i][j]); - if (index >= 0 && index < nlocal) ncount += 2; + if (tag[i] != dihedral_atom2[i][j]) continue; + m = atom->map(dihedral_atom1[i][j]); + if (m < 0 || m >= nlocal) nsend++; + m = atom->map(dihedral_atom4[i][j]); + if (m < 0 || m >= nlocal) nsend++; } } int *proclist; - memory->create(proclist,ncount,"special:proclist"); + memory->create(proclist,nsend,"special:proclist"); PairRvous *inbuf = (PairRvous *) - memory->smalloc((bigint) ncount*sizeof(PairRvous),"special:inbuf"); - + memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf"); + // setup input buf to rendezvous comm - // one datum for each owned atom: datum = proc, atomID - // sent to owner of atomID - // one datum for each 1-4 partner: datum = atomID, ID - // sent to owner of atomID - // two datums for each dihedral 1-4 endatoms : datum = atomID, ID - // sent to owner of atomID - m = 0; + nsend = 0; for (i = 0; i < nlocal; i++) { - proc = hashlittle(&tag[i],sizeof(tagint),0) % nprocs; - proclist[m] = proc; - //inbuf[m].me = me; - inbuf[m].atomID = tag[i]; - inbuf[m].partnerID = 0; - m++; - - for (j = 0; j < nspecial[i][2]; j++) { - proclist[m] = proc; - //inbuf[m].me = -1; - inbuf[m].atomID = tag[i]; - inbuf[m].partnerID = onefour[i][j]; - m++; - } - for (j = 0; j < num_dihedral[i]; j++) { - index = atom->map(dihedral_atom2[i][j]); - if (index < 0 || index >= nlocal) continue; + if (tag[i] != dihedral_atom2[i][j]) continue; - proclist[m] = hashlittle(&dihedral_atom1[i][j],sizeof(tagint),0) % nprocs; - //inbuf[m].me = -2; - inbuf[m].atomID = dihedral_atom1[i][j]; - inbuf[m].partnerID = dihedral_atom4[i][j]; - m++; + m = atom->map(dihedral_atom1[i][j]); + if (m < 0 || m >= nlocal) { + proclist[nsend] = dihedral_atom1[i][j] % nprocs; + inbuf[nsend].atomID = dihedral_atom1[i][j]; + inbuf[nsend].partnerID = dihedral_atom4[i][j]; + nsend++; + } - proclist[m] = hashlittle(&dihedral_atom4[i][j],sizeof(tagint),0) % nprocs; - //inbuf[m].me = -2; - inbuf[m].atomID = dihedral_atom4[i][j]; - inbuf[m].partnerID = dihedral_atom1[i][j]; - m++; + m = atom->map(dihedral_atom4[i][j]); + if (m < 0 || m >= nlocal) { + proclist[nsend] = dihedral_atom4[i][j] % nprocs; + inbuf[nsend].atomID = dihedral_atom4[i][j]; + inbuf[nsend].partnerID = dihedral_atom1[i][j]; + nsend++; + } } } @@ -1049,26 +1122,81 @@ void Special::dihedral_trim() // when done: each atom has atom ID of owning atom of its body char *buf; - int nreturn = comm->rendezvous(ncount,proclist, + int nreturn = comm->rendezvous(nsend,proclist, (char *) inbuf,sizeof(PairRvous), - rendezvous_trim,buf,sizeof(PairRvous), + rendezvous_pairs,buf,sizeof(PairRvous), (void *) this); PairRvous *outbuf = (PairRvous *) buf; memory->destroy(proclist); memory->sfree(inbuf); + // flag all onefour atoms to keep + + int max = 0; + for (i = 0; i < nlocal; i++) + max = MAX(max,nspecial[i][2]); + MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world); + + int **flag; + memory->create(flag,nlocal,maxall,"special:flag"); + + for (i = 0; i < nlocal; i++) + for (j = 0; j < nspecial[i][2]; j++) + flag[i][j] = 0; + // reset nspecial[2] and onefour for all owned atoms based on output info - for (i = 0; i < nlocal; i++) nspecial[i][2] = 0; + for (i = 0; i < nlocal; i++) { + for (j = 0; j < num_dihedral[i]; j++) { + if (tag[i] != dihedral_atom2[i][j]) continue; + + m = atom->map(dihedral_atom1[i][j]); + if (m >= 0 && m < nlocal) { + for (k = 0; k < nspecial[m][2]; k++) + if (onefour[m][k] == dihedral_atom4[i][j]) { + flag[m][k] = 1; + break; + } + } + + m = atom->map(dihedral_atom4[i][j]); + if (m >= 0 && m < nlocal) { + for (k = 0; k < nspecial[m][2]; k++) + if (onefour[m][k] == dihedral_atom1[i][j]) { + flag[m][k] = 1; + break; + } + } + } + } for (m = 0; m < nreturn; m++) { i = atom->map(outbuf[m].atomID); - onefour[i][nspecial[i][2]++] = outbuf[m].partnerID; + for (k = 0; k < nspecial[i][2]; k++) + if (onefour[i][k] == outbuf[m].partnerID) { + flag[i][k] = 1; + break; + } } memory->destroy(outbuf); + // use flag values to compress onefour list for each atom + + for (i = 0; i < nlocal; i++) { + j = 0; + while (j < nspecial[i][2]) { + if (flag[i][j] == 0) { + onefour[i][j] = onefour[i][nspecial[i][2]-1]; + flag[i][j] = flag[i][nspecial[i][2]-1]; + nspecial[i][2]--; + } else j++; + } + } + + memory->destroy(flag); + // if no dihedrals are defined, delete all 1-4 neighs } else { @@ -1093,8 +1221,8 @@ void Special::dihedral_trim() /* ---------------------------------------------------------------------- process data for atoms assigned to me in rendezvous decomposition - inbuf = list of N PairRvous datums - outbuf = empty + inbuf = list of N IDRvous datums + no outbuf ------------------------------------------------------------------------- */ int Special::rendezvous_ids(int n, char *inbuf, @@ -1109,7 +1237,6 @@ int Special::rendezvous_ids(int n, char *inbuf, memory->create(procowner,n,"special:procowner"); memory->create(atomIDs,n,"special:atomIDs"); - // NOTE: when to free these vectors IDRvous *in = (IDRvous *) inbuf; @@ -1120,13 +1247,12 @@ int Special::rendezvous_ids(int n, char *inbuf, // store rendezvous data in Special class - sptr->ncount = n; + sptr->nrvous = n; sptr->procowner = procowner; sptr->atomIDs = atomIDs; - proclist = NULL; - outbuf = NULL; - + // flag = 0: no 2nd irregular comm needed in comm->rendezvous + flag = 0; return 0; } @@ -1138,7 +1264,7 @@ int Special::rendezvous_ids(int n, char *inbuf, outbuf = same list of N PairRvous datums, routed to different procs ------------------------------------------------------------------------- */ -int Special::rendezvous_1234(int n, char *inbuf, +int Special::rendezvous_pairs(int n, char *inbuf, int &flag, int *&proclist, char *&outbuf, void *ptr) { @@ -1153,10 +1279,10 @@ int Special::rendezvous_1234(int n, char *inbuf, // hash atom IDs stored in rendezvous decomposition - int ncount = sptr->ncount; + int nrvous = sptr->nrvous; tagint *atomIDs = sptr->atomIDs; - for (int i = 0; i < ncount; i++) + for (int i = 0; i < nrvous; i++) atom->map_one(atomIDs[i],i); // proclist = owner of atomID in caller decomposition @@ -1172,7 +1298,6 @@ int Special::rendezvous_1234(int n, char *inbuf, } outbuf = inbuf; - // NOTE: set out = in flag // re-create atom map @@ -1180,148 +1305,12 @@ int Special::rendezvous_1234(int n, char *inbuf, atom->nghost = 0; atom->map_set(); + // flag = 1: outbuf = inbuf + flag = 1; return n; } -/* ---------------------------------------------------------------------- - process data for atoms assigned to me in rendezvous decomposition - inbuf = list of N PairRvous datums - create outbuf = list of Nout PairRvous datums -------------------------------------------------------------------------- */ - -int Special::rendezvous_trim(int n, char *inbuf, - int &flag, int *&proclist, char *&outbuf, - void *ptr) -{ - int i,j,m; - - /* - Special *sptr = (Special *) ptr; - Atom *atom = sptr->atom; - Memory *memory = sptr->memory; - - // clear atom map so it can be here as a hash table - // faster than an STL map for large atom counts - - atom->map_clear(); - - // initialize hash - // ncount = number of atoms assigned to me - // key = atom ID - // value = index into Ncount-length data structure - - PairRvous *in = (PairRvous *) inbuf; - //std::map hash; - tagint id; - - int ncount = 0; - for (i = 0; i < n; i++) - if (in[i].me >= 0) - //hash[in[i].atomID] = ncount++; - atom->map_one(in[i].atomID,ncount++); - - // procowner = caller proc that owns each atom - // atomID = ID of each rendezvous atom I own - // npartner = # of 1-3 partners for each atom I own - - int *procowner,*npartner; - tagint *atomID; - memory->create(procowner,ncount,"special:procowner"); - memory->create(atomID,ncount,"special:atomID"); - memory->create(npartner,ncount,"special:npartner"); - for (m = 0; m < ncount; m++) npartner[m] = 0; - - for (i = 0; i < n; i++) { - //m = hash.find(in[i].atomID)->second; - m = atom->map(in[i].atomID); - if (in[i].me >= 0) { - procowner[m] = in[i].me; - atomID[m] = in[i].atomID; - } else if (in[i].me == -1) npartner[m]++; - } - - int max = 0; - for (m = 0; m < ncount; m++) max = MAX(max,npartner[m]); - - // partner = list of 1-3 or 1-4 partners for each atom I own - - int **partner; - memory->create(partner,ncount,max,"special:partner"); - for (m = 0; m < ncount; m++) npartner[m] = 0; - - for (i = 0; i < n; i++) { - if (in[i].me >= 0 || in[i].me == -2) continue; - //m = hash.find(in[i].atomID)->second; - m = atom->map(in[i].atomID); - partner[m][npartner[m]++] = in[i].partnerID; - } - - // flag = 1 if partner is in an actual angle or in a dihedral - - int **flag; - memory->create(flag,ncount,max,"special:flag"); - - for (i = 0; i < ncount; i++) - for (j = 0; j < npartner[i]; j++) - flag[i][j] = 0; - - tagint actual; - for (i = 0; i < n; i++) { - if (in[i].me != -2) continue; - actual = in[i].partnerID; - //m = hash.find(in[i].atomID)->second; - m = atom->map(in[i].atomID); - for (j = 0; j < npartner[m]; j++) - if (partner[m][j] == actual) { - flag[m][j] = 1; - break; - } - } - - // pass list of PairRvous datums back to comm->rendezvous - - int nout = 0; - for (m = 0; m < ncount; m++) nout += npartner[m]; - - memory->create(proclist,nout,"special:proclist"); - PairRvous *out = (PairRvous *) - memory->smalloc((bigint) nout*sizeof(PairRvous),"special:out"); - - nout = 0; - for (m = 0; m < ncount; m++) - for (j = 0; j < npartner[m]; j++) { - if (flag[m][j] == 0) continue; - proclist[nout] = procowner[m]; - out[nout].atomID = atomID[m]; - out[nout].partnerID = partner[m][j]; - nout++; - } - - outbuf = (char *) out; - - // clean up - // Comm::rendezvous will delete proclist and out (outbuf) - - memory->destroy(procowner); - memory->destroy(atomID); - memory->destroy(npartner); - memory->destroy(partner); - memory->destroy(flag); - - // re-create atom map - - atom->map_init(0); - atom->nghost = 0; - atom->map_set(); - - */ - - //return nout; - flag = 2; - return 0; -} - /* ---------------------------------------------------------------------- allow fixes to alter special list currently, only fix drude does this diff --git a/src/special.h b/src/special.h index 772ba613ac..d02a8522f6 100644 --- a/src/special.h +++ b/src/special.h @@ -31,7 +31,7 @@ class Special : protected Pointers { // data used by rendezvous callback methods - int ncount; + int nrvous; tagint *atomIDs; int *procowner; @@ -44,6 +44,8 @@ class Special : protected Pointers { tagint atomID,partnerID; }; + // private methods + void atom_owners(); void onetwo_build_newton(); void onetwo_build_newton_off(); @@ -60,8 +62,7 @@ class Special : protected Pointers { // callback functions for rendezvous communication static int rendezvous_ids(int, char *, int &, int *&, char *&, void *); - static int rendezvous_1234(int, char *, int &, int *&, char *&, void *); - static int rendezvous_trim(int, char *, int &, int *&, char *&, void *); + static int rendezvous_pairs(int, char *, int &, int *&, char *&, void *); }; }