From cd9aaec0b6f772daaf52d0f2a6256a2e1948d3d5 Mon Sep 17 00:00:00 2001
From: Abdullah Alperen <alperena@msu.edu>
Date: Fri, 9 Nov 2018 11:07:25 -0500
Subject: [PATCH] PuReMD: add preprocessors for Neutral Territory

---
 PuReMD/src/allocate.c       | 19 +++++++++++---
 PuReMD/src/basic_comm.c     |  5 +++-
 PuReMD/src/basic_comm.h     |  4 +++
 PuReMD/src/box.c            |  2 ++
 PuReMD/src/comm_tools.c     | 12 ++++++++-
 PuReMD/src/comm_tools.h     |  4 +++
 PuReMD/src/forces.c         | 33 +++++++++++++++++++++++-
 PuReMD/src/grid.c           | 10 ++++++--
 PuReMD/src/init_md.c        |  8 ++++++
 PuReMD/src/io_tools.c       | 15 +++++++++--
 PuReMD/src/linear_solvers.c | 50 +++++++++++++++++++++++++++++--------
 PuReMD/src/neighbors.c      |  4 +++
 PuReMD/src/reax_types.h     | 27 +++++++++++++++-----
 13 files changed, 166 insertions(+), 27 deletions(-)

diff --git a/PuReMD/src/allocate.c b/PuReMD/src/allocate.c
index c91fae55..77b4164d 100644
--- a/PuReMD/src/allocate.c
+++ b/PuReMD/src/allocate.c
@@ -735,19 +735,24 @@ int  Allocate_MPI_Buffers( mpi_datatypes *mpi_data, int est_recv,
                                       comm );
     }
 
-    /* Note that in_nt_buffers are allocated in comm_tools.c */
-    
+#if defined(NEUTRAL_TERRITORY)
     /* Neutral Territory out buffers */
     for ( i = 0; i < REAX_MAX_NT_NBRS; ++i )
     {
+        //TODO: add an exact number
+        /* in buffers */
+        mpi_data->in_nt_buffer[i] = (void*) 
+            scalloc( my_nbrs[i].est_recv, sizeof(real), "mpibuf:in_nt_buffer", comm );
+
         mpi_buf = &( mpi_data->out_nt_buffers[i] );
         /* allocate storage for the neighbor processor i */
         mpi_buf->index = (int*)
                          scalloc( my_nbrs[i].est_send, sizeof(int), "mpibuf:nt_index", comm );
         mpi_buf->out_atoms = (void*)
-                             scalloc( my_nbrs[i].est_send, sizeof(boundary_atom), "mpibuf:nt_out_atoms",
+                             scalloc( my_nbrs[i].est_send, sizeof(real), "mpibuf:nt_out_atoms",
                                       comm );
     }
+#endif
 
     return SUCCESS;
 }
@@ -768,7 +773,7 @@ void Deallocate_MPI_Buffers( mpi_datatypes *mpi_data )
         sfree( mpi_buf->out_atoms, "mpibuf:out_atoms" );
     }
 
-
+#if defined(NEUTRAL_TERRITORY)
     for ( i = 0; i < REAX_MAX_NT_NBRS; ++i )
     {
         sfree( mpi_data->in_nt_buffer[i], "in_nt_buffer" );
@@ -777,6 +782,7 @@ void Deallocate_MPI_Buffers( mpi_datatypes *mpi_data )
         sfree( mpi_buf->index, "mpibuf:nt_index" );
         sfree( mpi_buf->out_atoms, "mpibuf:nt_out_atoms" );
     }
+#endif
 }
 
 
@@ -1028,6 +1034,8 @@ void ReAllocate( reax_system *system, control_params *control,
                 break;
             }
         }
+
+#if defined(NEUTRAL_TERRITORY)
         // also check individual outgoing Neutral Territory buffers
         mpi_flag = 0;
         for ( p = 0; p < REAX_MAX_NT_NBRS; ++p )
