From 3f7a446e0b5ae3f1885684126a6e1528621d4736 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Wed, 9 Sep 2020 00:24:16 -0400
Subject: [PATCH] Revert "PG-PuReMD: remove unused members for structs involved
 in MPI comms for atom exchanges."

This reverts commit 6a1385e9027e1d0dd59fd917b8f3af4f01c2146c.
---
 PG-PuReMD/src/comm_tools.c        | 13 ++++++
 PG-PuReMD/src/cuda/cuda_forces.cu | 12 +++++-
 PG-PuReMD/src/forces.c            | 11 +++++
 PG-PuReMD/src/geo_tools.c         |  6 +++
 PG-PuReMD/src/grid.c              |  1 +
 PG-PuReMD/src/init_md.c           | 68 +++++++++++++++++++------------
 PG-PuReMD/src/io_tools.c          | 12 +++---
 PG-PuReMD/src/reax_types.h        | 22 +++++++++-
 8 files changed, 111 insertions(+), 34 deletions(-)

diff --git a/PG-PuReMD/src/comm_tools.c b/PG-PuReMD/src/comm_tools.c
index a5dd4052..3371568b 100644
--- a/PG-PuReMD/src/comm_tools.c
+++ b/PG-PuReMD/src/comm_tools.c
@@ -297,7 +297,10 @@ void Update_Comm( reax_system * const system )
 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 );
@@ -311,7 +314,10 @@ static void Pack_MPI_Atom( mpi_atom * const matm, const reax_atom * const ratm,
 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 );
@@ -432,7 +438,10 @@ static void Pack_Boundary_Atom( boundary_atom * const matm,
         reax_atom const * 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 );
 }
 
@@ -441,7 +450,11 @@ 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 );
 }
 
