/*----------------------------------------------------------------------
  PuReMD - Purdue ReaxFF Molecular Dynamics Program

  Copyright (2010) Purdue University
  Hasan Metin Aktulga, haktulga@cs.purdue.edu
  Joseph Fogarty, jcfogart@mail.usf.edu
  Sagar Pandit, pandit@usf.edu
  Ananth Y Grama, ayg@cs.purdue.edu

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of
  the License, or (at your option) any later version.

  This program 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.
  See the GNU General Public License for more details:
  <http://www.gnu.org/licenses/>.
  ----------------------------------------------------------------------*/

#include "reax_types.h"

#include "comm_tools.h"

#include "grid.h"
#include "reset_tools.h"
#include "tool_box.h"
#include "vector.h"


#if defined(NEUTRAL_TERRITORY)
void Setup_NT_Comm( reax_system * const system, control_params * const control,
                 mpi_datatypes * const mpi_data )
{
    int i, d;
    real bndry_cut;
    neighbor_proc *nbr_pr;
    simulation_box *my_box;
    ivec nbr_coords, nbr_recv_coords;
    ivec r[12] = {
        {0, 0, -1}, // -z
        {0, 0, +1}, // +z
        {0, -1, 0}, // -y
        {-1, -1, 0}, // -x-y
        {-1, 0, 0}, // -x
        {-1, +1, 0},  // -x+y

        {0, 0, +1}, // +z
        {0, 0, -1}, // -z
        {0, +1, 0}, // +y
        {+1, +1, 0}, // +x+y
        {+1, 0, 0}, // +x
        {+1, -1, 0}  // +x-y
    };
    my_box = &system->my_box;
    bndry_cut = system->bndry_cuts.ghost_cutoff;
    system->num_nt_nbrs = MAX_NT_NBRS;

    /* identify my neighbors */
    for ( i = 0; i < system->num_nt_nbrs; ++i )
    {
        nbr_pr = &system->my_nt_nbrs[i];
        ivec_Sum( nbr_coords, system->my_coords, r[i] ); /* actual nbr coords */
        MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_coords, &nbr_pr->rank );
        
        /* set the rank of the neighbor processor in the receiving direction */
        ivec_Sum( nbr_recv_coords, system->my_coords, r[i + 6] ); /* actual nbr coords */
        MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_recv_coords, &nbr_pr->receive_rank );

        for ( d = 0; d < 3; ++d )
        {
            /* determine the boundary area with this nbr */
            if ( r[i][d] < 0 )
            {
                nbr_pr->bndry_min[d] = my_box->min[d];
                nbr_pr->bndry_max[d] = my_box->min[d] + bndry_cut;
            }
            else if ( r[i][d] > 0 )
            {
                nbr_pr->bndry_min[d] = my_box->max[d] - bndry_cut;
                nbr_pr->bndry_max[d] = my_box->max[d];
            }
            else
            {
                nbr_pr->bndry_min[d] = my_box->min[d];
                nbr_pr->bndry_max[d] = my_box->max[d];
            }

            /* determine if it is a periodic neighbor */
            if ( nbr_coords[d] < 0 )
            {
                nbr_pr->prdc[d] = -1;
            }
            else if ( nbr_coords[d] >= control->procs_by_dim[d] )
            {
                nbr_pr->prdc[d] = 1;
            }
            else
            {
                nbr_pr->prdc[d] = 0;
            }
        }

    }
}


