/* ----------------------------------------------------------------------
    This is the

    ██╗     ██╗ ██████╗  ██████╗  ██████╗ ██╗  ██╗████████╗███████╗
    ██║     ██║██╔════╝ ██╔════╝ ██╔════╝ ██║  ██║╚══██╔══╝██╔════╝
    ██║     ██║██║  ███╗██║  ███╗██║  ███╗███████║   ██║   ███████╗
    ██║     ██║██║   ██║██║   ██║██║   ██║██╔══██║   ██║   ╚════██║
    ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║  ██║   ██║   ███████║
    ╚══════╝╚═╝ ╚═════╝  ╚═════╝  ╚═════╝ ╚═╝  ╚═╝   ╚═╝   ╚══════╝®

    DEM simulation engine, released by
    DCS Computing Gmbh, Linz, Austria
    http://www.dcs-computing.com, office@dcs-computing.com

    LIGGGHTS® is part of CFDEM®project:
    http://www.liggghts.com | http://www.cfdem.com

    Core developer and main author:
    Christoph Kloss, christoph.kloss@dcs-computing.com

    LIGGGHTS® is open-source, distributed under the terms of the GNU Public
    License, version 2 or later. It is distributed in the hope that it will
    be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
    of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
    received a copy of the GNU General Public License along with LIGGGHTS®.
    If not, see http://www.gnu.org/licenses . See also top-level README
    and LICENSE files.

    LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
    the producer of the LIGGGHTS® software and the CFDEM®coupling software
    See http://www.cfdem.com/terms-trademark-policy for details.

-------------------------------------------------------------------------
    Contributing author and copyright for this file:

    Christoph Kloss (DCS Computing GmbH, Linz)
    Christoph Kloss (JKU Linz)
    Arno Mayrhofer (DCS Computing GmbH, Linz)

    This file is from LAMMPS, but has been modified. Copyright for
    modification:

    Copyright 2012-     DCS Computing GmbH, Linz
    Copyright 2009-2012 JKU Linz

    Copyright of original file:
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov

    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing author (triclinic) : Pieter in 't Veld (SNL)
------------------------------------------------------------------------- */

#include <mpi.h>
#include <cmath>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include "comm_brick.h"
#include "universe.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "domain_wedge.h"
#include "math_extra_liggghts.h"
#include "neighbor.h"
#include "group.h"
#include "modify.h"
#include "fix.h"
#include "compute.h"
#include "output.h"
#include "dump.h"
#include "math_extra.h"
#include "error.h"
#include "memory.h"

using namespace LAMMPS_NS;

#define BUFFACTOR 1.5
#define BUFMIN 1000
#define BUFEXTRA 1000
#define BIG 1.0e20

enum{SINGLE,MULTI};               // same as in Comm

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

CommBrick::CommBrick(LAMMPS *lmp) :
  Comm(lmp),
  sendnum(NULL), recvnum(NULL), sendproc(NULL), recvproc(NULL),
  size_forward_recv(NULL),
  size_reverse_send(NULL), size_reverse_recv(NULL),
  slablo(NULL), slabhi(NULL), multilo(NULL), multihi(NULL),
  cutghostmulti(NULL), pbc_flag(NULL), pbc(NULL), firstrecv(NULL),
  sendlist(NULL), maxsendlist(NULL)
{
  style = 0;
  layout = LAYOUT_UNIFORM;
  pbc_flag = NULL;
  init_buffers();
}

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

CommBrick::~CommBrick()
{
  free_swap();
  if (style == MULTI) {
    free_multi();
  }

  if (sendlist) for (int i = 0; i < maxswap; i++) memory->destroy(sendlist[i]);
  memory->sfree(sendlist);
  memory->destroy(maxsendlist);
}

/* ---------------------------------------------------------------------- */
//IMPORTANT: we *MUST* pass "*oldcomm" to the Comm initializer here, as
//           the code below *requires* that the (implicit) copy constructor
//           for Comm is run and thus creating a shallow copy of "oldcomm".
//           The call to Comm::copy_arrays() then converts the shallow copy
//           into a deep copy of the class with the new layout.

CommBrick::CommBrick(LAMMPS *lmp, Comm *oldcomm) : Comm(*oldcomm)
{
  if (oldcomm->get_layout() == LAYOUT_TILED)
    error->all(FLERR,"Cannot change to comm_style brick from tiled layout");

  style = 0;
  layout = oldcomm->get_layout();
  Comm::copy_arrays(oldcomm);
  init_buffers();
}

CommBrick::CommBrick(Comm *main_comm, ParallelBase *pb) :
  Comm(*main_comm)
{
  style = 0;
  is_subcomm = true;
  pb_ = pb;
  main_comm->register_subcomm(this);
  Comm::copy_pointers(main_comm);
  init_buffers();
  subcomms_.clear();
}

/* ----------------------------------------------------------------------
   initialize comm buffers and other data structs local to CommBrick
------------------------------------------------------------------------- */

void CommBrick::init_buffers()
{
  multilo = multihi = NULL;
  cutghostmulti = NULL;

  // bufextra = max size of one exchanged atom
  //          = allowed overflow of sendbuf in exchange()
  // atomvec, fix reset these 2 maxexchange values if needed
  // only necessary if their size > BUFEXTRA

  maxexchange = maxexchange_atom + maxexchange_fix;
  bufextra = maxexchange + BUFEXTRA;

  maxsend = BUFMIN;
  memory->create(buf_send,maxsend+bufextra,"comm:buf_send");
  maxrecv = BUFMIN;
  memory->create(buf_recv,maxrecv,"comm:buf_recv");

  nswap = 0;
  maxswap = 6;
  allocate_swap(maxswap);

  sendlist = (int **) memory->smalloc(maxswap*sizeof(int *),"comm:sendlist");
  memory->create(maxsendlist,maxswap,"comm:maxsendlist");
  for (int i = 0; i < maxswap; i++) {
    maxsendlist[i] = BUFMIN;
    memory->create(sendlist[i],BUFMIN,"comm:sendlist[i]");
  }
}

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

void CommBrick::init()
{
  Comm::init();

  // memory for multi-style communication

  if (mode == MULTI && multilo == NULL) {
    allocate_multi(maxswap);
    memory->create(cutghostmulti,atom->ntypes+1,3,"comm:cutghostmulti");
  }
  if (mode == SINGLE && multilo) {
    free_multi();
    memory->destroy(cutghostmulti);
  }
}

/* ----------------------------------------------------------------------
   setup spatial-decomposition communication patterns
   function of neighbor cutoff(s) & cutghostuser & current box size
   single mode sets slab boundaries (slablo,slabhi) based on max cutoff
   multi mode sets type-dependent slab boundaries (multilo,multihi)
------------------------------------------------------------------------- */

