}

}

/* ----------------------------------------------------------------------

   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

   atoms must be remapped to be inside simulation box before this is called

   for triclinic:

     atoms must be in lamda coords (0-1) before this is called

------------------------------------------------------------------------- */

void Comm::irregular()

{

  // clear global->local map since atoms move to new procs

  // zero out ghosts so map_set() at end will operate only on local atoms

  // exchange() doesn't need to zero ghosts b/c borders()

  //   is called right after and it zeroes ghosts and calls map_set()

  if (map_style) atom->map_clear();

  atom->nghost = 0;

  // 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()

  // if irregular_lookup() returns me, due to round-off

  //   in triclinic x2lamda(), then keep atom and don't send

  // 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]) {

      proclist[nsendatom] = irregular_lookup(x[i]);

      if (proclist[nsendatom] != me) {

	if (nsend > maxsend) grow_send(nsend,1);

	sizes[nsendatom] = avec->pack_exchange(i,&buf_send[nsend]);

	nsend += sizes[nsendatom];

	nsendatom++;

	avec->copy(nlocal-1,i);

	nlocal--;

      } else i++;

    } 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<int>

      (procgrid[0] * (x[0]-boxlo[0]) / (boxhi[0]-boxlo[0]));

    loc[1] = static_cast<int>

      (procgrid[1] * (x[1]-boxlo[1]) / (boxhi[1]-boxlo[1]));

    loc[2] = static_cast<int>

      (procgrid[2] * (x[2]-boxlo[2]) / (boxhi[2]-boxlo[2]));

  } else {

    loc[0] = static_cast<int> (procgrid[0] * x[0]);

    loc[1] = static_cast<int> (procgrid[1] * x[1]);

    loc[2] = static_cast<int> (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

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 ghost_velocity;               // 1 if ghost atoms have velocity, 0 if not

  int maxforward_fix;               // comm sizes called from Fix,Pair

  int maxreverse_fix;

  int maxforward_pair;

  int maxreverse_pair;

  double cutghost[3];               // cutoffs used for acquiring ghost atoms

  double cutghostuser;              // user-specified ghost cutoff

  int ***grid2proc;                 // which proc owns i,j,k loc in 3d grid

  Comm(class LAMMPS *);

  ~Comm();

  void forward_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 style;                        // single vs multi-type comm

  int nswap;                        // # of swaps to perform

  int need[3];                      // procs I need atoms from in each dim

  int triclinic;                    // 0 if domain is orthog, 1 if triclinic

  int maxswap;                      // max # of swaps memory is allocated for

  int size_forward;                 // # of per-atom datums in forward 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

  double cutghostuser;              // user-specified ghost cutoff

  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 bordergroup;                  // only communicate this group in borders

  int *firstrecv;                   // where to put 1st recv atom in each swap

  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,

  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 *);

};

}

#include "string.h"

#include "displace_atoms.h"

#include "atom.h"

#include "modify.h"

#include "domain.h"

#include "lattice.h"

#include "comm.h"

#include "irregular.h"

#include "group.h"

#include "random_park.h"

#include "error.h"

  if (domain->box_exist == 0) 

    error->all("Displace_atoms command before simulation box is defined");

  if (narg < 2) error->all("Illegal displace_atoms command");

  // init entire system since comm->irregular is done

  // comm::init needs neighbor::init needs pair::init needs kspace::init, etc

  if (comm->me == 0 && screen)

    fprintf(screen,"System init for displace_atoms ...\n");

  lmp->init();

  if (modify->nfix_restart_peratom) 

    error->all("Cannot displace_atoms after "

	       "reading restart file with per-atom info");

  if (comm->me == 0 && screen) fprintf(screen,"Displacing atoms ...\n");

  // move atoms randomly

  // makes atom result independent of what proc owns it via random->reset()

  if (style == RANDOM) {

    RanPark *random = new RanPark(lmp,1);

  // move atoms back inside simulation box and to new processors

  // use remap() instead of pbc() in case atoms moved a long distance

  // use irregular() instead of exchange() in case atoms moved a long distance

  // use irregular() in case atoms moved a long distance

  double **x = atom->x;

  int *image = atom->image;

  if (domain->triclinic) domain->x2lamda(atom->nlocal);

  domain->reset_box();

  comm->setup();

  comm->irregular();

  Irregular *irregular = new Irregular(lmp);

  irregular->migrate_atoms();

  delete irregular;

  if (domain->triclinic) domain->lamda2x(atom->nlocal);

  // check if any atoms were lost

#include "string.h"

#include "displace_box.h"

#include "atom.h"

#include "modify.h"

#include "domain.h"

#include "lattice.h"

#include "comm.h"

#include "irregular.h"

#include "group.h"

#include "error.h"

  if (domain->box_exist == 0) 

    error->all("Displace_box command before simulation box is defined");

  if (narg < 2) error->all("Illegal displace_box command");

  // init entire system since comm->irregular is done

  // comm::init needs neighbor::init needs pair::init needs kspace::init, etc

  if (comm->me == 0 && screen)

    fprintf(screen,"System init for displace_box ...\n");

  lmp->init();

  if (modify->nfix_restart_peratom) 

    error->all("Cannot displace_box after "

	       "reading restart file with per-atom info");

  if (comm->me == 0 && screen) fprintf(screen,"Displacing box ...\n");

  }

  // move atoms back inside simulation box and to new processors

  // use remap() instead of pbc() in case box

  //   moved a long distance relative to atoms

  // use irregular() instead of exchange() in case box

  //   moved a long distance relative to atoms

  // use remap() instead of pbc()

  //   in case box moved a long distance relative to atoms

  // use irregular() in case box moved a long distance relative to atoms

  double **x = atom->x;

  int *image = atom->image;

  if (domain->triclinic) domain->x2lamda(atom->nlocal);

  domain->reset_box();

  comm->setup();

  comm->irregular();

  Irregular *irregular = new Irregular(lmp);

  irregular->migrate_atoms();

  delete irregular;

  if (domain->triclinic) domain->lamda2x(atom->nlocal);

  // clean up

#include "atom.h"

#include "update.h"

#include "comm.h"

#include "irregular.h"

#include "domain.h"

#include "lattice.h"

#include "force.h"

  rfix = NULL;

  flip = 0;

  if (force_reneighbor) irregular = new Irregular(lmp);

  else irregular = NULL;

  TWOPI = 8.0*atan(1.0);

}

  delete [] set;

  delete [] rfix;

  delete irregular;

  // reset domain's h_rate = 0.0, since this fix may have made it non-zero

  double *h_rate = domain->h_rate;

int FixDeform::setmask()

{

  int mask = 0;

  mask |= PRE_EXCHANGE;

  if (force_reneighbor) mask |= PRE_EXCHANGE;

  mask |= END_OF_STEP;

  return mask;

}

  for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);

  domain->x2lamda(atom->nlocal);

  comm->irregular();

  irregular->migrate_atoms();

  domain->lamda2x(atom->nlocal);

  flip = 0;