diff --git a/PuReMD/src/bond_orders.c b/PuReMD/src/bond_orders.c
index b1e9b228511c0a7cd32e4b7fbfaa8ca44ce03fae..da07f839e2771860abcce8099aa2017796cdf759 100644
--- a/PuReMD/src/bond_orders.c
+++ b/PuReMD/src/bond_orders.c
@@ -802,7 +802,7 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
 }
 
 
-int compare_bonds( const void *p1, const void *p2 )
+static int compare_bonds( const void *p1, const void *p2 )
 {
     return ((bond_data *)p1)->nbr - ((bond_data *)p2)->nbr;
 }
diff --git a/PuReMD/src/comm_tools.c b/PuReMD/src/comm_tools.c
index 0b399825a58054fdfa6aa07f5351d4f221f45c55..9c2728f86ae1885b9abc8b4157c64cd207c5686c 100644
--- a/PuReMD/src/comm_tools.c
+++ b/PuReMD/src/comm_tools.c
@@ -53,7 +53,7 @@ void Setup_NT_Comm( reax_system* system, control_params* control,
     bndry_cut = system->bndry_cuts.ghost_cutoff;
 
     /* identify my neighbors */
-    system->num_nt_nbrs = MAX_NT_NBRS;
+    system->num_nt_nbrs = REAX_MAX_NT_NBRS;
     for ( i = 0; i < system->num_nt_nbrs; ++i )
     {
         ivec_Sum( nbr_coords, system->my_coords, r[i] ); /* actual nbr coords */
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 30a169bcdbc1402cad9c4e65d17382b2836835cf..dadb988eedf2d143d02fd26ec430aabc14b91dd5 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -55,6 +55,13 @@
 
 interaction_function Interaction_Functions[NUM_INTRS];
 
+
+static int compare_bonds( const void *p1, const void *p2 )
+{
+    return ((bond_data *)p1)->nbr - ((bond_data *)p2)->nbr;
+}
+
+
 void Dummy_Interaction( reax_system *system, control_params *control,
                         simulation_data *data, storage *workspace,
                         reax_list **lists, output_controls *out_control )
@@ -531,6 +538,11 @@ void Init_Forces( reax_system *system, control_params *control,
 
     if ( far_nbrs->format == FULL_LIST )
     {
+
+        for( i = 0; i < system->N; ++i )
+            qsort( &bonds->bond_list[Start_Index(i, bonds)],
+                    Num_Entries(i, bonds), sizeof(bond_data), compare_bonds );
+
         /* set sym_index for bonds list (far_nbrs full list) */
         for ( i = 0; i < system->N; ++i )
         {
diff --git a/PuReMD/src/torsion_angles.c b/PuReMD/src/torsion_angles.c
index 8d9949dadee37483c189e16ac79aad2c27e0b03e..4719aa33c2d79b6c07636afb69cd8dfa5857c63f 100644
--- a/PuReMD/src/torsion_angles.c
+++ b/PuReMD/src/torsion_angles.c
@@ -342,7 +342,8 @@ void Torsion_Angles( reax_system *system, control_params *control,
                                                  fbp->V2 * exp_tor1 * (1.0 - cos2omega) +
                                                  fbp->V3 * (1.0 + cos3omega) );
 
-                                    data->my_en.e_tor += e_tor = fn10 * sin_ijk * sin_jkl * CV;
+                                    e_tor = fn10 * sin_ijk * sin_jkl * CV;
+                                    data->my_en.e_tor += e_tor;
 
                                     dfn11 = (-p_tor3 * exp_tor3_DjDk +
                                              (p_tor3 * exp_tor3_DjDk - p_tor4 * exp_tor4_DjDk) *
@@ -375,9 +376,8 @@ void Torsion_Angles( reax_system *system, control_params *control,
 
                                     /* 4-body conjugation energy */
                                     fn12 = exp_cot2_ij * exp_cot2_jk * exp_cot2_kl;
-                                    data->my_en.e_con += e_con =
-                                                             fbp->p_cot1 * fn12 *
-                                                             (1.0 + (SQR(cos_omega) - 1.0) * sin_ijk * sin_jkl);
+                                    e_con = fbp->p_cot1 * fn12 * (1.0 + (SQR(cos_omega) - 1.0) * sin_ijk * sin_jkl);
+                                    data->my_en.e_con += e_con;
 
                                     Cconj = -2.0 * fn12 * fbp->p_cot1 * p_cot2 *
                                             (1.0 + (SQR(cos_omega) - 1.0) * sin_ijk * sin_jkl);