void CommBrick::setup()
{
  // cutghost[] = max distance at which ghost atoms need to be acquired
  // for orthogonal:
  //   cutghost is in box coords = neigh->cutghost in all 3 dims
  // for triclinic:
  //   neigh->cutghost = distance between tilted planes in box coords
  //   cutghost is in lamda coords = distance between those planes
  // for multi:
  //   cutghostmulti = same as cutghost, only for each atom type

  double *prd,*sublo,*subhi;

  double cut;
  if (pb_)
    cut = pb_->get_cutoff();
  else
  {
    cut = MAX(neighbor->cutneighmax,cutghostuser);

    int nfix = modify->nfix;
    Fix **fix = modify->fix;

    for(int ifix = 0; ifix < nfix; ifix++)
    {
        
        if(!fix[ifix]->use_rad_for_cut_neigh_and_ghost())
          continue;

        double cut_fix = fix[ifix]->extend_cut_ghost();
        cut = MAX(cut,cut_fix + neighbor->skin);
        
    }
  }

  if (domain->triclinic == 0) {
    prd = domain->prd;
    sublo = domain->sublo;
    subhi = domain->subhi;
    
    cutghost[0] = cutghost[1] = cutghost[2] = cut;

    if (mode == MULTI) {
      double *cuttype = neighbor->cuttype;
      for (int i = 1; i <= atom->ntypes; i++) {
        cut = 0.0;
        if (cutusermulti) cut = cutusermulti[i];
        cutghostmulti[i][0] = MAX(cut,cuttype[i]);
        cutghostmulti[i][1] = MAX(cut,cuttype[i]);
        cutghostmulti[i][2] = MAX(cut,cuttype[i]);
      }
    }

  } else {
    prd = domain->prd_lamda;
    sublo = domain->sublo_lamda;
    subhi = domain->subhi_lamda;
    double *h_inv = domain->h_inv;
    double length0,length1,length2;
    length0 = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]);
    cutghost[0] = cut * length0;
    length1 = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]);
    cutghost[1] = cut * length1;
    length2 = h_inv[2];
    cutghost[2] = cut * length2;

    if (mode == MULTI) {
      double *cuttype = neighbor->cuttype;
      for (int i = 1; i <= atom->ntypes; i++) {
        cut = 0.0;
        if (cutusermulti) cut = cutusermulti[i];
        cutghostmulti[i][0] = length0 * MAX(cut,cuttype[i]);
        cutghostmulti[i][1] = length1 * MAX(cut,cuttype[i]);
        cutghostmulti[i][2] = length2 * MAX(cut,cuttype[i]);
      }
    }
  }

  // recvneed[idim][0/1] = # of procs away I recv atoms from, within cutghost
  //   0 = from left, 1 = from right
  //   do not cross non-periodic boundaries, need[2] = 0 for 2d
  // sendneed[idim][0/1] = # of procs away I send atoms to
  //   0 = to left, 1 = to right
  //   set equal to recvneed[idim][1/0] of neighbor proc
  // maxneed[idim] = max procs away any proc recvs atoms in either direction
  // layout = UNIFORM = uniform sized sub-domains:
  //   maxneed is directly computable from sub-domain size
  //     limit to procgrid-1 for non-PBC
  //   recvneed = maxneed except for procs near non-PBC
  //   sendneed = recvneed of neighbor on each side
  // layout = NONUNIFORM = non-uniform sized sub-domains:
  //   compute recvneed via updown() which accounts for non-PBC
  //   sendneed = recvneed of neighbor on each side
  //   maxneed via Allreduce() of recvneed

  int *periodicity = domain->periodicity;
  int left,right;

  if (layout == LAYOUT_UNIFORM) {
    maxneed[0] = static_cast<int> (cutghost[0] * procgrid[0] / prd[0]) + 1;
    maxneed[1] = static_cast<int> (cutghost[1] * procgrid[1] / prd[1]) + 1;
    maxneed[2] = static_cast<int> (cutghost[2] * procgrid[2] / prd[2]) + 1;
    if (domain->dimension == 2) maxneed[2] = 0;
    if (!periodicity[0]) maxneed[0] = MIN(maxneed[0],procgrid[0]-1);
    if (!periodicity[1]) maxneed[1] = MIN(maxneed[1],procgrid[1]-1);
    if (!periodicity[2]) maxneed[2] = MIN(maxneed[2],procgrid[2]-1);

    if (!periodicity[0]) {
      recvneed[0][0] = MIN(maxneed[0],myloc[0]);
      recvneed[0][1] = MIN(maxneed[0],procgrid[0]-myloc[0]-1);
      left = myloc[0] - 1;
      if (left < 0) left = procgrid[0] - 1;
      sendneed[0][0] = MIN(maxneed[0],procgrid[0]-left-1);
      right = myloc[0] + 1;
      if (right == procgrid[0]) right = 0;
      sendneed[0][1] = MIN(maxneed[0],right);
    } else recvneed[0][0] = recvneed[0][1] =
             sendneed[0][0] = sendneed[0][1] = maxneed[0];

    if (!periodicity[1]) {
      recvneed[1][0] = MIN(maxneed[1],myloc[1]);
      recvneed[1][1] = MIN(maxneed[1],procgrid[1]-myloc[1]-1);
      left = myloc[1] - 1;
      if (left < 0) left = procgrid[1] - 1;
      sendneed[1][0] = MIN(maxneed[1],procgrid[1]-left-1);
      right = myloc[1] + 1;
      if (right == procgrid[1]) right = 0;
      sendneed[1][1] = MIN(maxneed[1],right);
    } else recvneed[1][0] = recvneed[1][1] =
             sendneed[1][0] = sendneed[1][1] = maxneed[1];

    if (!periodicity[2]) {
      recvneed[2][0] = MIN(maxneed[2],myloc[2]);
      recvneed[2][1] = MIN(maxneed[2],procgrid[2]-myloc[2]-1);
      left = myloc[2] - 1;
      if (left < 0) left = procgrid[2] - 1;
      sendneed[2][0] = MIN(maxneed[2],procgrid[2]-left-1);
      right = myloc[2] + 1;
      if (right == procgrid[2]) right = 0;
      sendneed[2][1] = MIN(maxneed[2],right);
    } else recvneed[2][0] = recvneed[2][1] =
             sendneed[2][0] = sendneed[2][1] = maxneed[2];

  } else {
    recvneed[0][0] = updown(0,0,myloc[0],prd[0],periodicity[0],xsplit);
    recvneed[0][1] = updown(0,1,myloc[0],prd[0],periodicity[0],xsplit);
    left = myloc[0] - 1;
    if (left < 0) left = procgrid[0] - 1;
    sendneed[0][0] = updown(0,1,left,prd[0],periodicity[0],xsplit);
    right = myloc[0] + 1;
    if (right == procgrid[0]) right = 0;
    sendneed[0][1] = updown(0,0,right,prd[0],periodicity[0],xsplit);

    recvneed[1][0] = updown(1,0,myloc[1],prd[1],periodicity[1],ysplit);
    recvneed[1][1] = updown(1,1,myloc[1],prd[1],periodicity[1],ysplit);
    left = myloc[1] - 1;
    if (left < 0) left = procgrid[1] - 1;
    sendneed[1][0] = updown(1,1,left,prd[1],periodicity[1],ysplit);
    right = myloc[1] + 1;
    if (right == procgrid[1]) right = 0;
    sendneed[1][1] = updown(1,0,right,prd[1],periodicity[1],ysplit);

    if (domain->dimension == 3) {
      recvneed[2][0] = updown(2,0,myloc[2],prd[2],periodicity[2],zsplit);
      recvneed[2][1] = updown(2,1,myloc[2],prd[2],periodicity[2],zsplit);
      left = myloc[2] - 1;
      if (left < 0) left = procgrid[2] - 1;
      sendneed[2][0] = updown(2,1,left,prd[2],periodicity[2],zsplit);
      right = myloc[2] + 1;
      if (right == procgrid[2]) right = 0;
      sendneed[2][1] = updown(2,0,right,prd[2],periodicity[2],zsplit);
    } else recvneed[2][0] = recvneed[2][1] =
             sendneed[2][0] = sendneed[2][1] = 0;

    int all[6];
    MPI_Allreduce(&recvneed[0][0],all,6,MPI_INT,MPI_MAX,world);
    maxneed[0] = MAX(all[0],all[1]);
    maxneed[1] = MAX(all[2],all[3]);
    maxneed[2] = MAX(all[4],all[5]);
  }

  // allocate comm memory

  nswap = 2 * (maxneed[0]+maxneed[1]+maxneed[2]);
  if (nswap > maxswap) grow_swap(nswap);

  dw_ = dynamic_cast<DomainWedge*>(domain);
  if(dw_)
  {
      ia = dw_->index_axis();
      iphi = dw_->index_phi();
      dw_->n1(nleft);
      dw_->n2(nright);
      dw_->center(c);
      double cutghmax = MathExtraLiggghts::max(cutghost[0],cutghost[1],cutghost[2]);
      pleft[0]  = c[0] - nleft[0]  * (use_gran_opt() ? (cutghmax / 2. + neighbor->skin/2.) : cutghmax);
      pleft[1]  = c[1] - nleft[1]  * (use_gran_opt() ? (cutghmax / 2. + neighbor->skin/2.) : cutghmax);
      pright[0] = c[0] - nright[0] * (use_gran_opt() ? (cutghmax / 2. + neighbor->skin/2.) : cutghmax);
      pright[1] = c[1] - nright[1] * (use_gran_opt() ? (cutghmax / 2. + neighbor->skin/2.) : cutghmax);
      
  }

  // setup parameters for each exchange:
  // sendproc = proc to send to at each swap
  // recvproc = proc to recv from at each swap
  // for mode SINGLE:
  //   slablo/slabhi = boundaries for slab of atoms to send at each swap
  //   use -BIG/midpt/BIG to insure all atoms included even if round-off occurs
  //   if round-off, atoms recvd across PBC can be < or > than subbox boundary
  //   note that borders() only loops over subset of atoms during each swap
  //   treat all as PBC here, non-PBC is handled in borders() via r/s need[][]
  // for mode MULTI:
  //   multilo/multihi is same, with slablo/slabhi for each atom type
  // pbc_flag: 0 = nothing across a boundary, 1 = something across a boundary
  // pbc = -1/0/1 for PBC factor in each of 3/6 orthogonal/triclinic dirs
  // for triclinic, slablo/hi and pbc_border will be used in lamda (0-1) coords
  // 1st part of if statement is sending to the west/south/down
  // 2nd part of if statement is sending to the east/north/up

  int dim,ineed;

  int iswap = 0;
  for (dim = 0; dim < 3; dim++) {
    double slab_offset;
    
    if (pb_)
      slab_offset = pb_->get_slab_offset(dim);
    else
      slab_offset = use_gran_opt() ? (cutghost[dim] / 2. + neighbor->skin/2.) : cutghost[dim];
    for (ineed = 0; ineed < 2*maxneed[dim]; ineed++) {
      
      if(dw_ && dim != ia && dim != iphi)
      {
        
        nswap--;
        continue;
      }

      pbc_flag[iswap] = 0;
      pbc[iswap][0] = pbc[iswap][1] = pbc[iswap][2] =
        pbc[iswap][3] = pbc[iswap][4] = pbc[iswap][5] = 0;

      if (ineed % 2 == 0) {
        sendproc[iswap] = procneigh[dim][0];
        recvproc[iswap] = procneigh[dim][1];
        if (mode == SINGLE) {
          if (ineed < 2) slablo[iswap] = -BIG;
          else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
          slabhi[iswap] = sublo[dim] + slab_offset; 
        } else {
          for (int i = 1; i <= atom->ntypes; i++) {
            if (ineed < 2) multilo[iswap][i] = -BIG;
            else multilo[iswap][i] = 0.5 * (sublo[dim] + subhi[dim]);
            multihi[iswap][i] = sublo[dim] + cutghostmulti[i][dim];
          }
        }
        if (myloc[dim] == 0) {
          pbc_flag[iswap] = 1;
          pbc[iswap][dim] = 1;
          if (domain->triclinic) {
            if (dim == 1) pbc[iswap][5] = 1;
            else if (dim == 2) pbc[iswap][4] = pbc[iswap][3] = 1;
          }
        }

      } else {
        sendproc[iswap] = procneigh[dim][1];
        recvproc[iswap] = procneigh[dim][0];
        if (style == SINGLE) {
          slablo[iswap] = subhi[dim] - slab_offset; 
          if (ineed < 2) slabhi[iswap] = BIG;
          else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
        } else {
          for (int i = 1; i <= atom->ntypes; i++) {
            multilo[iswap][i] = subhi[dim] - cutghostmulti[i][dim];
            if (ineed < 2) multihi[iswap][i] = BIG;
            else multihi[iswap][i] = 0.5 * (sublo[dim] + subhi[dim]);
          }
        }
        if (myloc[dim] == procgrid[dim]-1) {
          pbc_flag[iswap] = 1;
          pbc[iswap][dim] = -1;
          if (domain->triclinic) {
            if (dim == 1) pbc[iswap][5] = -1;
            else if (dim == 2) pbc[iswap][4] = pbc[iswap][3] = -1;
          }
        }
      }

      iswap++;
    }
  }
}

