PuReMD-old: clean-up and comment initialization optimizations. WIP: remove...

PuReMD-old: clean-up and comment initialization optimizations. WIP: remove redundant computation of bonds list (not working for full far nbrs list). Improve hbond list computation for full far nbr list case.

PuReMD/src/allocate.c

@@ -145,7 +145,6 @@ void DeAllocate_Workspace( control_params *control, storage *workspace )
     sfree( workspace->Clp, "Clp" );
     sfree( workspace->vlpex, "vlpex" );
     sfree( workspace->bond_mark, "bond_mark" );
-    sfree( workspace->done_after, "done_after" );
 
     /* CM storage */
     sfree( workspace->Hdia_inv, "Hdia_inv" );
@@ -304,7 +303,6 @@ int Allocate_Workspace( reax_system *system, control_params *control,
     workspace->Clp = smalloc( total_real, "Clp", comm );
     workspace->vlpex = smalloc( total_real, "vlpex", comm );
     workspace->bond_mark = scalloc( total_cap, sizeof(int), "bond_mark", comm );
-    workspace->done_after = scalloc( total_cap, sizeof(int), "done_after", comm );
 
     /* CM storage */
     workspace->Hdia_inv = scalloc( total_cap, sizeof(real), "Hdia_inv", comm );

PuReMD/src/bond_orders.c

@@ -662,6 +662,10 @@ void Add_dBond_to_Forces( int i, int pj,
 }
 
 
+/* Compute the bond order term between atoms i and j,
+ * and if this term exceeds the cutoff bo_cut, then adds
+ * BOTH atoms the bonds list (i.e., compute term once
+ * and copy to avoid redundant computation) */
 int BOp( storage *workspace, reax_list *bonds, real bo_cut,
          int i, int btop_i, int j, ivec *rel_box, real d, rvec *dvec,
          int far_nbr_list_format, single_body_parameters *sbp_i,
@@ -683,21 +687,169 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
         C12 = twbp->p_bo1 * pow( d / twbp->r_s, twbp->p_bo2 );
         BO_s = (1.0 + bo_cut) * exp( C12 );
     }
-    else BO_s = C12 = 0.0;
+    else
+    {
+        C12 = 0.0;
+        BO_s = 0.0;
+    }
 
     if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
     {
         C34 = twbp->p_bo3 * pow( d / twbp->r_p, twbp->p_bo4 );
         BO_pi = exp( C34 );
     }
-    else BO_pi = C34 = 0.0;
+    else
+    {
+        C34 = 0.0;
+        BO_pi = 0.0;
+    }
 
     if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
     {
         C56 = twbp->p_bo5 * pow( d / twbp->r_pp, twbp->p_bo6 );
         BO_pi2 = exp( C56 );
     }
-    else BO_pi2 = C56 = 0.0;
+    else
+    {
+        C56 = 0.0;
+        BO_pi2 = 0.0;
+    }
+
+    /* Initially BO values are the uncorrected ones, page 1 */
+    BO = BO_s + BO_pi + BO_pi2;
+
+    if ( BO >= bo_cut )
+    {
+        /****** bonds i-j and j-i ******/
+        ibond = &bonds->bond_list[btop_i];
+        btop_j = End_Index( j, bonds );
+        jbond = &bonds->bond_list[btop_j];
+
+        ibond->nbr = j;
+        ibond->d = d;
+        rvec_Copy( ibond->dvec, *dvec );
+        ivec_Copy( ibond->rel_box, *rel_box );
+        ibond->dbond_index = btop_i;
+        ibond->sym_index = btop_j;
+        jbond->nbr = i;
+        jbond->d = d;
+        rvec_Scale( jbond->dvec, -1.0, *dvec );
+        ivec_Scale( jbond->rel_box, -1.0, *rel_box );
+        jbond->dbond_index = btop_i;
+        jbond->sym_index = btop_i;
+
+        Set_End_Index( j, btop_j + 1, bonds );
+        
+        bo_ij = &ibond->bo_data;
+        bo_ij->BO = BO;
+        bo_ij->BO_s = BO_s;
+        bo_ij->BO_pi = BO_pi;
+        bo_ij->BO_pi2 = BO_pi2;
+        bo_ji = &jbond->bo_data;
+        bo_ji->BO = BO;
+        bo_ji->BO_s = BO_s;
+        bo_ji->BO_pi = BO_pi;
+        bo_ji->BO_pi2 = BO_pi2;
+
+        /* Bond Order page2-3, derivative of total bond order prime */
+        Cln_BOp_s = twbp->p_bo2 * C12 / r2;
+        Cln_BOp_pi = twbp->p_bo4 * C34 / r2;
+        Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2;
+
+        /* Only dln_BOp_xx wrt. dr_i is stored here, note that
+         * dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
+        rvec_Scale( bo_ij->dln_BOp_s, -1.0 * bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
+        rvec_Scale( bo_ij->dln_BOp_pi, -1.0 * bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
+        rvec_Scale( bo_ij->dln_BOp_pi2, -1.0 * bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
+        rvec_Scale( bo_ji->dln_BOp_s, -1.0, bo_ij->dln_BOp_s );
+        rvec_Scale( bo_ji->dln_BOp_pi, -1.0, bo_ij->dln_BOp_pi );
+        rvec_Scale( bo_ji->dln_BOp_pi2, -1.0, bo_ij->dln_BOp_pi2 );
+
+        /* Only dBOp wrt. dr_i is stored here, note that
+         * dBOp/dr_i = -dBOp/dr_j and all others are 0 */
+        rvec_Scale( bo_ij->dBOp, -1.0 * (bo_ij->BO_s * Cln_BOp_s 
+                    + bo_ij->BO_pi * Cln_BOp_pi 
+                    + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
+        rvec_Scale( bo_ji->dBOp, -1.0, bo_ij->dBOp );
+
+        rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
+        rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );
+
+        bo_ij->BO_s -= bo_cut;
+        bo_ij->BO -= bo_cut;
+        workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
+        bo_ij->Cdbo = 0.0;
+        bo_ij->Cdbopi = 0.0;
+        bo_ij->Cdbopi2 = 0.0;
+        bo_ji->BO_s -= bo_cut;
+        bo_ji->BO -= bo_cut;
+        workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
+        bo_ji->Cdbo = 0.0;
+        bo_ji->Cdbopi = 0.0;
+        bo_ji->Cdbopi2 = 0.0;
+
+        return 1;
+    }
+
+    return 0;
+}
+
+
+/* Compute the bond order term between atoms i and j,
+ * and if this term exceeds the cutoff bo_cut, then adds
+ * to the bond list according to the following convention:
+ *   * if the far neighbor list is store in half format,
+ *      add BOTH atoms to each other's portion of the bond list
+ *   * if the far neighbor list is store in full format,
+ *      add atom i to atom j's bonds list ONLY */
+int BOp_redundant( storage *workspace, reax_list *bonds, real bo_cut,
+         int i, int btop_i, int j, ivec *rel_box, real d, rvec *dvec,
+         int far_nbr_list_format, single_body_parameters *sbp_i,
+         single_body_parameters *sbp_j, two_body_parameters *twbp )
+{
+    real r2, C12, C34, C56;
+    real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
+    real BO, BO_s, BO_pi, BO_pi2;
+    bond_data *ibond;
+    bond_order_data *bo_ij;
+    int btop_j;
+    bond_data *jbond;
+    bond_order_data *bo_ji;
+
+    r2 = SQR(d);
+
+    if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 )
+    {
+        C12 = twbp->p_bo1 * pow( d / twbp->r_s, twbp->p_bo2 );
+        BO_s = (1.0 + bo_cut) * exp( C12 );
+    }
+    else
+    {
+        C12 = 0.0;
+        BO_s = 0.0;
+    }
+
+    if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
+    {
+        C34 = twbp->p_bo3 * pow( d / twbp->r_p, twbp->p_bo4 );
+        BO_pi = exp( C34 );
+    }
+    else
+    {
+        C34 = 0.0;
+        BO_pi = 0.0;
+    }
+
+    if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
+    {
+        C56 = twbp->p_bo5 * pow( d / twbp->r_pp, twbp->p_bo6 );
+        BO_pi2 = exp( C56 );
+    }
+    else
+    {
+        C56 = 0.0;
+        BO_pi2 = 0.0;
+    }
 
     /* Initially BO values are the uncorrected ones, page 1 */
     BO = BO_s + BO_pi + BO_pi2;
@@ -730,14 +882,14 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
             Set_End_Index( j, btop_j + 1, bonds );
         }
         