static int Sort_Neutral_Territory( reax_system *system, int dir, mpi_out_data *out_bufs, int write )
{
    int i, cnt;
    reax_atom *atoms;
    neighbor_proc *nbr_pr;

    cnt = 0;
    atoms = system->my_atoms;
    /* place each atom into the appropriate outgoing list */
    nbr_pr = &system->my_nt_nbrs[dir];

    for ( i = 0; i < system->n; ++i )
    {
        if ( nbr_pr->bndry_min[0] <= atoms[i].x[0]
                && atoms[i].x[0] < nbr_pr->bndry_max[0]
                && nbr_pr->bndry_min[1] <= atoms[i].x[1]
                && atoms[i].x[1] < nbr_pr->bndry_max[1]
                && nbr_pr->bndry_min[2] <= atoms[i].x[2]
                && atoms[i].x[2] < nbr_pr->bndry_max[2] )
        {
            if ( write )
            {
                out_bufs[dir].index[out_bufs[dir].cnt] = i;
                out_bufs[dir].cnt++;
            }
            else
            {
                cnt++;
            }
        }
    }

    return cnt;
}


void Estimate_NT_Atoms( reax_system * const system, mpi_datatypes * const mpi_data )
{
    int d;
    mpi_out_data *out_bufs;
    neighbor_proc *nbr;

    out_bufs = mpi_data->out_nt_buffers;

    for ( d = 0; d < 6; ++d )
    {
        /* count the number of atoms in each processor's outgoing list */
        nbr = &system->my_nt_nbrs[d];
        nbr->est_send = Sort_Neutral_Territory( system, d, out_bufs, 0 );

        /* estimate the space needed based on the count above */
        nbr->est_send = MAX( MIN_SEND, nbr->est_send * SAFER_ZONE_NT );

        /* allocate the estimated space */
        out_bufs[d].index = scalloc( 2 * nbr->est_send, sizeof(int),
                "Estimate_NT_Atoms::out_bufs[d].index", MPI_COMM_WORLD );
        out_bufs[d].out_atoms = scalloc( 2 * nbr->est_send, sizeof(real),
                "Estimate_NT_Atoms::out_bufs[d].out_atoms", MPI_COMM_WORLD );

        /* sort the atoms to their outgoing buffers */
        // TODO: to call or not to call?
        //Sort_Neutral_Territory( system, d, out_bufs, 1 );
    }
}
#endif


void Check_MPI_Error( int code, const char * const filename, int line )
{
    int len;
    char err_msg[MPI_MAX_ERROR_STRING];

    if ( code != MPI_SUCCESS )
    {
        MPI_Error_string( code, err_msg, &len );

        fprintf( stderr, "[ERROR] MPI error\n" );
        fprintf( stderr, "    [INFO] At line %d in file %.*s\n",
                line, (int) strnlen(filename, MAX_STR), filename );
        fprintf( stderr, "    [INFO] Error code %d\n", code );
        fprintf( stderr, "    [INFO] Error message: %.*s\n", len, err_msg );
        MPI_Abort( MPI_COMM_WORLD, RUNTIME_ERROR );
    }
}


void Setup_Comm( reax_system * const system, control_params * const control,
        mpi_datatypes * const mpi_data )
{
    int i, d;
    const real bndry_cut = system->bndry_cuts.ghost_cutoff;
    neighbor_proc *nbr_pr;
    simulation_box * const my_box = &system->my_box;
    ivec nbr_coords;
    ivec r[6] = {
        { -1, 0, 0}, { 1, 0, 0}, // -x, +x
        {0, -1, 0}, {0, 1, 0}, // -y, +y
        {0, 0, -1}, {0, 0, 1}, // -z, +z
    };

    /* identify my neighbors */
    system->num_nbrs = MAX_NBRS;
    for ( i = 0; i < system->num_nbrs; ++i )
    {
        ivec_Sum( nbr_coords, system->my_coords, r[i] ); /* actual nbr coords */
        nbr_pr = &system->my_nbrs[i];
        MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_coords, &nbr_pr->rank );

        for ( d = 0; d < 3; ++d )
        {
            /* determine the boundary area with this nbr */
            if ( r[i][d] < 0 )
            {
                nbr_pr->bndry_min[d] = my_box->min[d];
                nbr_pr->bndry_max[d] = my_box->min[d] + bndry_cut;
            }
            else if ( r[i][d] > 0 )
            {
                nbr_pr->bndry_min[d] = my_box->max[d] - bndry_cut;
                nbr_pr->bndry_max[d] = my_box->max[d];
            }
            else
            {
                nbr_pr->bndry_min[d] = NEG_INF;
                nbr_pr->bndry_max[d] = NEG_INF;
            }

            /* determine if it is a periodic neighbor */
            if ( nbr_coords[d] < 0 )
            {
                nbr_pr->prdc[d] = -1;
            }
            else if ( nbr_coords[d] >= control->procs_by_dim[d] )
            {
                nbr_pr->prdc[d] = 1;
            }
            else
            {
                nbr_pr->prdc[d] = 0;
            }
        }
    }
}