/* ----------------------------------------------------------------------
   walk up/down the extent of nearby processors in dim and dir
   loc = myloc of proc to start at
   dir = 0/1 = walk to left/right
   do not cross non-periodic boundaries
   is not called for z dim in 2d
   return how many procs away are needed to encompass cutghost away from loc
------------------------------------------------------------------------- */

int CommBrick::updown(int dim, int dir, int loc,
                 double prd, int periodicity, double *split)
{
  int index,count;
  double frac,delta;

  if (dir == 0) {
    frac = cutghost[dim]/prd;
    index = loc - 1;
    delta = 0.0;
    count = 0;
    while (delta < frac) {
      if (index < 0) {
        if (!periodicity) break;
        index = procgrid[dim] - 1;
      }
      count++;
      delta += split[index+1] - split[index];
      index--;
    }

  } else {
    frac = cutghost[dim]/prd;
    index = loc + 1;
    delta = 0.0;
    count = 0;
    while (delta < frac) {
      if (index >= procgrid[dim]) {
        if (!periodicity) break;
        index = 0;
      }
      count++;
      delta += split[index+1] - split[index];
      index++;
    }
  }

  return count;
}

/* ----------------------------------------------------------------------
   forward communication of atom coords every timestep
   other per-atom attributes may also be sent via pack/unpack routines
------------------------------------------------------------------------- */

void CommBrick::forward_comm(int dummy)
{
  MPI_Request request;
  AtomVec *avec = atom->avec;
  double **x = atom->x;
  double *buf;

  double nrecv_scale[2] = {0., 0.};
  if (pb_) {
    pb_->get_nrecv_scale(OPERATION_COMM_FORWARD, nrecv_scale);

    // check if there is anything to transfer at all
    if (nrecv_scale[0] == 0)
      return;
  }

  // exchange data with another proc
  // if other proc is self, just copy
  // if comm_x_only set, exchange or copy directly to x, don't unpack

  for (int iswap = 0; iswap < nswap; iswap++) {
    if (sendproc[iswap] != me) {
      if (pb_) {
        const int size_recv = (size_forward_recv[iswap]*nrecv_scale[0])/nrecv_scale[1];
        if (size_recv)
          MPI_Irecv(buf_recv,size_recv,MPI_DOUBLE,
                    recvproc[iswap],0,world,&request);
        const int n = pb_->pack_comm(OPERATION_COMM_FORWARD, sendnum[iswap], sendlist[iswap],
                                     buf_send, pbc_flag[iswap], pbc[iswap]);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
        if (size_recv) MPI_Wait(&request,MPI_STATUS_IGNORE);
        pb_->unpack_comm(OPERATION_COMM_FORWARD, recvnum[iswap], firstrecv[iswap], buf_recv);
      } else if (comm_x_only) {
        if (size_forward_recv[iswap]) {
          buf = x[firstrecv[iswap]];
          MPI_Irecv(buf,size_forward_recv[iswap],MPI_DOUBLE,
                    recvproc[iswap],0,world,&request);
        }
        const int n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
                                      buf_send,pbc_flag[iswap],pbc[iswap]);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
        if (size_forward_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      } else if (ghost_velocity) {
        if (size_forward_recv[iswap])
          MPI_Irecv(buf_recv,size_forward_recv[iswap],MPI_DOUBLE,
                    recvproc[iswap],0,world,&request);
        const int n = avec->pack_comm_vel(sendnum[iswap],sendlist[iswap],
                                          buf_send,pbc_flag[iswap],pbc[iswap]);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
        if (size_forward_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
        avec->unpack_comm_vel(recvnum[iswap],firstrecv[iswap],buf_recv);
      } else {
        if (size_forward_recv[iswap])
          MPI_Irecv(buf_recv,size_forward_recv[iswap],MPI_DOUBLE,
                    recvproc[iswap],0,world,&request);
        const int n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
                                      buf_send,pbc_flag[iswap],pbc[iswap]);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
        if (size_forward_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
        avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_recv);
      }

    } else {
      if (pb_) {
        pb_->pack_comm(OPERATION_COMM_FORWARD, sendnum[iswap],sendlist[iswap],
                       buf_send,pbc_flag[iswap],pbc[iswap]);
        pb_->unpack_comm(OPERATION_COMM_FORWARD, recvnum[iswap], firstrecv[iswap], buf_send);
      } else if (comm_x_only) {
        if (sendnum[iswap])
          avec->pack_comm(sendnum[iswap],sendlist[iswap],
                          x[firstrecv[iswap]],pbc_flag[iswap],pbc[iswap]);
      } else if (ghost_velocity) {
        avec->pack_comm_vel(sendnum[iswap],sendlist[iswap],
                            buf_send,pbc_flag[iswap],pbc[iswap]);
        avec->unpack_comm_vel(recvnum[iswap],firstrecv[iswap],buf_send);
      } else {
        avec->pack_comm(sendnum[iswap],sendlist[iswap],
                        buf_send,pbc_flag[iswap],pbc[iswap]);
        avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send);
      }
    }
  }
}