diff --git a/PG-PuReMD/src/cuda/cuda_forces.cu b/PG-PuReMD/src/cuda/cuda_forces.cu
index 0a69e733..43ea11bf 100644
--- a/PG-PuReMD/src/cuda/cuda_forces.cu
+++ b/PG-PuReMD/src/cuda/cuda_forces.cu
@@ -151,7 +151,7 @@ CUDA_GLOBAL void k_init_hindex( reax_atom *my_atoms, int N )
 CUDA_GLOBAL void k_init_hbond_indices( reax_atom * atoms, single_body_parameters *sbp,
         int *hbonds, int *max_hbonds, int *indices, int *end_indices, int N )
 {
-    int i, hindex;
+    int i, hindex, my_hbonds;
 
     i = blockIdx.x * blockDim.x  + threadIdx.x;
 
@@ -165,14 +165,17 @@ CUDA_GLOBAL void k_init_hbond_indices( reax_atom * atoms, single_body_parameters
     if ( sbp[ atoms[i].type ].p_hbond == H_ATOM
             || sbp[ atoms[i].type ].p_hbond == H_BONDING_ATOM )
     {
+        my_hbonds = hbonds[i];
         indices[hindex] = max_hbonds[i];
         end_indices[hindex] = indices[hindex] + hbonds[i];
     }
     else
     {
+        my_hbonds = 0;
         indices[hindex] = 0;
         end_indices[hindex] = 0;
     }
+    atoms[i].num_hbonds = my_hbonds;
 }
 
 
@@ -742,6 +745,11 @@ CUDA_GLOBAL void k_init_bonds( reax_atom *my_atoms, single_body_parameters *sbp,
     num_bonds = btop_i - Start_Index( i, &bond_list );
     num_hbonds = ihb_top - Start_Index( atom_i->Hindex, &hbond_list );
 
+    /* copy (h)bond info to atom structure
+     * (needed for atom ownership transfer via MPI) */
+    my_atoms[i].num_bonds = num_bonds;
+    my_atoms[i].num_hbonds = num_hbonds;
+
     /* reallocation checks */
     if ( num_bonds > max_bonds[i] )
     {
@@ -1466,7 +1474,7 @@ int Cuda_Init_Forces( reax_system *system, control_params *control,
     int renbr, blocks, ret, realloc_bonds, realloc_hbonds, realloc_cm;
     static int dist_done = FALSE, cm_done = FALSE, bonds_done = FALSE;
 #if defined(LOG_PERFORMANCE)
-    real time;
+    double time;
     
     time = Get_Time( );
 #endif
diff --git a/PG-PuReMD/src/forces.c b/PG-PuReMD/src/forces.c
index 8b79e576..29e220c7 100644
--- a/PG-PuReMD/src/forces.c
+++ b/PG-PuReMD/src/forces.c
@@ -2093,6 +2093,17 @@ void Estimate_Storages( reax_system * const system, control_params * const contr
     system->total_cm_entries = MAX( system->total_cm_entries, MIN_CAP * MIN_CM_ENTRIES );
     system->total_thbodies = MAX( system->total_thbodies * SAFE_ZONE, MIN_3BODIES );
 
+    /* duplicate info in atom structs in case of
+     * ownership transfer across processor boundaries */
+    for ( i = 0; i < system->n; ++i )
+    {
+        system->my_atoms[i].num_bonds = system->bonds[i];
+    }
+    for ( i = 0; i < system->N; ++i )
+    {
+        system->my_atoms[i].num_hbonds = system->hbonds[i];
+    }
+
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "[INFO] p%d @ estimate storages: total_cm_entries = %d, total_thbodies = %d\n",
             system->my_rank, system->total_cm_entries, system->total_thbodies );
diff --git a/PG-PuReMD/src/geo_tools.c b/PG-PuReMD/src/geo_tools.c
index a2340304..3b3b0c1a 100644
--- a/PG-PuReMD/src/geo_tools.c
+++ b/PG-PuReMD/src/geo_tools.c
@@ -133,7 +133,10 @@ void Read_Geo_File( const char * const geo_file, reax_system * const system,
             rvec_MakeZero( atom->s );
             rvec_MakeZero( atom->t );
 
+            atom->imprt_id = -1;
             atom->Hindex = -1;
+            atom->num_bonds = 0;
+            atom->num_hbonds = 0;
 
             top++;
         }
@@ -401,7 +404,10 @@ void Read_PDB_File( const char * const pdb_file, reax_system * const system,
                 rvec_MakeZero( atom->s );
                 rvec_MakeZero( atom->t );
 
+                atom->imprt_id = -1;
                 atom->Hindex = -1;
+                atom->num_bonds = 0;
+                atom->num_hbonds = 0;
 
                 top++;
             }
diff --git a/PG-PuReMD/src/grid.c b/PG-PuReMD/src/grid.c
index d91470fd..87df81d1 100644
--- a/PG-PuReMD/src/grid.c
+++ b/PG-PuReMD/src/grid.c
@@ -672,6 +672,7 @@ void Reorder_My_Atoms( reax_system * const system, storage * const workspace )
             old_id = gc->atoms[l];
             old_atom = &system->my_atoms[old_id];
             memcpy( &new_atoms[top], old_atom, sizeof(reax_atom) );
+            new_atoms[top].imprt_id = -1;
             ++top;
         }
         g->end[ index_grid_3d(x, y, z, g) ] = top;
diff --git a/PG-PuReMD/src/init_md.c b/PG-PuReMD/src/init_md.c
index 4ec939ab..0aa321f4 100644
--- a/PG-PuReMD/src/init_md.c
+++ b/PG-PuReMD/src/init_md.c
@@ -443,37 +443,46 @@ void Init_MPI_Datatypes( reax_system * const system, storage * const workspace,
     /* mpi_atom */
     block[0] = 1;
     block[1] = 1;
-    block[2] = sizeof(m_sample[0].name) / sizeof(char);
-    block[3] = 3;
-    block[4] = 3;
-    block[5] = 3;
-    block[6] = 4;
-    block[7] = 4;
+    block[2] = 1;
+    block[3] = 1;
+    block[4] = 1;
+    block[5] = sizeof(m_sample[0].name) / sizeof(char);
+    block[6] = 3;
+    block[7] = 3;
+    block[8] = 3;
+    block[9] = 4;
+    block[10] = 4;
 
     MPI_Get_address( &m_sample[0], &base );
     MPI_Get_address( &m_sample[0].orig_id, &disp[0] );
-    MPI_Get_address( &m_sample[0].type, &disp[1] );
-    MPI_Get_address( &m_sample[0].name, &disp[2] );
-    MPI_Get_address( &m_sample[0].x, &disp[3] );
-    MPI_Get_address( &m_sample[0].v, &disp[4] );
-    MPI_Get_address( &m_sample[0].f_old, &disp[5] );
-    MPI_Get_address( &m_sample[0].s, &disp[6] );
-    MPI_Get_address( &m_sample[0].t, &disp[7] );
-    for ( i = 0; i < 8; ++i )
+    MPI_Get_address( &m_sample[0].imprt_id, &disp[1] );
+    MPI_Get_address( &m_sample[0].type, &disp[2] );
+    MPI_Get_address( &m_sample[0].num_bonds, &disp[3] );
+    MPI_Get_address( &m_sample[0].num_hbonds, &disp[4] );
+    MPI_Get_address( &m_sample[0].name, &disp[5] );
+    MPI_Get_address( &m_sample[0].x, &disp[6] );
+    MPI_Get_address( &m_sample[0].v, &disp[7] );
+    MPI_Get_address( &m_sample[0].f_old, &disp[8] );
+    MPI_Get_address( &m_sample[0].s, &disp[9] );
+    MPI_Get_address( &m_sample[0].t, &disp[10] );
+    for ( i = 0; i < 11; ++i )
     {
         disp[i] = MPI_Aint_diff( disp[i], base );
     }
 
     type[0] = MPI_INT;
     type[1] = MPI_INT;
-    type[2] = MPI_CHAR;
-    type[3] = MPI_DOUBLE;
-    type[4] = MPI_DOUBLE;
-    type[5] = MPI_DOUBLE;
+    type[2] = MPI_INT;
+    type[3] = MPI_INT;
+    type[4] = MPI_INT;
+    type[5] = MPI_CHAR;
     type[6] = MPI_DOUBLE;
     type[7] = MPI_DOUBLE;
+    type[8] = MPI_DOUBLE;
+    type[9] = MPI_DOUBLE;
+    type[10] = MPI_DOUBLE;
 
-    ret = MPI_Type_create_struct( 8, block, disp, type, &temp_type );
+    ret = MPI_Type_create_struct( 11, block, disp, type, &temp_type );
     Check_MPI_Error( ret, __FILE__, __LINE__ );
     /* account for struct padding after members */
     ret = MPI_Type_create_resized( temp_type, 0, sizeof(mpi_atom),
@@ -487,22 +496,31 @@ void Init_MPI_Datatypes( reax_system * const system, storage * const workspace,
     /* boundary_atom */
     block[0] = 1;
     block[1] = 1;
-    block[2] = 3;
+    block[2] = 1;
+    block[3] = 1;
+    block[4] = 1;
+    block[5] = 3;
 
     MPI_Get_address( &b_sample[0], &base );
     MPI_Get_address( &b_sample[0].orig_id, &disp[0] );
-    MPI_Get_address( &b_sample[0].type, &disp[1] );
-    MPI_Get_address( &b_sample[0].x, &disp[2] );
-    for ( i = 0; i < 3; ++i )
+    MPI_Get_address( &b_sample[0].imprt_id, &disp[1] );
+    MPI_Get_address( &b_sample[0].type, &disp[2] );
+    MPI_Get_address( &b_sample[0].num_bonds, &disp[3] );
+    MPI_Get_address( &b_sample[0].num_hbonds, &disp[4] );
+    MPI_Get_address( &b_sample[0].x, &disp[5] );
+    for ( i = 0; i < 6; ++i )
     {
         disp[i] = MPI_Aint_diff( disp[i], base );
     }
 
     type[0] = MPI_INT;
     type[1] = MPI_INT;
-    type[2] = MPI_DOUBLE;
+    type[2] = MPI_INT;
+    type[3] = MPI_INT;
+    type[4] = MPI_INT;
+    type[5] = MPI_DOUBLE;
 
-    ret = MPI_Type_create_struct( 3, block, disp, type, &temp_type );
+    ret = MPI_Type_create_struct( 6, block, disp, type, &temp_type );
     Check_MPI_Error( ret, __FILE__, __LINE__ );
     /* account for struct padding after members */
     ret = MPI_Type_create_resized( temp_type, 0, sizeof(boundary_atom),
diff --git a/PG-PuReMD/src/io_tools.c b/PG-PuReMD/src/io_tools.c
index 45e46f51..d97655b7 100644
--- a/PG-PuReMD/src/io_tools.c
+++ b/PG-PuReMD/src/io_tools.c
@@ -742,9 +742,9 @@ void Print_My_Ext_Atoms( reax_system *system )
 
     for ( i = 0; i < system->N; ++i )
     {
-        fprintf( fh, "p%-2d %-5d %2d %24.15e%24.15e%24.15e\n",
+        fprintf( fh, "p%-2d %-5d imprt%-5d %2d %24.15e%24.15e%24.15e\n",
                  system->my_rank, system->my_atoms[i].orig_id,
-                 system->my_atoms[i].type,
+                 system->my_atoms[i].imprt_id, system->my_atoms[i].type,
                  system->my_atoms[i].x[0],
                  system->my_atoms[i].x[1],
                  system->my_atoms[i].x[2] );
@@ -851,12 +851,12 @@ void Print_Symmetric_Sparse( reax_system *system, sparse_matrix *A, char *fname
             aj = &system->my_atoms[A->j[j]];
 
             fprintf( f, "%d %d %.15e\n",
-                    ai->orig_id, aj->orig_id, A->val[j] );
+                    ai->renumber, aj->renumber, A->val[j] );
 
-            if ( A->j[j] < system->n && ai->orig_id != aj->orig_id )
+            if ( A->j[j] < system->n && ai->renumber != aj->renumber )
             {
                 fprintf( f, "%d %d %.15e\n",
-                         aj->orig_id, ai->orig_id, A->val[j] );
+                         aj->renumber, ai->renumber, A->val[j] );
             }
         }
     }
@@ -884,7 +884,7 @@ void Print_Linear_System( reax_system *system, control_params *control,
     {
         ai = &system->my_atoms[i];
         fprintf( out, "%6d%2d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
-                 ai->orig_id, ai->type, ai->x[0], ai->x[1], ai->x[2],
+                 ai->renumber, ai->type, ai->x[0], ai->x[1], ai->x[2],
                  workspace->s[i], workspace->b_s[i],
                  workspace->t[i], workspace->b_t[i] );
     }
diff --git a/PG-PuReMD/src/reax_types.h b/PG-PuReMD/src/reax_types.h
index d372c158..194970ac 100644
--- a/PG-PuReMD/src/reax_types.h
+++ b/PG-PuReMD/src/reax_types.h
@@ -633,10 +633,16 @@ struct mpi_atom
 {
     /* atom serial number as given in the geo file */
     int orig_id;
+    /* local atom ID on neighbor processor ??? */
+    int imprt_id;
     /* non-negative integer used to indicate atom type,
      * as identified by short element string in force field file (single
      * body parameters section) */
     int type;
+    /* num. bonds associated with atom */
+    int num_bonds;
+    /* num. hydrogren bonds associated with atom */
+    int num_hbonds;
     /* atom name as given in the geo file */
     char name[MAX_ATOM_NAME_LEN];
     /* atomic position, 3D */
@@ -659,10 +665,16 @@ struct boundary_atom
 {
     /* atom serial number as given in the geo file */
     int orig_id;
+    /* local atom ID on neighbor processor ??? */
+    int imprt_id;
     /* non-negative integer used to indicate atom type,
      * as identified by short element string in force field file (single
      * body parameters section) */
     int type;
+    /* num. bonds associated with atom */
+    int num_bonds;
+    /* num. hydrogren bonds associated with atom */
+    int num_hbonds;
     /* atomic position, 3D */
     rvec x;
 };
@@ -1046,11 +1058,13 @@ struct reax_interaction
 };
 
 
-/* struct containing atom-specific information */
+/**/
 struct reax_atom
 {
     /* atom serial number as given in the geo file */
     int orig_id;
+    /* local atom ID on neighbor processor ??? */
+    int imprt_id;
     /* non-negative integer used to indicate atom type,
      * as identified by short element string in force field file (single
      * body parameters section) */
@@ -1076,6 +1090,12 @@ struct reax_atom
     /* unique non-negative integer index of atom if it is a hydrogen atom,
      * -1 otherwise */
     int Hindex;
+    /* num. bonds associated with atom */
+    int num_bonds;
+    /* num. hydrogren bonds associated with atom */
+    int num_hbonds;
+    /* ??? */
+    int renumber;
 #if defined(NEUTRAL_TERRITORY)
     /**/
     int nt_dir;
-- 
GitLab