void Update_Comm( reax_system * const system )
{
    int i, d;
    const real bndry_cut = system->bndry_cuts.ghost_cutoff;
    neighbor_proc *nbr_pr;
    simulation_box * const my_box = &system->my_box;
    ivec r[6] = {
        { -1, 0, 0}, { 1, 0, 0}, // -x, +x
        {0, -1, 0}, {0, 1, 0}, // -y, +y
        {0, 0, -1}, {0, 0, 1}, // -z, +z
    };

    /* identify my neighbors */
    for ( i = 0; i < system->num_nbrs; ++i )
    {
        nbr_pr = &system->my_nbrs[i];

        for ( d = 0; d < 3; ++d )
        {
            /* determine the boundary area with this nbr */
            if ( r[i][d] < 0 )
            {
                nbr_pr->bndry_min[d] = my_box->min[d];
                nbr_pr->bndry_max[d] = my_box->min[d] + bndry_cut;
            }
            else if ( r[i][d] > 0 )
            {
                nbr_pr->bndry_min[d] = my_box->max[d] - bndry_cut;
                nbr_pr->bndry_max[d] = my_box->max[d];
            }
            else
            {
                nbr_pr->bndry_min[d] = NEG_INF;
                nbr_pr->bndry_max[d] = NEG_INF;
            }
        }
    }
}


/********************* ATOM TRANSFER ***********************/

/***************** PACK & UNPACK ATOMS *********************/
static void Pack_MPI_Atom( mpi_atom * const matm, const reax_atom * const ratm, int i )
{
    matm->orig_id = ratm->orig_id;
    matm->imprt_id = i;
    matm->type = ratm->type;
    matm->num_bonds = ratm->num_bonds;
    matm->num_hbonds = ratm->num_hbonds;
    strncpy( matm->name, ratm->name, sizeof(matm->name) - 1 );
    matm->name[sizeof(matm->name) - 1] = '\0';
    rvec_Copy( matm->x, ratm->x );
    rvec_Copy( matm->v, ratm->v );
    rvec_Copy( matm->f_old, ratm->f_old );
    memcpy( matm->s, ratm->s, sizeof(rvec4) ); //rvec_Copy( matm->s, ratm->s );
    memcpy( matm->t, ratm->t, sizeof(rvec4) ); //rvec_Copy( matm->t, ratm->t );
}


static void Unpack_MPI_Atom( reax_atom * const ratm, const mpi_atom * const matm )
{
    ratm->orig_id = matm->orig_id;
    ratm->imprt_id = matm->imprt_id;
    ratm->type = matm->type;
    ratm->num_bonds = matm->num_bonds;
    ratm->num_hbonds = matm->num_hbonds;
    strncpy( ratm->name, matm->name, sizeof(ratm->name) - 1 );
    ratm->name[sizeof(ratm->name) - 1] = '\0';
    rvec_Copy( ratm->x, matm->x );
    rvec_Copy( ratm->v, matm->v );
    rvec_Copy( ratm->f_old, matm->f_old );
    memcpy( ratm->s, matm->s, sizeof(rvec4) ); //rvec_Copy( ratm->s, matm->s );
    memcpy( ratm->t, matm->t, sizeof(rvec4) ); //rvec_Copy( ratm->t, matm->t );
}