/* ----------------------------------------------------------------------
   reverse communication of forces on atoms every timestep
   other per-atom attributes may also be sent via pack/unpack routines
------------------------------------------------------------------------- */

void CommBrick::reverse_comm()
{
  int n;
  MPI_Request request;
  AtomVec *avec = atom->avec;
  double **f = atom->f;
  double *buf;

  double nrecv_scale[2] = {0., 0.};
  if (pb_) {
    pb_->get_nrecv_scale(OPERATION_COMM_REVERSE, nrecv_scale);

    // check if there is anything to transfer at all
    if (nrecv_scale[0] == 0)
      return;
  }

  // exchange data with another proc
  // if other proc is self, just copy
  // if comm_f_only set, exchange or copy directly from f, don't pack

  for (int iswap = nswap-1; iswap >= 0; iswap--) {
    if (sendproc[iswap] != me) {
      if (pb_) {
        const int size_recv = (size_reverse_recv[iswap]*nrecv_scale[0])/nrecv_scale[1];
        if (size_recv)
          MPI_Irecv(buf_recv,size_recv,MPI_DOUBLE,
                    sendproc[iswap],0,world,&request);
        n = pb_->pack_reverse(OPERATION_COMM_REVERSE, recvnum[iswap],firstrecv[iswap],buf_send);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap],0,world);
        if (size_recv) MPI_Wait(&request,MPI_STATUS_IGNORE);
        pb_->unpack_reverse(OPERATION_COMM_REVERSE, sendnum[iswap],sendlist[iswap],buf_recv);
      } else if (comm_f_only) {
        if (size_reverse_recv[iswap])
          MPI_Irecv(buf_recv,size_reverse_recv[iswap],MPI_DOUBLE,
                    sendproc[iswap],0,world,&request);
        if (size_reverse_send[iswap]) {
          buf = f[firstrecv[iswap]];
          MPI_Send(buf,size_reverse_send[iswap],MPI_DOUBLE,
                   recvproc[iswap],0,world);
        }
        if (size_reverse_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      } else {
        if (size_reverse_recv[iswap])
          MPI_Irecv(buf_recv,size_reverse_recv[iswap],MPI_DOUBLE,
                    sendproc[iswap],0,world,&request);
        n = avec->pack_reverse(recvnum[iswap],firstrecv[iswap],buf_send);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap],0,world);
        if (size_reverse_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      }
      if (!pb_)
        avec->unpack_reverse(sendnum[iswap],sendlist[iswap],buf_recv);

    } else {
      if (pb_) {
        pb_->pack_reverse(OPERATION_COMM_REVERSE, recvnum[iswap],firstrecv[iswap],buf_send);
        pb_->unpack_reverse(OPERATION_COMM_REVERSE, sendnum[iswap],sendlist[iswap],buf_send);
      } else if (comm_f_only) {
        if (sendnum[iswap])
          avec->unpack_reverse(sendnum[iswap],sendlist[iswap],
                               f[firstrecv[iswap]]);
      } else {
        avec->pack_reverse(recvnum[iswap],firstrecv[iswap],buf_send);
        avec->unpack_reverse(sendnum[iswap],sendlist[iswap],buf_send);
      }
    }
  }
}

/* ----------------------------------------------------------------------
   exchange: move atoms to correct processors
   atoms exchanged with all 6 stencil neighbors
   send out atoms that have left my box, receive ones entering my box
   atoms will be lost if not inside a stencil proc's box
     can happen if atom moves outside of non-periodic bounary
     or if atom moves more than one proc away
   this routine called before every reneighboring
   for triclinic, atoms must be in lamda coords (0-1) before exchange is called
------------------------------------------------------------------------- */

void CommBrick::exchange()
{
  int i,m,nsend,nrecv,nrecv1,nrecv2,nlocal;
  double lo = 0., hi = 0., value;
  double **x;
  double *sublo,*subhi;
  MPI_Request request;
  AtomVec *avec = atom->avec;

  if (pb_)
    pb_->clear_exchange();
    
  else
  {
    // clear global->local map for owned and ghost atoms
    // b/c atoms migrate to new procs in exchange() and
    // new ghosts are created in borders()
    // map_set() is done at end of borders()
    // clear ghost count and any ghost bonus data internal to AtomVec

    if (map_style) atom->map_clear();
    atom->nghost = 0;
    atom->avec->clear_bonus();

    // insure send buf is large enough for single atom
    // bufextra = max size of one atom = allowed overflow of sendbuf
    // fixes can change per-atom size requirement on-the-fly

    int bufextra_old = bufextra;
    maxexchange = maxexchange_atom + maxexchange_fix;
    bufextra = maxexchange + BUFEXTRA;
    if (bufextra > bufextra_old)
      memory->grow(buf_send,maxsend+bufextra,"comm:buf_send");
  }

  // subbox bounds for orthogonal or triclinic

  if (domain->triclinic == 0) {
    sublo = domain->sublo;
    subhi = domain->subhi;
  } else {
    sublo = domain->sublo_lamda;
    subhi = domain->subhi_lamda;
  }

  //Record exchange events in advance if necessary
  exchangeEventsRecorder();

  // loop over dimensions

  const int dimension = domain->dimension;

  for (int dim = 0; dim < dimension; dim++) {

    // fill buffer with atoms leaving my box, using < and >=
    // when atom is deleted, fill it in with last atom

    if (pb_)
      nsend = pb_->pack_exchange(dim);
    else
    {
      x = atom->x;
      lo = sublo[dim];
      hi = subhi[dim];
      nlocal = atom->nlocal;
      i = nsend = 0;

      while (i < nlocal) {
        if (x[i][dim] < lo || x[i][dim] >= hi) {
          if (nsend > maxsend) grow_send(nsend,1);
          nsend += avec->pack_exchange(i,&buf_send[nsend]);
          avec->copy(nlocal-1,i,1);
          nlocal--;
        } else i++;
      }
      atom->nlocal = nlocal;
    }

    // send/recv atoms in both directions
    // send size of message first so receiver can realloc buf_recv if needed
    // if 1 proc in dimension, no send/recv
    //   set nrecv = 0 so buf_send atoms will be lost
    // if 2 procs in dimension, single send/recv
    // if more than 2 procs in dimension, send/recv to both neighbors

    if (procgrid[dim] == 1) nrecv = 0;
    else {
      MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][0],0,
                   &nrecv1,1,MPI_INT,procneigh[dim][1],0,world,
                   MPI_STATUS_IGNORE);
      nrecv = nrecv1;
      if (procgrid[dim] > 2) {
        MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][1],0,
                     &nrecv2,1,MPI_INT,procneigh[dim][0],0,world,
                     MPI_STATUS_IGNORE);
        nrecv += nrecv2;
      }
      if (nrecv > maxrecv) grow_recv(nrecv);

      MPI_Irecv(buf_recv,nrecv1,MPI_DOUBLE,procneigh[dim][1],0,
                world,&request);
      MPI_Send(buf_send,nsend,MPI_DOUBLE,procneigh[dim][0],0,world);
      MPI_Wait(&request,MPI_STATUS_IGNORE);

      if (procgrid[dim] > 2) {
        MPI_Irecv(&buf_recv[nrecv1],nrecv2,MPI_DOUBLE,procneigh[dim][0],0,
                  world,&request);
        MPI_Send(buf_send,nsend,MPI_DOUBLE,procneigh[dim][1],0,world);
        MPI_Wait(&request,MPI_STATUS_IGNORE);
      }
    }

    if (pb_)
      pb_->unpack_exchange(nrecv, dim);
    else
    {
      // check incoming atoms to see if they are in my box
      // if so, add to my list
      // box check is only for this dimension,
      //   atom may be passed to another proc in later dims

      m = 0;
      while (m < nrecv) {
        value = buf_recv[m+dim+1];
        if (value >= lo && value < hi) m += avec->unpack_exchange(&buf_recv[m]);
        else m += static_cast<int> (buf_recv[m]);
      }
    }
  }

  if (pb_)
    pb_->post_exchange();
  else
  {
    
    bigint tmp = atom->nlocal;
    MPI_Allreduce(&tmp,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);

    if (atom->firstgroupname) atom->first_reorder();
  }
}