-        bo_ij = &( ibond->bo_data );
+        bo_ij = &ibond->bo_data;
         bo_ij->BO = BO;
         bo_ij->BO_s = BO_s;
         bo_ij->BO_pi = BO_pi;
         bo_ij->BO_pi2 = BO_pi2;
         if ( far_nbr_list_format == HALF_LIST )
         {
-            bo_ji = &( jbond->bo_data );
+            bo_ji = &jbond->bo_data;
             bo_ji->BO = BO;
             bo_ji->BO_s = BO_s;
             bo_ji->BO_pi = BO_pi;
@@ -750,7 +902,7 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
         Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2;
 
         /* Only dln_BOp_xx wrt. dr_i is stored here, note that
-           dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
+         * dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
         rvec_Scale( bo_ij->dln_BOp_s, -1.0 * bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
         rvec_Scale( bo_ij->dln_BOp_pi, -1.0 * bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
         rvec_Scale( bo_ij->dln_BOp_pi2, -1.0 * bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
@@ -762,7 +914,7 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
         }
 
         /* Only dBOp wrt. dr_i is stored here, note that
-           dBOp/dr_i = -dBOp/dr_j and all others are 0 */
+         * dBOp/dr_i = -dBOp/dr_j and all others are 0 */
         rvec_Scale( bo_ij->dBOp, -1.0 * (bo_ij->BO_s * Cln_BOp_s 
                     + bo_ij->BO_pi * Cln_BOp_pi 
                     + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );

PuReMD/src/bond_orders.h

@@ -24,6 +24,7 @@
 
 #include "reax_types.h"
 
+
 typedef struct
 {
     real C1dbo, C2dbo, C3dbo;
@@ -32,29 +33,42 @@ typedef struct
     real C1dDelta, C2dDelta, C3dDelta;
 } dbond_coefficients;
 
+
 #ifdef TEST_FORCES
 void Get_dBO( reax_system*, reax_list**, int, int, real, rvec* );
+
 void Get_dBOpinpi2( reax_system*, reax_list**,
-                    int, int, real, real, rvec*, rvec* );
+        int, int, real, real, rvec*, rvec* );
 
 void Add_dBO( reax_system*, reax_list**, int, int, real, rvec* );
+
 void Add_dBOpinpi2( reax_system*, reax_list**,
-                    int, int, real, real, rvec*, rvec* );
+        int, int, real, real, rvec*, rvec* );
 
 void Add_dBO_to_Forces( reax_system*, reax_list**, int, int, real );