/* Count number of atoms to be placed in egress MPI buffers.
 * The ownership of these atoms is to be transferred to other MPI processes. */
static void Count_Transfer_Atoms( reax_system const * const system,
        int start, int end, int dim, mpi_out_data * const out_bufs,
        int *cnt1, int *cnt2  )
{
    int i, d;
    simulation_box *my_box;

    my_box = &system->my_box;

    /* count atoms in the appropriate outgoing lists */
    for ( i = start; i < end; ++i )
    {
        for ( d = dim; d < 3; ++d )
        {
            if ( system->my_atoms[i].x[d] < my_box->min[d] )
            {
                ++(*cnt1);
                break;
            }
            else if ( system->my_atoms[i].x[d] >= my_box->max[d] )
            {
                ++(*cnt2);
                break;
            }
        }
    }
}


/* Populate egress MPI buffers with atoms (of reax_atom type) whose
 * ownership is being transferred to other MPI processes
 *
 * Note: the "sort" refers to dividing the atoms into their appropriate buffers */
static void Sort_Transfer_Atoms( reax_system * const system, int start, int end,
        int dim, mpi_out_data * const out_bufs, mpi_datatypes * const mpi_data )
{
    int i, d, out_cnt;
    simulation_box * const my_box = &system->my_box;
    mpi_atom *out_buf;

    /* place each atom into the appropriate outgoing list */
    for ( i = start; i < end; ++i )
    {
        for ( d = dim; d < 3; ++d )
        {
            if ( system->my_atoms[i].x[d] < my_box->min[d] )
            {
                out_cnt = out_bufs[2 * d].cnt++;
                out_buf = out_bufs[2 * d].out_atoms;
                Pack_MPI_Atom( &out_buf[out_cnt], &system->my_atoms[i], i );
                system->my_atoms[i].orig_id = -1;
                break;
            }
            else if ( system->my_atoms[i].x[d] >= my_box->max[d] )
            {
                out_cnt = out_bufs[2 * d + 1].cnt++;
                out_buf = out_bufs[2 * d + 1].out_atoms;
                Pack_MPI_Atom( &out_buf[out_cnt], &system->my_atoms[i], i );
                system->my_atoms[i].orig_id = -1;
                break;
            }
        }
    }
}


/*********************** UNPACKER **************************/
static void Unpack_Transfer_Message( reax_system * const system,
        int end, void * const mpi_buffer,
        int cnt, neighbor_proc * const nbr, int dim )
{
    int i;
    real dx;
    mpi_atom * const src = (mpi_atom *) mpi_buffer;

    for ( i = 0; i < cnt; ++i )
    {
        Unpack_MPI_Atom( &system->my_atoms[end + i], &src[i] );
    }

    /* adjust coordinates of recved atoms if nbr is a periodic one */
    if ( nbr->prdc[dim] )
    {
        dx = nbr->prdc[dim] * system->big_box.box_norms[dim];

        for ( i = 0; i < cnt; ++i )
        {
            system->my_atoms[end + i].x[dim] += dx;
        }
    }
}


/************** EXCHANGE BOUNDARY ATOMS *****************/

/************ PACK & UNPACK BOUNDARY ATOMS **************/
static void Pack_Boundary_Atom( boundary_atom * const matm,
        const reax_atom * const ratm, int i )
{
    matm->orig_id = ratm->orig_id;
    matm->imprt_id = i;
    matm->type = ratm->type;
    matm->num_bonds = ratm->num_bonds;
    matm->num_hbonds = ratm->num_hbonds;
    rvec_Copy( matm->x, ratm->x );
}