/* ----------------------------------------------------------------------
   borders: list nearby atoms to send to neighboring procs at every timestep
   one list is created for every swap that will be made
   as list is made, actually do swaps
   this does equivalent of a forward_comm(), so don't need to explicitly
     call forward_comm() on reneighboring timestep
   this routine is called before every reneighboring
   for triclinic, atoms must be in lamda coords (0-1) before borders is called
------------------------------------------------------------------------- */

void CommBrick::borders()
{
  int i,n,iswap,dim,ineed,twoneed;
  int nsend,nrecv,sendflag,nfirst,nlast,ngroup;
  double lo = 0.0,hi = 0.0;
  double **x;
  double *buf, *mlo = NULL, *mhi = NULL;
  MPI_Request request;
  AtomVec *avec = atom->avec;

  // do swaps over all 3 dimensions

  nfirst = 0;
  iswap = 0;
  smax = rmax = 0;

  for (dim = 0; dim < 3; dim++) {
    nlast = 0;
    twoneed = 2*maxneed[dim];
    for (ineed = 0; ineed < twoneed; ineed++) {

      if(dw_ && dim != ia && dim != iphi)
      {
        
        continue;
      }

      // find atoms within slab boundaries lo/hi using <= and >=
      // check atoms between nfirst and nlast
      //   for first swaps in a dim, check owned and ghost
      //   for later swaps in a dim, only check newly arrived ghosts
      // store sent atom indices in sendlist for use in future timesteps

      x = atom->x;
      if (mode == SINGLE) {
        lo = slablo[iswap];
        hi = slabhi[iswap];
      } else {
        mlo = multilo[iswap];
        mhi = multihi[iswap];
      }
      if (ineed % 2 == 0) {
        nfirst = nlast;
        if (pb_)
          nlast = pb_->get_ntotal();
        else
          nlast = atom->nlocal + atom->nghost;
      }

      nsend = 0;

      // sendflag = 0 if I do not send on this swap
      // sendneed test indicates receiver no longer requires data
      // e.g. due to non-PBC or non-uniform sub-domains

      if (ineed/2 >= sendneed[dim][ineed % 2]) sendflag = 0;
      else sendflag = 1;

      // find send atoms according to SINGLE vs MULTI
      // all atoms eligible versus only atoms in bordergroup
      // can only limit loop to bordergroup for first sends (ineed < 2)
      // on these sends, break loop in two: owned (in group) and ghost

      if (sendflag) {
        if (!bordergroup || ineed >= 2) {
          if (bordergroup && pb_)
            error->all(FLERR, "Group command not implemented for ParallelBase");
          if (mode == SINGLE) {
            
            if(dw_ && dim == iphi)
            {
                for (i = nfirst; i < nlast; i++)
                  if (decide_wedge(i,dim,lo,hi,ineed)) {
                    if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                    sendlist[iswap][nsend++] = i;
                  }
            }
            else 
            {
                for (i = nfirst; i < nlast; i++)
                  if (decide(i,dim,lo,hi,ineed)) {
                    if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                    sendlist[iswap][nsend++] = i;
                  }
            }
          } else {
            if (pb_)
              error->all(FLERR, "Group command not implemented for ParallelBase");
            for (i = nfirst; i < nlast; i++) {
              const int itype = get_type(i);
              if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
              }
            }
          }

        } else {
          if (style == SINGLE) {
            
            if(dw_ && dim == iphi)
            {
                ngroup = atom->nfirst;
                for (i = 0; i < ngroup; i++)
                  if (decide_wedge(i,dim,lo,hi,ineed)) {
                    if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                    sendlist[iswap][nsend++] = i;
                  }
                for (i = atom->nlocal; i < nlast; i++)
                  if (decide_wedge(i,dim,lo,hi,ineed)) {
                    if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                    sendlist[iswap][nsend++] = i;
                  }
            }
            else 
            {
                ngroup = atom->nfirst;
                for (i = 0; i < ngroup; i++)
                  if (decide(i,dim,lo,hi,ineed)) {
                    if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                    sendlist[iswap][nsend++] = i;
                  }
                for (i = atom->nlocal; i < nlast; i++)
                  if (decide(i,dim,lo,hi,ineed)) {
                    if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                    sendlist[iswap][nsend++] = i;
                  }
            }
          } else {
            ngroup = atom->nfirst;
            for (i = 0; i < ngroup; i++) {
              const int itype = get_type(i);
              if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
              }
            }
            for (i = atom->nlocal; i < nlast; i++) {
              const int itype = get_type(i);
              if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
              }
            }
          }
        }
      }

      // pack up list of border atoms

      if (nsend*size_border > maxsend) grow_send(nsend*size_border,0);
      if (pb_)
        n = pb_->pack_comm(OPERATION_COMM_BORDERS,
                           nsend,sendlist[iswap],buf_send,
                           pbc_flag[iswap],pbc[iswap]);
      else {
        if (ghost_velocity)
          n = avec->pack_border_vel(nsend,sendlist[iswap],buf_send,
                                    pbc_flag[iswap],pbc[iswap]);
        else
          n = avec->pack_border(nsend,sendlist[iswap],buf_send,
                                pbc_flag[iswap],pbc[iswap]);
      }

      // swap atoms with other proc
      // no MPI calls except SendRecv if nsend/nrecv = 0
      // put incoming ghosts at end of my atom arrays
      // if swapping with self, simply copy, no messages

      if (sendproc[iswap] != me) {
        MPI_Sendrecv(&nsend,1,MPI_INT,sendproc[iswap],0,
                     &nrecv,1,MPI_INT,recvproc[iswap],0,world,
                     MPI_STATUS_IGNORE);
        if (nrecv*size_border > maxrecv) grow_recv(nrecv*size_border);
        if (nrecv) MPI_Irecv(buf_recv,nrecv*size_border,MPI_DOUBLE,
                             recvproc[iswap],0,world,&request);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
        if (nrecv) MPI_Wait(&request,MPI_STATUS_IGNORE);
        buf = buf_recv;
      } else {
        nrecv = nsend;
        buf = buf_send;
      }

      // unpack buffer

      if (pb_)
        pb_->unpack_comm(OPERATION_COMM_BORDERS,
                         nrecv, pb_->get_ntotal(), buf);
      else {
        if (ghost_velocity)
          avec->unpack_border_vel(nrecv,atom->nlocal+atom->nghost,buf);
        else
          avec->unpack_border(nrecv,atom->nlocal+atom->nghost,buf);
      }

      // set all pointers & counters

      smax = MAX(smax,nsend);
      rmax = MAX(rmax,nrecv);
      sendnum[iswap] = nsend;
      recvnum[iswap] = nrecv;
      size_forward_recv[iswap] = nrecv*size_forward;
      size_reverse_send[iswap] = nrecv*size_reverse;
      size_reverse_recv[iswap] = nsend*size_reverse;
      firstrecv[iswap] = pb_ ? pb_->get_ntotal() : atom->nlocal + atom->nghost;
      if (pb_)
        pb_->add_nghost(nrecv);
      else
        atom->nghost += nrecv;
      iswap++;
    }
  }

  // insure send/recv buffers are long enough for all forward & reverse comm

  int max = MAX(maxforward*smax,maxreverse*rmax);
  if (max > maxsend) grow_send(max,0);
  max = MAX(maxforward*rmax,maxreverse*smax);
  if (max > maxrecv) grow_recv(max);

  // reset global->local map

  if (pb_)
    pb_->generate_map();
  else if (map_style) atom->map_set();

  exchangeEventsCorrector();

}

