From a544b061dba11033df1453a7204641712fa76d2d Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Sun, 13 Mar 2016 17:48:24 -0400
Subject: [PATCH] sPuReMD: add cutoff-based ILU factor generation.
---
puremd_rc_1003/sPuReMD/forces.c | 51 +++++++++++++++++++++++++++-----
puremd_rc_1003/sPuReMD/init_md.c | 6 ++++
puremd_rc_1003/sPuReMD/mytypes.h | 26 ++++++++--------
3 files changed, 63 insertions(+), 20 deletions(-)
diff --git a/puremd_rc_1003/sPuReMD/forces.c b/puremd_rc_1003/sPuReMD/forces.c
index d9348b33..cf7bd414 100644
--- a/puremd_rc_1003/sPuReMD/forces.c
+++ b/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 );
diff --git a/puremd_rc_1003/sPuReMD/init_md.c b/puremd_rc_1003/sPuReMD/init_md.c
index d6f15a1e..fb95ee22 100644
--- a/puremd_rc_1003/sPuReMD/init_md.c
+++ b/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",
diff --git a/puremd_rc_1003/sPuReMD/mytypes.h b/puremd_rc_1003/sPuReMD/mytypes.h
index 5618418c..9bef9190 100644
--- a/puremd_rc_1003/sPuReMD/mytypes.h
+++ b/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;
--
GitLab