static void Unpack_Boundary_Atom( reax_atom * const ratm,
        const boundary_atom * const matm )
{
    ratm->orig_id = matm->orig_id;
    ratm->imprt_id = matm->imprt_id;
    ratm->type = matm->type;
    ratm->num_bonds = matm->num_bonds;
    ratm->num_hbonds = matm->num_hbonds;
    //ratm->renumber = offset + matm->imprt_id;
    rvec_Copy( ratm->x, matm->x );
}


void Count_Boundary_Atoms( reax_system const * const system,
        int start, int end, int dim, mpi_out_data * const out_bufs,
        int *cnt1, int *cnt2 )
{
    int i;

    /* place each atom into the appropriate outgoing list */
    for ( i = 0; i < end; ++i )
    {
        if ( system->my_nbrs[2 * dim].bndry_min[dim] <= system->my_atoms[i].x[dim]
                && system->my_atoms[i].x[dim] < system->my_nbrs[2 * dim].bndry_max[dim] )
        {
            ++(*cnt1);
        }

        if ( system->my_nbrs[2 * dim + 1].bndry_min[dim] <= system->my_atoms[i].x[dim]
                && system->my_atoms[i].x[dim] < system->my_nbrs[2 * dim + 1].bndry_max[dim] )
        {
            ++(*cnt2);
        }
    }
}


void Sort_Boundary_Atoms( reax_system * const system, int start, int end,
        int dim, mpi_out_data * const out_bufs, mpi_datatypes * const mpi_data )
{
    int i, p, out_cnt;
    boundary_atom *out_buf;
    neighbor_proc *nbr_pr;

    /* place each atom into the appropriate egress buffer */
    for ( i = 0; i < end; ++i )
    {
        for ( p = 2 * dim; p < 2 * dim + 2; ++p )
        {
            nbr_pr = &system->my_nbrs[p];
            out_buf = (boundary_atom *) out_bufs[p].out_atoms;

            if ( nbr_pr->bndry_min[dim] <= system->my_atoms[i].x[dim]
                    && system->my_atoms[i].x[dim] < nbr_pr->bndry_max[dim] )
            {
                out_cnt = out_bufs[p].cnt++;
                out_bufs[p].index[out_cnt] = i;
                Pack_Boundary_Atom( &out_buf[out_cnt], &system->my_atoms[i], i );
            }
        }
    }
}


/* Copy received atoms out of MPI ingress buffer */
void Unpack_Exchange_Message( reax_system * const system, int end,
        void * const mpi_buffer, int cnt,
        neighbor_proc * const nbr, int dim )
{
    int i;
    real dx;
    const boundary_atom * const src = (boundary_atom *) mpi_buffer;

    if ( end + cnt > system->N )
    {
        system->my_atoms = srealloc( system->my_atoms,
                sizeof(reax_atom) * (end + cnt), "Unpack_Exchange_Message" );
        system->N = end + cnt;
    }

    if ( nbr->prdc[dim] )
    {
        /* adjust coordinates of recved atoms if nbr is a periodic one */
        dx = nbr->prdc[dim] * system->big_box.box_norms[dim];

        for ( i = 0; i < cnt; ++i )
        {
            Unpack_Boundary_Atom( &system->my_atoms[end + i], &src[i] );
            system->my_atoms[end + i].x[dim] += dx;
        }
    }
    else
    {
        for ( i = 0; i < cnt; ++i )
        {
            Unpack_Boundary_Atom( &system->my_atoms[end + i], &src[i] );
        }
    }

    /* record the atoms recv'd from this nbr */
    nbr->atoms_str = end;
    nbr->atoms_cnt = cnt;
}


static void Count_Position_Updates( reax_system const * const system,
        int start, int end, int dim, mpi_out_data * const out_bufs,
        int *cnt1, int *cnt2 )
{
    *cnt1 = out_bufs[2 * dim].cnt;
    *cnt2 = out_bufs[2 * dim + 1].cnt;
}


