From d6406285e9e97a3ca152b5bbde45b6f5edcbe72d Mon Sep 17 00:00:00 2001
From: Abdullah Alperen <alperena@msu.edu>
Date: Thu, 10 Jan 2019 10:50:02 -0500
Subject: [PATCH] PuReMD: fix bond_list

---
 PuReMD/src/forces.c | 44 +++++++++++++++++++++++++-------------------
 1 file changed, 25 insertions(+), 19 deletions(-)

diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 4588515d..e57adbb2 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -1146,23 +1146,23 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
             j = far_nbrs->far_nbr_list.nbr[pj];
             atom_j = &system->my_atoms[j];
 
-            if ( far_nbrs->far_nbr_list.d[pj] <= control->nonb_cut )
-//            if ( i <= j && far_nbrs->far_nbr_list.d[pj] <= control->nonb_cut )
+//            if ( far_nbrs->far_nbr_list.d[pj] <= control->nonb_cut )
+            if ( i <= j && far_nbrs->far_nbr_list.d[pj] <= control->nonb_cut )
             {
                 type_j = atom_j->type;
                 sbp_j = &system->reax_param.sbp[type_j];
                 twbp = &system->reax_param.tbp[type_i][type_j];
 
-                if ( BOp_redundant( workspace, bonds, control->bo_cut,
-                            i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
-                            &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
-                            &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
-                            sbp_i, sbp_j, twbp ) )
-//                if ( BOp( workspace, bonds, control->bo_cut,
+//                if ( BOp_redundant( workspace, bonds, control->bo_cut,
 //                            i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
 //                            &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
 //                            &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
 //                            sbp_i, sbp_j, twbp ) )
+                if ( BOp( workspace, bonds, control->bo_cut,
+                            i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+                            &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+                            &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
+                            sbp_i, sbp_j, twbp ) )
                 {
                     num_bonds += 2;
                     ++btop_i;
@@ -1171,7 +1171,7 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
                      * to search for it's bonded neighbors later */
                     if ( workspace->bond_mark[j] == 1000 )
                     {
-                        workspace->bond_mark[j] = 1;
+                        workspace->bond_mark[j] = 101;
                         q[ push++ ] = j;
                     }
                 }
@@ -1185,6 +1185,7 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
     for ( k = 0; k < push; ++k )
     {
         i = q[k];
+        workspace->bond_mark[i] -= 100;
         atom_i = &system->my_atoms[i];
         type_i = atom_i->type;
         btop_i = End_Index( i, bonds );
@@ -1201,35 +1202,40 @@ static void Init_Bond_Full( reax_system *system, control_params *control,
             {
                 continue;
             }
+/*            if ( i>j && workspace->bond_mark[i] > workspace->bond_mark[j] 
+                    && workspace->bond_mark[j] == 1000 )
+            {
+                continue;
+            }*/
 
             atom_j = &system->my_atoms[j];
 
-            if ( far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
-//            if ( i <= j && far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
+//            if ( far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
+            if (  workspace->bond_mark[j] > 100 && far_nbrs->far_nbr_list.d[pj] <= control->bond_cut )
             {
                 type_j = atom_j->type;
                 sbp_j = &system->reax_param.sbp[type_j];
                 twbp = &system->reax_param.tbp[type_i][type_j];
 
-                if ( BOp_redundant( workspace, bonds, control->bo_cut,
-                            i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
-                            &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
-                            &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
-                            sbp_i, sbp_j, twbp ) )
-//                if ( BOp( workspace, bonds, control->bo_cut,
+//                if ( BOp_redundant( workspace, bonds, control->bo_cut,
 //                            i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
 //                            &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
 //                            &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
 //                            sbp_i, sbp_j, twbp ) )
+                if ( BOp( workspace, bonds, control->bo_cut,
+                            i, btop_i, far_nbrs->far_nbr_list.nbr[pj],
+                            &far_nbrs->far_nbr_list.rel_box[pj], far_nbrs->far_nbr_list.d[pj],
+                            &far_nbrs->far_nbr_list.dvec[pj], far_nbrs->format,
+                            sbp_i, sbp_j, twbp ) )
                 {
                     num_bonds += 2;
                     ++btop_i;
 
                     if ( workspace->bond_mark[j] == 1000 )
                     {
-                        workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
+                        workspace->bond_mark[j] = workspace->bond_mark[i] + 100;
 
-                        if ( workspace->bond_mark[j] < 4 )
+                        if ( workspace->bond_mark[i] < 3 )
                         {
                             q[ push++ ] = j;
                         }
-- 
GitLab