@@ -1041,6 +1049,7 @@ void ReAllocate( reax_system *system, control_params *control,
             }
         }
     }
+#endif
 
     if ( mpi_flag )
     {
@@ -1065,12 +1074,14 @@ void ReAllocate( reax_system *system, control_params *control,
             total_send += nbr_pr->est_send;
         }
 
+#if defined(NEUTRAL_TERRITORY)
         for ( p = 0; p < REAX_MAX_NT_NBRS; ++p )
         {
             nbr_pr   = &( system->my_nt_nbrs[p] );
             nbr_data = &( mpi_data->out_nt_buffers[p] );
             nbr_pr->est_send = MAX( nbr_data->cnt * SAFER_ZONE, MIN_SEND );
         }
+#endif
 
 #if defined(DEBUG_FOCUS)
         fprintf( stderr, "p%d: reallocating mpi_buf: recv=%d send=%d total=%dMB\n",
diff --git a/PuReMD/src/basic_comm.c b/PuReMD/src/basic_comm.c
index 35e04dc4..834034ed 100644
--- a/PuReMD/src/basic_comm.c
+++ b/PuReMD/src/basic_comm.c
@@ -129,6 +129,7 @@ void Dist( reax_system* system, mpi_datatypes *mpi_data,
 }
 
 
+#if defined(NEUTRAL_TERRITORY)
 void Dist_NT( reax_system* system, mpi_datatypes *mpi_data,
         void *buf, MPI_Datatype type, int scale, dist_packer pack )
 {
@@ -189,6 +190,7 @@ void Dist_NT( reax_system* system, mpi_datatypes *mpi_data,
     fprintf( stderr, "p%d dist: done\n", system->my_rank );
 #endif
 }
+#endif
 
 
 void real_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
@@ -299,7 +301,7 @@ void Coll( reax_system* system, mpi_datatypes *mpi_data,
 #endif
 }
 
-
+#if defined(NEUTRAL_TERRITORY)
 void Coll_NT( reax_system* system, mpi_datatypes *mpi_data,
         void *buf, MPI_Datatype type, int scale, coll_unpacker unpack )
 {
@@ -361,6 +363,7 @@ void Coll_NT( reax_system* system, mpi_datatypes *mpi_data,
     fprintf( stderr, "p%d coll: done\n", system->my_rank );
 #endif
 }
+#endif
 
 
 #endif /*PURE_REAX*/
diff --git a/PuReMD/src/basic_comm.h b/PuReMD/src/basic_comm.h
index 31afb564..8c914973 100644
--- a/PuReMD/src/basic_comm.h
+++ b/PuReMD/src/basic_comm.h
@@ -29,13 +29,17 @@ void real_packer( void*, mpi_out_data* );
 void rvec_packer( void*, mpi_out_data* );
 void rvec2_packer( void*, mpi_out_data* );
 void Dist(reax_system*, mpi_datatypes*, void*, MPI_Datatype, int, dist_packer);
+#if defined(NEUTRAL_TERRITORY)
 void Dist_NT(reax_system*, mpi_datatypes*, void*, MPI_Datatype, int, dist_packer);
+#endif
 
 void real_unpacker( void*, void*, mpi_out_data* );
 void rvec_unpacker( void*, void*, mpi_out_data* );
 void rvec2_unpacker( void*, void*, mpi_out_data* );
 void Coll( reax_system*, mpi_datatypes*, void*, MPI_Datatype, int, coll_unpacker );
+#if defined(NEUTRAL_TERRITORY)
 void Coll_NT( reax_system*, mpi_datatypes*, void*, MPI_Datatype, int, coll_unpacker );
+#endif
 
 real Parallel_Norm( real*, int, MPI_Comm );
 real Parallel_Dot( real*, real*, int, MPI_Comm );
diff --git a/PuReMD/src/box.c b/PuReMD/src/box.c
index 212d755a..d81a96db 100644
--- a/PuReMD/src/box.c
+++ b/PuReMD/src/box.c
@@ -285,7 +285,9 @@ void Setup_Environment( reax_system *system, control_params *control,
     Setup_My_Box( system, control );
     Setup_My_Ext_Box( system, control );
     Setup_Comm( system, control, mpi_data );
+#if defined(NEUTRAL_TERRITORY)
     Setup_NT_Comm( system, control, mpi_data );
+#endif
 
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "p%d coord: %d %d %d\n",
diff --git a/PuReMD/src/comm_tools.c b/PuReMD/src/comm_tools.c
index 18fb2ded..5646b9cf 100644
--- a/PuReMD/src/comm_tools.c
+++ b/PuReMD/src/comm_tools.c
@@ -25,7 +25,7 @@
 #include "tool_box.h"
 #include "vector.h"
 
-
+#if defined(NEUTRAL_TERRITORY)
 void Setup_NT_Comm( reax_system* system, control_params* control,
                  mpi_datatypes *mpi_data )
 {
@@ -95,8 +95,10 @@ void Setup_NT_Comm( reax_system* system, control_params* control,
 
     }
 }
+#endif
 
 
+#if defined(NEUTRAL_TERRITORY)
 int Sort_Neutral_Territory( reax_system *system, int dir, mpi_out_data *out_bufs, int write )
 {
     int i, cnt;
@@ -131,8 +133,10 @@ int Sort_Neutral_Territory( reax_system *system, int dir, mpi_out_data *out_bufs
 
     return cnt;
 }
+#endif
 
 
+#if defined(NEUTRAL_TERRITORY)
 void Init_Neutral_Territory( reax_system* system, mpi_datatypes *mpi_data )
 {
     int d, end, cnt;
@@ -169,7 +173,10 @@ void Init_Neutral_Territory( reax_system* system, mpi_datatypes *mpi_data )
         end += cnt;
     }
 }
+#endif
 
+
+#if defined(NEUTRAL_TERRITORY)
 void Estimate_NT_Atoms( reax_system *system, mpi_datatypes *mpi_data )
 {
     int d;
@@ -196,6 +203,7 @@ void Estimate_NT_Atoms( reax_system *system, mpi_datatypes *mpi_data )
         //Sort_Neutral_Territory( system, d, out_bufs, 1 );
     }
 }
+#endif
 
 
 void Setup_Comm( reax_system* system, control_params* control,
@@ -822,7 +830,9 @@ void Comm_Atoms( reax_system *system, control_params *control,
 
         Bin_Boundary_Atoms( system );
         
+#if defined(NEUTRAL_TERRITORY)
         Init_Neutral_Territory( system, mpi_data );
+#endif
     }
     else
     {
diff --git a/PuReMD/src/comm_tools.h b/PuReMD/src/comm_tools.h
index 787e120d..c333fa0c 100644
--- a/PuReMD/src/comm_tools.h
+++ b/PuReMD/src/comm_tools.h
@@ -25,12 +25,16 @@
 #include "reax_types.h"
 
 void Setup_Comm( reax_system*, control_params*, mpi_datatypes* );
+#if defined(NEUTRAL_TERRITORY)
 void Setup_NT_Comm( reax_system*, control_params*, mpi_datatypes* );
+#endif
 void Update_Comm( reax_system* );
 
 void Sort_Boundary_Atoms( reax_system*, int, int, int, mpi_out_data* );
 void Estimate_Boundary_Atoms( reax_system*, int, int, int, mpi_out_data* );
+#if defined(NEUTRAL_TERRITORY)
 void Estimate_NT_Atoms( reax_system*, mpi_datatypes* );
+#endif
 void Unpack_Exchange_Message( reax_system*, int, void*, int,
                               neighbor_proc*, int );
 void Unpack_Estimate_Message( reax_system*, int, void*, int,
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index c710e2ae..7b66fa32 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -346,10 +346,12 @@ void Init_Forces( reax_system *system, control_params *control,
     int jhb_top;
     int start_j, end_j;
     int btop_j;
+#if defined(NEUTRAL_TERRITORY)
     int total_cnt[6];
     int bin[6];
     int total_sum[6];
     int nt_flag;
+#endif
 
     far_nbrs = lists[FAR_NBRS];
     bonds = lists[BONDS];
@@ -371,8 +373,11 @@ void Init_Forces( reax_system *system, control_params *control,
     num_hbonds = 0;
     btop_i = 0;
     renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
+#if defined(NEUTRAL_TERRITORY)
     nt_flag = 1;
+#endif
 
+#if defined(NEUTRAL_TERRITORY)
     if( renbr )
     {
         for ( i = 0; i < 6; ++i )
@@ -405,6 +410,7 @@ void Init_Forces( reax_system *system, control_params *control,
         }
         H->NT = total_sum[5] + total_cnt[5];
     }
+#endif
 
     for ( i = 0; i < system->N; ++i )
     {
@@ -429,12 +435,14 @@ void Init_Forces( reax_system *system, control_params *control,
             local = 1;
             cutoff = control->nonb_cut;
         }
+#if defined(NEUTRAL_TERRITORY)
         else if( atom_i->nt_dir != -1 )
         {
             local = 2;
             cutoff = control->nonb_cut;
             nt_flag = 0;
         }
+#endif
         else
         {
             local = 0;
@@ -524,6 +532,7 @@ void Init_Forces( reax_system *system, control_params *control,
                             && (j < system->n || atom_i->orig_id < atom_j->orig_id))
                       || far_nbrs->format == FULL_LIST )
                     {
+#if defined(NEUTRAL_TERRITORY)
                         if( atom_j->nt_dir != -1 || j < system->n )
                         {
                             if( j < system->n )
@@ -545,6 +554,18 @@ void Init_Forces( reax_system *system, control_params *control,
                             }
                             ++Htop;
                         }
+#else
+                        H->entries[Htop].j = j;
+                        if ( control->tabulate == 0 )
+                        {
+                            H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
+                        }
+                        else
+                        {
+                            H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
+                        }
+                        ++Htop;
+#endif
                     }
 
                     /* hydrogen bond lists */
@@ -574,6 +595,7 @@ void Init_Forces( reax_system *system, control_params *control,
                         }
                     }
                 }
+#if defined(NEUTRAL_TERRITORY)
                 else if ( local == 2 )
                 {
                     /* H matrix entry */
@@ -610,6 +632,7 @@ void Init_Forces( reax_system *system, control_params *control,
                     }
 
                 }
+#endif
 
                 /* uncorrected bond orders */
                 if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
@@ -638,10 +661,12 @@ void Init_Forces( reax_system *system, control_params *control,
             if ( ihb == 1 )
                 Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
         }
+#if defined(NEUTRAL_TERRITORY)
         else if ( local == 2 )
         {
             H->end[atom_i->pos] = Htop;
         }
+#endif
     }
 
     if ( far_nbrs->format == FULL_LIST )
@@ -977,15 +1002,17 @@ void Estimate_Storages( reax_system *system, control_params *control,
             ++(*matrix_dim);
             ihb = sbp_i->p_hbond;
         }
+#if defined(NEUTRAL_TERRITORY)
         else if ( atom_i->nt_dir != -1 )
         {
             local = 2;
             cutoff = control->nonb_cut;
             // TODO: Diag entries for NT atoms?
             //++(*Htop);
-            //++(*matrix_dim);
+            ++(*matrix_dim);
             ihb = -1;
         }
+#endif
         else
         {
             local = 0;
@@ -1015,7 +1042,9 @@ void Estimate_Storages( reax_system *system, control_params *control,
                             && (j < system->n || atom_i->orig_id < atom_j->orig_id))
                       || far_nbrs->format == FULL_LIST )
                     {
+#if defined(NEUTRAL_TERRITORY)
                         if( j < system->n || (system->my_atoms[j]).nt_dir != -1 )
+#endif
                         {
                             ++(*Htop);
                         }
@@ -1039,6 +1068,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
                     }
                 }
 
+#if defined(NEUTRAL_TERRITORY)
                 else if ( local == 2 )
                 {
 
@@ -1052,6 +1082,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
                         }
                     }
                 }
+#endif
 
                 /* uncorrected bond orders */
                 if ( nbr_pj->d <= control->bond_cut )
diff --git a/PuReMD/src/grid.c b/PuReMD/src/grid.c
index 1bd1b72a..ce251430 100644
--- a/PuReMD/src/grid.c
+++ b/PuReMD/src/grid.c
@@ -35,8 +35,8 @@ void Mark_GCells( reax_system* system, grid *g, ivec procs, MPI_Comm comm )
     ivec send_span, recv_span;
     ivec str_send, end_send;
     ivec str_recv, end_recv;
+#if defined(NEUTRAL_TERRITORY)
     ivec nt_str, nt_end;
-
     ivec dir[6] = {
         {0, 0, +1}, // +z
         {0, 0, -1}, // -z
@@ -45,6 +45,7 @@ void Mark_GCells( reax_system* system, grid *g, ivec procs, MPI_Comm comm )
         {+1, 0, 0}, // +x
         {+1, -1, 0}  // +x-y
     };
+#endif
 
     /* clear all gcell type info */
     for ( x = 0; x < g->ncells[0]; x++ )
@@ -61,6 +62,7 @@ void Mark_GCells( reax_system* system, grid *g, ivec procs, MPI_Comm comm )
                 ivec_MakeZero( g->cells[x][y][z].rel_box );
             }
     
+#if defined(NEUTRAL_TERRITORY)
     /* mark NT cells */
     for ( i = 0; i < 6; ++i )
     {
@@ -93,6 +95,7 @@ void Mark_GCells( reax_system* system, grid *g, ivec procs, MPI_Comm comm )
             }
         }
     }
+#endif
 
     /* loop over neighbors */
     for ( r[0] = -1; r[0] <= 1; ++r[0])
@@ -179,8 +182,11 @@ void Find_Neighbor_GridCells( grid *g, control_params *control )
                 gc = &(g->cells[ci[0]][ci[1]][ci[2]]);
                 top = 0;
                 //fprintf( stderr, "grid1: %d %d %d:\n", ci[0], ci[1], ci[2] );
-
+#if defined(NEUTRAL_TERRITORY)
                 if ( gc->type == NATIVE || ( gc->type >= NT_NBRS && gc->type < NT_NBRS + 6 ) )
+#else
+                if ( gc->type == NATIVE )
+#endif
                     gc->cutoff = control->vlist_cut;
                 else gc->cutoff = control->bond_cut;
 
diff --git a/PuReMD/src/init_md.c b/PuReMD/src/init_md.c
index a6058ec3..b8115820 100644
--- a/PuReMD/src/init_md.c
+++ b/PuReMD/src/init_md.c
@@ -155,7 +155,9 @@ int Init_System( reax_system *system, control_params *control,
             Estimate_Boundary_Atoms, Unpack_Estimate_Message, 1 );
     system->total_cap = MAX( (int)(system->N * SAFE_ZONE), MIN_CAP );
     Bin_Boundary_Atoms( system );
+#if defined(NEUTRAL_TERRITORY)
     Estimate_NT_Atoms( system, mpi_data );
+#endif
 
     //fprintf( stderr, "p%d SEND RECV SEND!\n", system->my_rank );
     //MPI_Barrier( mpi_data->world );
@@ -469,10 +471,12 @@ int Init_MPI_Datatypes( reax_system *system, storage *workspace,
     /* init mpi buffers  */
     mpi_data->in1_buffer = NULL;
     mpi_data->in2_buffer = NULL;
+#if defined(NEUTRAL_TERRITORY)
     for ( i = 0; i < REAX_MAX_NT_NBRS; ++i )
     {
         mpi_data->in_nt_buffer[i] = NULL;
     }
+#endif
 
     /* mpi_atom - [orig_id, imprt_id, type, num_bonds, num_hbonds, name,
        x, v, f_old, s, t] */
@@ -640,7 +644,11 @@ int  Init_Lists( reax_system *system, control_params *control,
     Estimate_Storages( system, control, lists,
             &Htop, hb_top, bond_top, &num_3body, comm, &matrix_dim );
 
+#if defined(NEUTRAL_TERRITORY)
     Allocate_Matrix( &(workspace->H), matrix_dim, Htop, cm_format, comm );
+#else
+    Allocate_Matrix( &(workspace->H), system->local_cap, Htop, comm );
+#endif
     workspace->L = NULL;
     workspace->U = NULL;
     workspace->H_spar_patt = NULL;
diff --git a/PuReMD/src/io_tools.c b/PuReMD/src/io_tools.c
index b22ef6bf..1abf9753 100644
--- a/PuReMD/src/io_tools.c
+++ b/PuReMD/src/io_tools.c
@@ -851,8 +851,18 @@ void Print_HBonds( reax_system *system, reax_list **lists,
     FILE *fout;
     reax_list *hbonds = lists[HBONDS];
 
-    sprintf( fname, "%s.hbonds.%d.%d", control->sim_name, step, system->my_rank );
-    fout = sfopen( fname, "w", "Print_HBonds" );
+    snprintf( fname, MAX_STR, "hbonds.%d.%d", step, system->my_rank );
+    
+    fprintf( stdout, "p%d, Round 1, Printing H bonds, numH = %d\n", system->my_rank, system->numH );
+    fflush( stdout );
+    MPI_Barrier( MPI_COMM_WORLD );
+    
+    //fout = sfopen( fname, "w", "Print_HBonds" );
+    fout = fopen( fname, "w" );
+
+    fprintf( stdout, "p%d, Round 2, Printing H bonds, numH = %d\n", system->my_rank, system->numH );
+    fflush( stdout );
+    MPI_Barrier( MPI_COMM_WORLD );
 
     for ( i = 0; i < system->numH; ++i )
     {
@@ -867,6 +877,7 @@ void Print_HBonds( reax_system *system, reax_list **lists,
         }
     }
 
+
     sfclose( fout, "Print_HBonds" );
 }
 
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index fd068c99..36a6b15e 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -1201,20 +1201,29 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
 
     t_start = Get_Time( );
     scale = sizeof(real) / sizeof(void);
+#if defined(NEUTRAL_TERRITORY)
+    Dist_NT( system, mpi_data, x, MPI_DOUBLE, scale, real_packer );
+#else
     Dist( system, mpi_data, x, MPI_DOUBLE, scale, real_packer );
-//    Dist_NT( system, mpi_data, x, MPI_DOUBLE, scale, real_packer );
+#endif
     t_comm += Get_Timing_Info( t_start );
 
     t_start = Get_Time( );
-    //Sparse_MatVec( H, x, workspace->q, system->N );
+#if defined(NEUTRAL_TERRITORY)
     Sparse_MatVec( H, x, workspace->q, H->NT );
+#else
+    Sparse_MatVec( H, x, workspace->q, system->N );
+#endif
     t_spmv += Get_Timing_Info( t_start );
 
     if ( H->format == SYM_HALF_MATRIX )
     {
         t_start = Get_Time( );
+#if defined(NEUTRAL_TERRITORY)
+        Coll_NT( system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker );
+#else
         Coll( system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker );
-//        Coll_NT( system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker );
+#endif
         t_comm += Get_Timing_Info( t_start );
     }
 
@@ -1226,13 +1235,19 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
     if( control->cm_solver_pre_comp_type == SAI_PC )
     {
         t_start = Get_Time( );
+#if defined(NEUTRAL_TERRITORY)
+        Dist_NT( system, mpi_data, workspace->r, MPI_DOUBLE, scale, real_packer );
+#else
         Dist( system, mpi_data, workspace->r, MPI_DOUBLE, scale, real_packer );
-//        Dist_NT( system, mpi_data, workspace->r, MPI_DOUBLE, scale, real_packer );
+#endif
         t_comm += Get_Timing_Info( t_start );
         
         t_start = Get_Time( );
-        //Sparse_MatVec( workspace->H_app_inv, workspace->r, workspace->d, system->n );
+#if defined(NEUTRAL_TERRITORY)
         Sparse_MatVec( workspace->H_app_inv, workspace->r, workspace->d, H->NT );
+#else
+        Sparse_MatVec( workspace->H_app_inv, workspace->r, workspace->d, system->n );
+#endif
         t_pa += Get_Timing_Info( t_start );
     }
 
@@ -1254,20 +1269,29 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
     for ( i = 0; i < control->cm_solver_max_iters && sqrt(sig_new) / b_norm > tol; ++i )
     {
         t_start = Get_Time( );
+#if defined(NEUTRAL_TERRITORY)
+        Dist_NT( system, mpi_data, workspace->d, MPI_DOUBLE, scale, real_packer );
+#else
         Dist( system, mpi_data, workspace->d, MPI_DOUBLE, scale, real_packer );
-//        Dist_NT( system, mpi_data, workspace->d, MPI_DOUBLE, scale, real_packer );
+#endif
         t_comm += Get_Timing_Info( t_start );
 
         t_start = Get_Time( );
-        //Sparse_MatVec( H, workspace->d, workspace->q, system->N );
+#if defined(NEUTRAL_TERRITORY)
         Sparse_MatVec( H, workspace->d, workspace->q, H->NT );
+#else
+        Sparse_MatVec( H, workspace->d, workspace->q, system->N );
+#endif
         t_spmv += Get_Timing_Info( t_start );
 
         if ( H->format == SYM_HALF_MATRIX )
         {
             t_start = Get_Time( );
+#if defined(NEUTRAL_TERRITORY)
+            Coll_NT(system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker);
+#else
             Coll(system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker);
-//            Coll_NT(system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker);
+#endif
             t_comm += Get_Timing_Info( t_start );
         }
 
@@ -1285,13 +1309,19 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
         if( control->cm_solver_pre_comp_type == SAI_PC )
         {
             t_start = Get_Time( );
+#if defined(NEUTRAL_TERRITORY)
+            Dist_NT( system, mpi_data, workspace->r, MPI_DOUBLE, scale, real_packer );
+#else
             Dist( system, mpi_data, workspace->r, MPI_DOUBLE, scale, real_packer );
-//            Dist_NT( system, mpi_data, workspace->r, MPI_DOUBLE, scale, real_packer );
+#endif
             t_comm += Get_Timing_Info( t_start );
 
             t_start = Get_Time( );
-            //Sparse_MatVec( workspace->H_app_inv, workspace->r, workspace->p, system->n );
+#if defined(NEUTRAL_TERRITORY)
             Sparse_MatVec( workspace->H_app_inv, workspace->r, workspace->p, H->NT );
+#else
+            Sparse_MatVec( workspace->H_app_inv, workspace->r, workspace->p, system->n );
+#endif
             t_pa += Get_Timing_Info( t_start );
         }
 
diff --git a/PuReMD/src/neighbors.c b/PuReMD/src/neighbors.c
index c62af733..69f8c653 100644
--- a/PuReMD/src/neighbors.c
+++ b/PuReMD/src/neighbors.c
@@ -104,6 +104,7 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
                 for (l = gci->str; l < gci->end; ++l )
                 {
                     atom1 = &(system->my_atoms[l]);
+#if defined(NEUTRAL_TERRITORY)
                     if( gci->type >= NT_NBRS && gci->type < NT_NBRS + 6 )
                     {
                         atom1->nt_dir = gci->type - NT_NBRS;
@@ -112,6 +113,7 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
                     {
                         atom1->nt_dir = -1;
                     }
+#endif
                     Set_Start_Index( l, num_far, far_nbrs );
                     //fprintf( stderr, "\tatom %d\n", atom1 );
 
@@ -213,6 +215,7 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists,
                 for ( l = gci->str; l < gci->end; ++l )
                 {
                     atom1 = &(system->my_atoms[l]);
+#if defined(NEUTRAL_TERRITORY)
                     if( gci->type >= NT_NBRS && gci->type < NT_NBRS + 6 )
                     {
                         atom1->nt_dir = gci->type - NT_NBRS;
@@ -221,6 +224,7 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists,
                     {
                         atom1->nt_dir = -1;
                     }
+#endif
                     //fprintf( stderr, "\tatom %d: ", l );
                     //tmp = num_far; tested = 0;
                     itr = 0;
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index b42dcd22..9957f148 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -217,12 +217,15 @@ typedef struct
     MPI_Datatype bond_view;
     MPI_Datatype angle_line;
     MPI_Datatype angle_view;
+    mpi_out_data out_buffers[REAX_MAX_NBRS];
 
-    mpi_out_data out_buffers[REAX_MAX_NT_NBRS];
-    mpi_out_data out_nt_buffers[REAX_MAX_NT_NBRS];
     void *in1_buffer;
     void *in2_buffer;
+
+#if defined(NEUTRAL_TERRITORY)
+    mpi_out_data out_nt_buffers[REAX_MAX_NT_NBRS];
     void *in_nt_buffer[REAX_MAX_NT_NBRS];
+#endif
 } mpi_datatypes;
 
 
@@ -435,8 +438,10 @@ typedef struct
     int num_bonds;
     int num_hbonds;
     int renumber;
+#if defined(NEUTRAL_TERRITORY)
     int nt_dir;
     int pos;
+#endif
 } reax_atom;
 
 
@@ -544,16 +549,20 @@ typedef struct
     int              n, N, bigN, numH;
     int              local_cap, total_cap, gcell_cap, Hcap;
     int              est_recv, est_trans, max_recved;
-    int              wsize, my_rank, num_nbrs, num_nt_nbrs;
+    int              wsize, my_rank, num_nbrs;
     ivec             my_coords;
     neighbor_proc    my_nbrs[REAX_MAX_NBRS];
-    neighbor_proc    my_nt_nbrs[REAX_MAX_NT_NBRS];
     int             *global_offset;
     simulation_box   big_box, my_box, my_ext_box;
     grid             my_grid;
     boundary_cutoff  bndry_cuts;
 
     reax_atom       *my_atoms;
+
+#if defined(NEUTRAL_TERRITORY)
+    int              num_nt_nbrs;
+    neighbor_proc    my_nt_nbrs[REAX_MAX_NT_NBRS];
+#endif
 } reax_system;
 
 
@@ -887,6 +896,7 @@ typedef struct
 } far_neighbor_data;
 
 
+#if defined(NEUTRAL_TERRITORY)
 typedef struct
 {
     int nbr;
@@ -894,6 +904,7 @@ typedef struct
     real d;
     rvec dvec;
 } nt_neighbor_data;
+#endif
 
 
 typedef struct
@@ -966,7 +977,10 @@ typedef struct
 {
     /* matrix storage format */
     int format;
-    int cap, n, m, NT;
+    int cap, n, m;
+#if defined(NEUTRAL_TERRITORY)
+    int NT;
+#endif
     int *start, *end;
     sparse_matrix_entry *entries;
 } sparse_matrix;
@@ -1080,7 +1094,6 @@ typedef union
     //dbond_data         *dbo_list;
     //dDelta_data        *dDelta_list;
     //far_neighbor_data  *far_nbr_list;
-    //nt_neighbor_data   *nt_nbr_list;
     //hbond_data         *hbond_list;
 } list_type;
 
@@ -1107,7 +1120,9 @@ typedef struct
     dbond_data *dbo_list;
     dDelta_data *dDelta_list;
     far_neighbor_data *far_nbr_list;
+#if defined(NEUTRAL_TERRITORY)
     nt_neighbor_data *nt_nbr_list;
+#endif
     hbond_data *hbond_list;
 } reax_list;
 
-- 
GitLab