static void Sort_Position_Updates( reax_system * const system, int start, int end,
        int dim, mpi_out_data * const out_bufs, mpi_datatypes * const mpi_data )
{
    int i, p;
    rvec *out;

    for ( p = 2 * dim; p < 2 * dim + 2; ++p )
    {
        out = (rvec*) out_bufs[p].out_atoms;

        for ( i = 0; i < out_bufs[p].cnt; ++i )
        {
            rvec_Copy( out[i], system->my_atoms[ out_bufs[p].index[i] ].x );
        }
    }
}


static void Unpack_Position_Updates( reax_system * const system, int end,
        void * const mpi_buffer, int cnt, neighbor_proc * const nbr, int dim )
{
    int i;
    const int start = nbr->atoms_str;
    real dx;
    rvec * const src = (rvec*) mpi_buffer;

    if ( nbr->prdc[dim] )
    {
        /* adjust coordinates of recved atoms if nbr is a periodic one */
        dx = nbr->prdc[dim] * system->big_box.box_norms[dim];

        for ( i = 0; i < cnt; ++i )
        {
            rvec_Copy( system->my_atoms[start + i].x, src[i] );
            system->my_atoms[start + i].x[dim] += dx;
        }
    }
    else
    {
        for ( i = 0; i < cnt; ++i )
        {
            rvec_Copy( system->my_atoms[start + i].x, src[i] );
        }
    }
}


int SendRecv( reax_system * const system, mpi_datatypes * const mpi_data,
        MPI_Datatype type, message_counter count_func,
        message_sorter sort_func, unpacker unpack, int clr )
{
    int d, cnt1, cnt2, start, end, ret;
    mpi_out_data * const out_bufs = mpi_data->out_buffers;
    const MPI_Comm comm = mpi_data->comm_mesh3D;
    MPI_Request req1, req2;
    MPI_Status stat1, stat2;
    neighbor_proc *nbr1, *nbr2;
    MPI_Aint extent, lower_bound;

    MPI_Type_get_extent( type, &lower_bound, &extent );

    if ( clr == TRUE )
    {
        Reset_Out_Buffers( mpi_data->out_buffers, system->num_nbrs );
    }
    start = 0;
    end = system->n;

    for ( d = 0; d < 3; ++d )
    {
        cnt1 = 0;
        cnt2 = 0;
        nbr1 = &system->my_nbrs[2 * d];
        nbr2 = &system->my_nbrs[2 * d + 1];

        if ( count_func != NULL )
        {
            count_func( system, start, end, d, out_bufs,
                    &cnt1, &cnt2 );
        }

        check_srealloc( &out_bufs[2 * d].out_atoms,
                &out_bufs[2 * d].out_atoms_size,
                cnt1 * MPI_Aint_add( lower_bound, extent ),
                "SendRecv::mpi_data->out_atoms" );
        check_srealloc( (void **) &out_bufs[2 * d].index,
                &out_bufs[2 * d].index_size,
                cnt1 * sizeof(int),
                "SendRecv::mpi_data->index" );

        check_srealloc( &out_bufs[2 * d + 1].out_atoms,
                &out_bufs[2 * d + 1].out_atoms_size,
                cnt2 * MPI_Aint_add( lower_bound, extent ),
                "SendRecv::mpi_data->out_atoms" );
        check_srealloc( (void **) &out_bufs[2 * d + 1].index,
                &out_bufs[2 * d + 1].index_size,
                cnt2 * sizeof(int),
                "SendRecv::mpi_data->index" );

        if ( sort_func != NULL )
        {
            sort_func( system, start, end, d, out_bufs, mpi_data );
        }
        start = end;

        /* send both messages in dimension d */
        ret = MPI_Isend( out_bufs[2 * d].out_atoms, out_bufs[2 * d].cnt, type,
                nbr1->rank, 2 * d, comm, &req1 );
        Check_MPI_Error( ret, __FILE__, __LINE__ );

        ret = MPI_Isend( out_bufs[2 * d + 1].out_atoms, out_bufs[2 * d + 1].cnt, type,
                nbr2->rank, 2 * d + 1, comm, &req2 );
        Check_MPI_Error( ret, __FILE__, __LINE__ );

        /* recv and unpack atoms from nbr1 in dimension d */
        ret = MPI_Probe( nbr1->rank, 2 * d + 1, comm, &stat1 );
        Check_MPI_Error( ret, __FILE__, __LINE__ );
        ret = MPI_Get_count( &stat1, type, &cnt1 );
        Check_MPI_Error( ret, __FILE__, __LINE__ );

        check_smalloc( &mpi_data->in1_buffer, &mpi_data->in1_buffer_size,
                cnt1 * MPI_Aint_add( lower_bound, extent ),
                "SendRecv::mpi_data->in1_buffer" );

        ret = MPI_Recv( mpi_data->in1_buffer, cnt1, type,
                nbr1->rank, 2 * d + 1, comm, &stat1 );
        Check_MPI_Error( ret, __FILE__, __LINE__ );

        if ( unpack != NULL )
        {
            unpack( system, end, mpi_data->in1_buffer, cnt1, nbr1, d );
            end += cnt1;
        }

        /* recv and unpack atoms from nbr2 in dimension d */
        ret = MPI_Probe( nbr2->rank, 2 * d, comm, &stat2 );
        Check_MPI_Error( ret, __FILE__, __LINE__ );
        ret = MPI_Get_count( &stat2, type, &cnt2 );
        Check_MPI_Error( ret, __FILE__, __LINE__ );

        check_smalloc( &mpi_data->in2_buffer, &mpi_data->in2_buffer_size,
                cnt2 * MPI_Aint_add( lower_bound, extent ),
                "SendRecv::mpi_data->in2_buffer" );

        ret = MPI_Recv( mpi_data->in2_buffer, cnt2, type,
                nbr2->rank, 2 * d, comm, &stat2 );
        Check_MPI_Error( ret, __FILE__, __LINE__ );

        if ( unpack != NULL )
        {
            unpack( system, end, mpi_data->in2_buffer, cnt2, nbr2, d );
            end += cnt2;
        }
    }

    return end;
}