+
 void Add_dBOpinpi2_to_Forces( reax_system*, reax_list**,
-                              int, int, real, real );
+        int, int, real, real );
 
 void Add_dDelta( reax_system*, reax_list**, int, real, rvec* );
+
 void Add_dDelta_to_Forces( reax_system *, reax_list**, int, real );
 #endif
 
 void Add_dBond_to_Forces( int, int, storage*, reax_list** );
+
 void Add_dBond_to_Forces_NPT( int, int, simulation_data*,
-                              storage*, reax_list** );
+        storage*, reax_list** );
+
 int BOp( storage*, reax_list*, real, int, int, int, ivec*, real, rvec*,
         int, single_body_parameters*, single_body_parameters*,
         two_body_parameters* );
+
+int BOp_redundant( storage*, reax_list*, real, int, int, int, ivec*, real, rvec*,
+        int, single_body_parameters*, single_body_parameters*,
+        two_body_parameters* );
+
 void BO( reax_system*, control_params*, simulation_data*,
          storage*, reax_list**, output_controls* );
+
 #endif

PuReMD/src/reax_types.h

@@ -1033,9 +1033,9 @@ typedef struct
     real *dDelta_lp, *dDelta_lp_temp;
     real *nlp, *nlp_temp, *Clp, *vlpex;
     rvec *dDeltap_self;
-    int *bond_mark, *done_after;
+    int *bond_mark;
 
-    /* QEq storage */
+    /* charge matrix storage */
     sparse_matrix *H;
     sparse_matrix *L;
     sparse_matrix *U;
@@ -1051,23 +1051,23 @@ typedef struct
     real *y, *g;
     real *hc, *hs;
     real **h, **v;
-    /* GMRES, PIPECG storage */
+    /* GMRES, PIPECG, PIPECR storage */
     real *z;
-    /* CG, PIPECG storage */
+    /* CG, PIPECG, PIPECR storage */
     real *d, *p, *q, *r;
-    /* PIPECG storage */
+    /* PIPECG, PIPECR storage */
     real *m, *n, *u, *w;
     /* dual-CG storage */
     rvec2 *r2, *d2, *q2, *p2;
     /* Taper */
-    real Tap[8]; //Tap7, Tap6, Tap5, Tap4, Tap3, Tap2, Tap1, Tap0;
+    real Tap[8];
 
     /* storage for analysis */
-    int  *mark, *old_mark;
+    int *mark, *old_mark;
     rvec *x_old;
 
     /* storage space for bond restrictions */
-    int  *restricted;
+    int *restricted;
     int **restricted_list;
 
     /* integrator */