/* ----------------------------------------------------------------------
   forward communication invoked by a Pair
------------------------------------------------------------------------- */

void CommBrick::forward_comm_pair(Pair *pair)
{
  int iswap,n;
  double *buf;
  MPI_Request request;

  for (iswap = 0; iswap < nswap; iswap++) {

    // pack buffer

    n = pair->pack_comm(sendnum[iswap],sendlist[iswap],
                        buf_send,pbc_flag[iswap],pbc[iswap]);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (recvnum[iswap])
        MPI_Irecv(buf_recv,n*recvnum[iswap],MPI_DOUBLE,
                  recvproc[iswap],0,world,&request);
      if (sendnum[iswap])
        MPI_Send(buf_send,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world);
      if (recvnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    pair->unpack_comm(recvnum[iswap],firstrecv[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   reverse communication invoked by a Pair
------------------------------------------------------------------------- */

void CommBrick::reverse_comm_pair(Pair *pair)
{
  int iswap,n;
  double *buf;
  MPI_Request request;

  for (iswap = nswap-1; iswap >= 0; iswap--) {

    // pack buffer

    n = pair->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (sendnum[iswap])
        MPI_Irecv(buf_recv,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,
                  world,&request);
      if (recvnum[iswap])
        MPI_Send(buf_send,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,world);
      if (sendnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    pair->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   forward communication invoked by a Fix
   n = constant number of datums per atom
------------------------------------------------------------------------- */

void CommBrick::forward_comm_fix(Fix *fix)
{
  int iswap,n;
  double *buf;
  MPI_Request request;

  for (iswap = 0; iswap < nswap; iswap++) {

    // pack buffer

    n = fix->pack_comm(sendnum[iswap],sendlist[iswap],
                       buf_send,pbc_flag[iswap],pbc[iswap]);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (recvnum[iswap])
        MPI_Irecv(buf_recv,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,
                  world,&request);
      if (sendnum[iswap])
        MPI_Send(buf_send,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world);
      if (recvnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    fix->unpack_comm(recvnum[iswap],firstrecv[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   reverse communication invoked by a Fix
   n = constant number of datums per atom
------------------------------------------------------------------------- */

void CommBrick::reverse_comm_fix(Fix *fix)
{
  int iswap,n;
  double *buf;
  MPI_Request request;

  for (iswap = nswap-1; iswap >= 0; iswap--) {

    // pack buffer

    n = fix->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (sendnum[iswap])
        MPI_Irecv(buf_recv,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,
                  world,&request);
      if (recvnum[iswap])
        MPI_Send(buf_send,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,world);
      if (sendnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    fix->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   forward communication invoked by a Fix
   n = total datums for all atoms, allows for variable number/atom
------------------------------------------------------------------------- */

void CommBrick::forward_comm_variable_fix(Fix *fix)
{
  int iswap,n;
  double *buf;
  MPI_Request request;
  MPI_Status status;

  for (iswap = 0; iswap < nswap; iswap++) {

    // pack buffer

    n = fix->pack_comm(sendnum[iswap],sendlist[iswap],
                       buf_send,pbc_flag[iswap],pbc[iswap]);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (recvnum[iswap])
        MPI_Irecv(buf_recv,maxrecv,MPI_DOUBLE,recvproc[iswap],0,
                  world,&request);
      if (sendnum[iswap])
        MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
      if (recvnum[iswap]) MPI_Wait(&request,&status);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    fix->unpack_comm(recvnum[iswap],firstrecv[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   forward communication invoked by a Compute
   nsize used only to set recv buffer limit
------------------------------------------------------------------------- */

void CommBrick::forward_comm_compute(Compute *compute)
{
  int iswap,n;
  double *buf;
  MPI_Request request;

  for (iswap = 0; iswap < nswap; iswap++) {

    // pack buffer

    n = compute->pack_comm(sendnum[iswap],sendlist[iswap],
                           buf_send,pbc_flag[iswap],pbc[iswap]);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (recvnum[iswap])
        MPI_Irecv(buf_recv,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,
                  world,&request);
      if (sendnum[iswap])
        MPI_Send(buf_send,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world);
      if (recvnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    compute->unpack_comm(recvnum[iswap],firstrecv[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   reverse communication invoked by a Compute
------------------------------------------------------------------------- */

void CommBrick::reverse_comm_compute(Compute *compute)
{
  int iswap,n;
  double *buf;
  MPI_Request request;

  for (iswap = nswap-1; iswap >= 0; iswap--) {

    // pack buffer

    n = compute->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (sendnum[iswap])
        MPI_Irecv(buf_recv,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,
                  world,&request);
      if (recvnum[iswap])
        MPI_Send(buf_send,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,world);
      if (sendnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    compute->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   forward communication invoked by a Dump
------------------------------------------------------------------------- */

void CommBrick::forward_comm_dump(Dump *dump)
{
  int iswap,n;
  double *buf;
  MPI_Request request;

  for (iswap = 0; iswap < nswap; iswap++) {

    // pack buffer

    n = dump->pack_comm(sendnum[iswap],sendlist[iswap],
                        buf_send,pbc_flag[iswap],pbc[iswap]);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (recvnum[iswap])
        MPI_Irecv(buf_recv,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,
                  world,&request);
      if (sendnum[iswap])
        MPI_Send(buf_send,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world);
      if (recvnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    dump->unpack_comm(recvnum[iswap],firstrecv[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   reverse communication invoked by a Dump
------------------------------------------------------------------------- */

void CommBrick::reverse_comm_dump(Dump *dump)
{
  int iswap,n;
  double *buf;
  MPI_Request request;

  for (iswap = nswap-1; iswap >= 0; iswap--) {

    // pack buffer

    n = dump->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send);

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (sendnum[iswap])
        MPI_Irecv(buf_recv,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,
                  world,&request);
      if (recvnum[iswap])
        MPI_Send(buf_send,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,world);
      if (sendnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    dump->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf);
  }
}

/* ----------------------------------------------------------------------
   forward communication of N values in per-atom array
------------------------------------------------------------------------- */

void CommBrick::forward_comm_array(int nsize, double **array)
{
  int i,j,k,m,iswap,last;
  double *buf;
  MPI_Request request;

  // insure send/recv bufs are big enough for nsize
  // based on smax/rmax from most recent borders() invocation

  if (nsize > maxforward) {
    maxforward = nsize;
    if (maxforward*smax > maxsend) grow_send(maxforward*smax,0);
    if (maxforward*rmax > maxrecv) grow_recv(maxforward*rmax);
  }

  for (iswap = 0; iswap < nswap; iswap++) {

    // pack buffer

    m = 0;
    for (i = 0; i < sendnum[iswap]; i++) {
      j = sendlist[iswap][i];
      for (k = 0; k < nsize; k++)
        buf_send[m++] = array[j][k];
    }

    // exchange with another proc
    // if self, set recv buffer to send buffer

    if (sendproc[iswap] != me) {
      if (recvnum[iswap])
        MPI_Irecv(buf_recv,nsize*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,
                  world,&request);
      if (sendnum[iswap])
        MPI_Send(buf_send,nsize*sendnum[iswap],MPI_DOUBLE,
                 sendproc[iswap],0,world);
      if (recvnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      buf = buf_recv;
    } else buf = buf_send;

    // unpack buffer

    m = 0;
    last = firstrecv[iswap] + recvnum[iswap];
    for (i = firstrecv[iswap]; i < last; i++)
      for (k = 0; k < nsize; k++)
        array[i][k] = buf[m++];
  }
}

/* ----------------------------------------------------------------------
   realloc the size of the iswap sendlist as needed with BUFFACTOR
------------------------------------------------------------------------- */

void CommBrick::grow_list(int iswap, int n)
{
  maxsendlist[iswap] = static_cast<int> (BUFFACTOR * n);
  memory->grow(sendlist[iswap],maxsendlist[iswap],"comm:sendlist[iswap]");
}

/* ----------------------------------------------------------------------
   realloc the buffers needed for swaps
------------------------------------------------------------------------- */

void CommBrick::grow_swap(int n)
{
  free_swap();
  allocate_swap(n);
  if (mode == MULTI) {
    free_multi();
    allocate_multi(n);
  }

  sendlist = (int **)
    memory->srealloc(sendlist,n*sizeof(int *),"comm:sendlist");
  memory->grow(maxsendlist,n,"comm:maxsendlist");
  for (int i = maxswap; i < n; i++) {
    maxsendlist[i] = BUFMIN;
    memory->create(sendlist[i],BUFMIN,"comm:sendlist[i]");
  }
  maxswap = n;
}

/* ----------------------------------------------------------------------
   allocation of swap info
------------------------------------------------------------------------- */

void CommBrick::allocate_swap(int n)
{
  memory->create(sendnum,n,"comm:sendnum");
  memory->create(recvnum,n,"comm:recvnum");
  memory->create(sendproc,n,"comm:sendproc");
  memory->create(recvproc,n,"comm:recvproc");
  memory->create(size_forward_recv,n,"comm:size");
  memory->create(size_reverse_send,n,"comm:size");
  memory->create(size_reverse_recv,n,"comm:size");
  memory->create(slablo,n,"comm:slablo");
  memory->create(slabhi,n,"comm:slabhi");
  memory->create(firstrecv,n,"comm:firstrecv");
  memory->create(pbc_flag,n,"comm:pbc_flag");
  memory->create(pbc,n,6,"comm:pbc");
}

/* ----------------------------------------------------------------------
   allocation of multi-type swap info
------------------------------------------------------------------------- */

void CommBrick::allocate_multi(int n)
{
  multilo = memory->create(multilo,n,atom->ntypes+1,"comm:multilo");
  multihi = memory->create(multihi,n,atom->ntypes+1,"comm:multihi");
}

/* ----------------------------------------------------------------------
   free memory for swaps
------------------------------------------------------------------------- */

void CommBrick::free_swap()
{
  memory->destroy(sendnum);
  memory->destroy(recvnum);
  memory->destroy(sendproc);
  memory->destroy(recvproc);
  memory->destroy(size_forward_recv);
  memory->destroy(size_reverse_send);
  memory->destroy(size_reverse_recv);
  memory->destroy(slablo);
  memory->destroy(slabhi);
  memory->destroy(firstrecv);
  memory->destroy(pbc_flag);
  memory->destroy(pbc);
}

/* ----------------------------------------------------------------------
   free memory for multi-type swaps
------------------------------------------------------------------------- */

void CommBrick::free_multi()
{
  memory->destroy(multilo);
  memory->destroy(multihi);
  multilo = multihi = NULL;
}

/* ----------------------------------------------------------------------
   return # of bytes of allocated memory
------------------------------------------------------------------------- */

bigint CommBrick::memory_usage()
{
  bigint bytes = 0;
  bytes += nprocs * sizeof(int);    // grid2proc
  for (int i = 0; i < nswap; i++)
    bytes += memory->usage(sendlist[i],maxsendlist[i]);
  bytes += memory->usage(buf_send,maxsend+bufextra);
  bytes += memory->usage(buf_recv,maxrecv);
  return bytes;
}

inline bool CommBrick::use_gran_opt()
{
    const bool has_radius = pb_ ? pb_->has_radius() : atom->radius != NULL;
    return (domain->triclinic == 0 && has_radius);
}

/* ----------------------------------------------------------------------
   decide if border element, optimization for granular
------------------------------------------------------------------------- */

inline bool CommBrick::decide(int i,int dim,double lo,double hi,int ineed)
{
    double **x = atom->x;
    const double pos = pb_ ? pb_->get_pos(i,dim) : x[i][dim];
    const double radius = pb_ ? pb_->get_radius(i) : atom->radius[i];
    const double offset = use_gran_opt() ? radius : 0.;

    if( ((ineed % 2 == 0) && pos >= lo && pos <= (hi + offset) ) ||
        ((ineed % 2 == 1) && pos >= (lo - offset) && pos <= hi )   )
        return true;

    return false;
}

/* ----------------------------------------------------------------------
   decide if border element for wedge case, optimization for granular
------------------------------------------------------------------------- */

inline bool CommBrick::decide_wedge(int i,int dim,double lo,double hi,int ineed)
{
    double **x = atom->x;
    double coo[2],d[2];
    coo[0] = pb_ ? pb_->get_pos(i, iphi) : x[i][iphi];
    coo[1] = pb_ ? pb_->get_pos(i, (iphi+1)%3) : x[i][(iphi+1)%3];
    const double radius = pb_ ? pb_->get_radius(i) : atom->radius[i];
    const double offset = use_gran_opt() ? radius : 0.;

    if (ineed % 2 == 0)
    {
        vectorSubtract2D(coo,pleft,d);
        if(vectorDot2D(d,nleft) >= -offset)
        {
            
            return true;
        }
    }
    
    else if (ineed % 2 == 1)
    {
        vectorSubtract2D(coo,pright,d);
        if(vectorDot2D(d,nright) >= -offset)
        {
            
            return true;
        }
    }
    
    return false;
}

/* ----------------------------------------------------------------------
   routine to determine Exchange Events
   similar to CommBrick:Exchange()
------------------------------------------------------------------------- */

void CommBrick::exchangeEventsRecorder()
{
    if(!exchangeEvents || nprocs==1 )
        return;

    exchangeEventsLocalId.clear();
    exchangeEventsReceivingProcess.clear();
    exchangeEventsGlobalProblemIds.clear();

    bool verbose = false; //true; //Developer to set here if needed for debugging
    int i;
    double **x;
    int    *tag;
    double *sublo,*subhi;
    double subloNeigh[3],subhiNeigh[3]; //boarders of 1st neighbor process (0=left if two), in each dim

    MPI_Status  status;
    MPI_Request request;

    //determine process boundaries
    //since called in CommBrick::Exchange, triclinic boundaries are already correct
    if (domain->triclinic == 0) {
      sublo = domain->sublo;
      subhi = domain->subhi;
    } else {
      sublo = domain->sublo_lamda;
      subhi = domain->subhi_lamda;
    }

    //Exchange info on sublo/subhi
    for(int dim=0;dim<3;dim++)
    {
       MPI_Send (&(sublo[dim]),      1, MPI_DOUBLE, procneigh[dim][1],0,world);          //to right
       MPI_Irecv(&(subloNeigh[dim]), 1, MPI_DOUBLE, procneigh[dim][0],0,world,&request); //from left
       MPI_Wait (&request,&status); //wait to receive

       MPI_Send (&(subhi[dim]),      1, MPI_DOUBLE, procneigh[dim][1],0,world);          //to right
       MPI_Irecv(&(subhiNeigh[dim]), 1, MPI_DOUBLE, procneigh[dim][0],0,world,&request); //from left
       MPI_Wait (&request,&status); //wait to receive
    }
    if(verbose)
       printf("[%d/%d]:CommBrick::exchangeEventsRecorder(); lo/hi: %g %g %g / %g %g %g, exchanging sublo/hi information with process %d %d %d: %g %g %g / %g %g %g \n",
               me, nprocs,
               sublo[0],sublo[1],sublo[2],subhi[0],subhi[1],subhi[2],
               procneigh[0][0], procneigh[1][0], procneigh[2][0],
               subloNeigh[0],subloNeigh[1],subloNeigh[2],
               subhiNeigh[0],subhiNeigh[1],subhiNeigh[2]);

    // loop over dimensions,
    x   = atom->x;
    tag = atom->tag;
    i = 0;
    while (i < atom->nlocal)
    {
        bool  willExchange    = false;
        bool  problemDetected = false;
//        if(verbose)
//           printf("[%d/%d]:CommBrick::exchangeEventsRecorder() looping  atom %d with pos: %g %g %g\n",
//                       me, nprocs, i, x[i][0], x[i][1], x[i][2]);

        for (int dim = 0; dim < 3; dim++)
        {
            if( procgrid[dim] == 1 ) //nothing to do if only one processor in this direction
              continue;

            double lo = sublo[dim];
            double hi = subhi[dim];
            double position = x[i][dim];
            bool   willExchangeThisDim = false;

            if(verbose)
              printf("[%d/%d]:CommBrick::exchangeEventsRecorder() looping dim: %d, atom %d with pos: %g\n",
                       me, nprocs, dim, i, position);

            if (position < lo || position >= hi)
            {
                if(willExchange) //exchange over multiple dims, has been already registered!
                {
                    exchangeEventsReceivingProcess.back() = -1; //reset receiving processor
                    problemDetected = true;
                    continue;
                }

                //Append results from this direction to overall event list
                //can have only one per local ID!
                willExchange        = true;
                willExchangeThisDim = true;
                exchangeEventsLocalId.push_back(i);
                exchangeEventsReceivingProcess.push_back(procneigh[dim][0]); //by default, use first neighbor, correct later

                if(verbose)
                printf("[%d/%d]:CommBrick::exchangeEventsRecorder() recording exchange for atom %d with pos: %g\n",
                       me, nprocs, i, position );
            }

            if( (procgrid[dim]>2) && willExchangeThisDim )
            {
                if(verbose)
                    printf("[%d/%d]: checking dim: %d because procgrid[dim] = %d \n",
                       me, nprocs, dim, procgrid[dim]);

                if( position < subloNeigh[dim] || position >= subhiNeigh[dim] ) //not on left neighbor
                {
                    exchangeEventsReceivingProcess.back() = procneigh[dim][1];
                    if(verbose)
                        printf("[%d/%d]: checking dim: %d, detected transfer to right process with id %d \n",
                               me, nprocs, dim, exchangeEventsReceivingProcess.back());
                }
            } //handle more than 2 processors in this dimension
        } //loop over dimensions

        //if problem detected, regist globalId for later determination of receiving process
        if(problemDetected)
           exchangeEventsGlobalProblemIds.push_back(tag[i]);

        i++;

    } //loop over atoms
}

/* ----------------------------------------------------------------------
   routine to determine the id of the process a particle has been
   transferred to in case of multi process dimensions crossing
------------------------------------------------------------------------- */

void CommBrick::exchangeEventsCorrector()
{
    if( !exchangeEvents || nprocs==1 )
        return;

    //Determine global count of problemIds, exit if none
    bool verbose = false; //true; //Developer to set here if needed for debugging
    int  global_sum = 0;
    int  local_sum = exchangeEventsGlobalProblemIds.size();
    MPI_Barrier(world);
    MPI_Allreduce(&local_sum, &global_sum,
                  1, MPI_INT, MPI_SUM, world
                 );
    if( global_sum==0 )
        return;

    //Collect information on problem IDs
    int * exchangeReceiveCounts;
    int * exchangeReceiveDisplacement;
    int * exchangeGlobalProblems;
    int * exchangeOwningProcess;
    memory->create(exchangeReceiveCounts,nprocs,"comm:exchangeReceiveCounts"); //TODO: do just once
    memory->create(exchangeReceiveDisplacement,nprocs,"comm:exchangeReceiveDisplacement");  //TODO
    memory->create(exchangeGlobalProblems,global_sum,"comm:exchangeGlobalProblems");
    memory->create(exchangeOwningProcess, global_sum,"comm:exchangeOwningProcess");

    MPI_Allgather(&local_sum,             1, MPI_INT,
                   exchangeReceiveCounts, 1, MPI_INT,
                   world
                 );

    exchangeReceiveDisplacement[0] = 0;
    for(int iPro=1; iPro<nprocs; iPro++)
        exchangeReceiveDisplacement[iPro] = exchangeReceiveCounts[iPro-1]
                                          + exchangeReceiveDisplacement[iPro-1];

    //Inform all CPUs about problem IDs
    MPI_Barrier(world);
    MPI_Allgatherv(&(exchangeEventsGlobalProblemIds.front()), local_sum, MPI_INT,
                   exchangeGlobalProblems, exchangeReceiveCounts, exchangeReceiveDisplacement,
                   MPI_INT, world
                  );

    if(verbose && me==0)
    {
        printf("**exchangeEventsCorrector: globalProblems %d with IDs: ",
                global_sum
              );
        for(int k=0;k<global_sum;k++)
          printf(" %d, ",exchangeGlobalProblems[k] );

        printf("; exchangeReceiveCounts:");

        for(int k=0;k<nprocs;k++)
          printf(" %d, ",
                exchangeReceiveCounts[k]
              );
        printf("\n");
    }

    //Loop through list of global problem IDs and detect if CPU owns this atom
    for(int iter=0; iter<global_sum; iter++)
    {
        exchangeOwningProcess[iter] = -1;
        int currLocal     = atom->map(exchangeGlobalProblems[iter]);
        if( (currLocal<0) || (currLocal >= atom->nlocal) ) //not owned
            continue;
        else
            exchangeOwningProcess[iter] = me;
    }

    MPI_Barrier(world);
    MPI_Allreduce(exchangeOwningProcess,
                  exchangeGlobalProblems, //re-use this container, is now (received) owning process id
                  global_sum, MPI_INT, MPI_MAX, world
                 );

    //Check if all problem IDs have been detected
    for(int j=0;j<global_sum;j++)
    {
       if(verbose && me==0)
           printf("**exchangeEventsCorrector: receivedOwningProcess[%d]: %d \n",
                   j, exchangeGlobalProblems[j]
                 );

       if(exchangeGlobalProblems[j]<0)
            error->one(FLERR,"CommBrick::exchangeEventsCorrector: Could not find particle. Must have left the domain. You must upgrade this check in order to handle this situation.");
    }

    //Fill in the process ids
    unsigned int checkCounter=0;
    for(int iter =  exchangeReceiveDisplacement[me];  //loop global list to get correct iter
            iter < (exchangeReceiveDisplacement[me]+exchangeReceiveCounts[me]);
            iter++
       )
       for(unsigned int j=0; j<exchangeEventsReceivingProcess.size(); j++)
         if(exchangeEventsReceivingProcess[j]==-1) //have invalid receiving process id
         {
            exchangeEventsReceivingProcess[j]=exchangeGlobalProblems[iter];
            checkCounter++;
         }

    if( checkCounter!=exchangeEventsGlobalProblemIds.size() )
        error->all(FLERR,"CommBrick::exchangeEventsCorrector: Problem when fill in corrected process ids! This is fatal.");

    //Clear memory
    memory->destroy(exchangeReceiveCounts);
    memory->destroy(exchangeReceiveDisplacement);
    memory->destroy(exchangeGlobalProblems);
    memory->destroy(exchangeOwningProcess);
}