sPuReMD: add cutoff-based ILU factor generation.

puremd_rc_1003/sPuReMD/forces.c

@@ -266,9 +266,9 @@ void Init_Forces( reax_system *system, control_params *control,
     int i, j, pj;
     int start_i, end_i;
     int type_i, type_j;
-    int Htop, btop_i, btop_j, num_bonds, num_hbonds;
+    int Htop, H_sp_top, btop_i, btop_j, num_bonds, num_hbonds;
     int ihb, jhb, ihb_top, jhb_top;
-    int flag;
+    int flag, flag_sp;
     real r_ij, r2, self_coef;
     real dr3gamij_1, dr3gamij_3, Tap;
     //real val, dif, base;
@@ -276,7 +276,7 @@ void Init_Forces( reax_system *system, control_params *control,
     real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
     real BO, BO_s, BO_pi, BO_pi2;
     real p_boc1, p_boc2;
-    sparse_matrix *H;
+    sparse_matrix *H, *H_sp;
     list *far_nbrs, *bonds, *hbonds;
     single_body_parameters *sbp_i, *sbp_j;
     two_body_parameters *twbp;
@@ -291,13 +291,19 @@ void Init_Forces( reax_system *system, control_params *control,
     hbonds = *lists + HBONDS;
 
     H = workspace->H;
+    H_sp = workspace->H_sp;
     Htop = 0;
+    H_sp_top = 0;
     num_bonds = 0;
     num_hbonds = 0;
     btop_i = btop_j = 0;
     p_boc1 = system->reaxprm.gp.l[0];
     p_boc2 = system->reaxprm.gp.l[1];
 
+    //TODO: add to control file, remove
+//    control->r_sp_cut = control->r_cut;
+    control->r_sp_cut = 4.0;
+
     for ( i = 0; i < system->N; ++i )
     {
         atom_i = &(system->atoms[i]);
@@ -305,6 +311,7 @@ void Init_Forces( reax_system *system, control_params *control,
         start_i = Start_Index(i, far_nbrs);
         end_i   = End_Index(i, far_nbrs);
         H->start[i] = Htop;
+        H_sp->start[i] = H_sp_top;
         btop_i = End_Index( i, bonds );
         sbp_i = &(system->reaxprm.sbp[type_i]);
         ihb = ihb_top = -1;
@@ -318,17 +325,31 @@ void Init_Forces( reax_system *system, control_params *control,
             atom_j = &(system->atoms[j]);
 
             flag = 0;
-            /*TODO: build H with smaller cutoff, use for preconditioning*/
+            flag_sp = 0;
             if ((data->step - data->prev_steps) % control->reneighbor == 0)
             {
-                if ( nbr_pj->d <= control->r_cut)
+                if ( nbr_pj->d <= control->r_cut )
+                {
                     flag = 1;
-                else flag = 0;
+                    if ( nbr_pj->d <= control->r_sp_cut )
+                    {
+                        flag_sp = 1;
+                    }
+                }
+                else
+                {
+                    flag = 0;
+                    flag_sp = 0;
+                }
             }
             else if ((nbr_pj->d = Sq_Distance_on_T3(atom_i->x, atom_j->x, &(system->box),
                                                     nbr_pj->dvec)) <= SQR(control->r_cut))
             {
-                nbr_pj->d = sqrt(nbr_pj->d);
+                if ( nbr_pj->d <= SQR(control->r_sp_cut))
+                {
+                    flag_sp = 1;
+                }
+                nbr_pj->d = SQRT( nbr_pj->d );
                 flag = 1;
             }
 
@@ -356,6 +377,14 @@ void Init_Forces( reax_system *system, control_params *control,
                 H->entries[Htop].val = self_coef * Tap * EV_to_KCALpMOL / dr3gamij_3;
                 ++Htop;
 
+                /* H_sp matrix entry */
+                if ( flag_sp )
+                {
+                    H_sp->entries[H_sp_top].j = j;
+                    H_sp->entries[H_sp_top].val = self_coef * Tap * EV_to_KCALpMOL / dr3gamij_3;
+                    ++H_sp_top;
+                }
+
                 /* hydrogen bond lists */
                 if ( control->hb_cut > 0 && (ihb == 1 || ihb == 2) &&
                         nbr_pj->d <= control->hb_cut )
@@ -514,6 +543,10 @@ void Init_Forces( reax_system *system, control_params *control,
         H->entries[Htop].val = system->reaxprm.sbp[type_i].eta;
         ++Htop;
 
+        H_sp->entries[H_sp_top].j = i;
+        H_sp->entries[H_sp_top].val = system->reaxprm.sbp[type_i].eta;
+        ++H_sp_top;
+
         Set_End_Index( i, btop_i, bonds );
         if ( ihb == 1 )
             Set_End_Index( workspace->hbond_index[i], ihb_top, hbonds );
@@ -521,8 +554,12 @@ void Init_Forces( reax_system *system, control_params *control,
         //     i, Start_Index( i, bonds ), End_Index( i, bonds ) );
     }
 
+//    printf("Htop = %d\n", Htop);
+//    printf("H_sp_top = %d\n", H_sp_top);
+
     // mark the end of j list
     H->start[i] = Htop;
+    H_sp->start[i] = H_sp_top;
     /* validate lists - decide if reallocation is required! */
     Validate_Lists( workspace, lists,
                     data->step, system->N, H->m, Htop, num_bonds, num_hbonds );

puremd_rc_1003/sPuReMD/init_md.c

@@ -243,6 +243,7 @@ void Init_Workspace( reax_system *system, control_params *control,
 
     /* QEq storage */
     workspace->H        = NULL;
+    workspace->H_sp     = NULL;
     workspace->L        = NULL;
     workspace->U        = NULL;
     workspace->droptol  = (real *) calloc( system->N, sizeof( real ) );
@@ -375,6 +376,11 @@ void Init_Lists( reax_system *system, control_params *control,
                             &Htop, hb_top, bond_top, &num_3body );
 
     Allocate_Matrix( &(workspace->H), system->N, Htop );
+    /* TODO: better estimate for H_sp?
+     *   If so, need to refactor Estimate_Storage_Sizes
+     *   to use various cut-off distances as parameters
+     *   (non-bonded, hydrogen, 3body, etc.) */
+    Allocate_Matrix( &(workspace->H_sp), system->N, Htop );
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "estimated storage - Htop: %d\n", Htop );
     fprintf( stderr, "memory allocated: H = %ldMB\n",

puremd_rc_1003/sPuReMD/mytypes.h

@@ -124,17 +124,17 @@ typedef real rvec[3];
 typedef int  ivec[3];
 typedef real rtensor[3][3];
 
-enum {NVE, NVT, NPT, sNPT, iNPT, ensNR, bNVT};
-enum {FAR_NBRS, NEAR_NBRS, THREE_BODIES, BONDS, OLD_BONDS,
-      HBONDS, DBO, DDELTA, LIST_N
+enum {NVE = 0, NVT = 1, NPT = 2, sNPT = 3, iNPT = 4, ensNR = 5, bNVT = 6};
+enum {FAR_NBRS = 0, NEAR_NBRS = 1, THREE_BODIES = 2, BONDS = 3, OLD_BONDS = 4,
+      HBONDS = 5, DBO = 6, DDELTA = 7, LIST_N = 8
      };
-enum {TYP_VOID, TYP_THREE_BODY, TYP_BOND, TYP_HBOND, TYP_DBO,
-      TYP_DDELTA, TYP_FAR_NEIGHBOR, TYP_NEAR_NEIGHBOR, TYP_N
+enum {TYP_VOID = 0, TYP_THREE_BODY = 1, TYP_BOND = 2, TYP_HBOND = 3, TYP_DBO = 4,
+      TYP_DDELTA = 5, TYP_FAR_NEIGHBOR = 6, TYP_NEAR_NEIGHBOR = 7, TYP_N = 8
      };
-enum {UNKNOWN, WATER};
-enum {NO_ANALYSIS, FRAGMENTS, REACTIONS, NUM_ANALYSIS};
-enum {WRITE_ASCII, WRITE_BINARY, RF_N};
-enum {XYZ, PDB, BGF, ASCII_RESTART, BINARY_RESTART, GF_N};
+enum {UNKNOWN = 0, WATER = 1};
+enum {NO_ANALYSIS = 0, FRAGMENTS = 1, REACTIONS = 2, NUM_ANALYSIS = 3};
+enum {WRITE_ASCII = 0, WRITE_BINARY = 1, RF_N = 2};
+enum {XYZ = 0, PDB = 1, BGF = 2, ASCII_RESTART = 3, BINARY_RESTART = 4, GF_N = 5};
 
 
 
@@ -410,7 +410,7 @@ typedef struct
     int reneighbor;
     real vlist_cut;
     real nbr_cut;
-    real r_cut, r_low; // upper and lower taper
+    real r_cut, r_sp_cut, r_low; // upper and lower taper
     real bo_cut;
     real thb_cut;
     real hb_cut;
@@ -653,9 +653,9 @@ typedef struct
  * See, e.g.,
  *   http://netlib.org/linalg/html_templates/node91.html#SECTION00931100000000000000
  *
- *   m: number of nonzeros (NNZ)
+ *   m: number of nonzeros (NNZ) ALLOCATED
  *   n: number of rows
- *   start: row pointer (last element contains NNZ + 1)
+ *   start: row pointer (last element contains ACTUAL NNZ)
  *   entries: (column index, val) pairs
  * */
 typedef struct
@@ -691,7 +691,7 @@ typedef struct
     rvec *dDeltap_self;
 
     /* QEq storage */
-    sparse_matrix *H, *L, *U;
+    sparse_matrix *H, *H_sp, *L, *U;
     real *droptol;
     real *w;
     real *Hdia_inv;