void Comm_Atoms( reax_system * const system, control_params * const control,
        simulation_data * const data, storage * const workspace,
        mpi_datatypes * const mpi_data, int renbr )
{
#if defined(LOG_PERFORMANCE)
    real t_start, t_elapsed;

    t_start = 0.0;
    t_elapsed = 0.0;

    if ( system->my_rank == MASTER_NODE )
    {
        t_start = Get_Time( );
    }
#endif

    if ( renbr == TRUE )
    {
        /* transfer ownership of any previous local atoms
         * which have moved outside of this processor's local sim. box */
        system->n = SendRecv( system, mpi_data, mpi_data->mpi_atom_type,
                &Count_Transfer_Atoms, &Sort_Transfer_Atoms,
                &Unpack_Transfer_Message, TRUE );

        Bin_My_Atoms( system, workspace );
        Reorder_My_Atoms( system, workspace );

        /* exchange ghost region info with neighbors */
        system->N = SendRecv( system, mpi_data, mpi_data->boundary_atom_type,
                &Count_Boundary_Atoms, &Sort_Boundary_Atoms,
                &Unpack_Exchange_Message, TRUE );

        Bin_Boundary_Atoms( system );
    }
    else
    {
        /* provide position updates of atoms to neighboring processors */
        SendRecv( system, mpi_data, mpi_data->mpi_rvec,
                &Count_Position_Updates, &Sort_Position_Updates,
                &Unpack_Position_Updates, FALSE );
    }

#if defined(LOG_PERFORMANCE)
    if ( system->my_rank == MASTER_NODE )
    {
        t_elapsed = Get_Elapsed_Time( t_start );
        data->timing.comm += t_elapsed;
    }
#endif
}