From d7897f19270d7e7db4dba9dbfffb279102680468 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Thu, 25 Oct 2018 05:37:56 -0700
Subject: [PATCH] PuReMD: refactor to allow selection of far range interaction
Verlet list format at runtime (half or full). Similar change for symmetric
charge matrix format (half or full). For SAI, toggle between half and full
depending on whether it is a refactoring step.
---
PuReMD/src/allocate.c | 39 ++++++--
PuReMD/src/allocate.h | 4 +-
PuReMD/src/bond_orders.c | 115 ++++++++++------------
PuReMD/src/bond_orders.h | 5 +-
PuReMD/src/forces.c | 188 +++++++++++++++++++-----------------
PuReMD/src/init_md.c | 55 +++++++----
PuReMD/src/integrate.c | 64 ++++++++++++
PuReMD/src/linear_solvers.c | 164 +++++++++++++++----------------
PuReMD/src/list.c | 4 +-
PuReMD/src/list.h | 2 +-
PuReMD/src/neighbors.c | 45 ++++-----
PuReMD/src/neighbors.h | 2 +-
PuReMD/src/nonbonded.c | 16 ++-
PuReMD/src/qEq.c | 4 +-
PuReMD/src/reax_defs.h | 126 ++++++++++++------------
PuReMD/src/reax_types.h | 26 ++++-
16 files changed, 488 insertions(+), 371 deletions(-)
diff --git a/PuReMD/src/allocate.c b/PuReMD/src/allocate.c
index ca67bea5..08f45586 100644
--- a/PuReMD/src/allocate.c
+++ b/PuReMD/src/allocate.c
@@ -394,8 +394,12 @@ int Allocate_Workspace( reax_system *system, control_params *control,
void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs,
MPI_Comm comm )
{
+ int format;
+
+ format = far_nbrs->format;
+
Delete_List( far_nbrs, comm );
- if (!Make_List( n, num_intrs, TYP_FAR_NEIGHBOR, far_nbrs, comm ))
+ if (!Make_List( n, num_intrs, TYP_FAR_NEIGHBOR, format, far_nbrs, comm ))
{
fprintf(stderr, "Problem in initializing far nbrs list. Terminating!\n");
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -403,7 +407,8 @@ void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs,
}
-int Allocate_Matrix( sparse_matrix **pH, int cap, int m, MPI_Comm comm )
+int Allocate_Matrix( sparse_matrix **pH, int cap, int m,
+ int format, MPI_Comm comm )
{
sparse_matrix *H;
@@ -412,6 +417,7 @@ int Allocate_Matrix( sparse_matrix **pH, int cap, int m, MPI_Comm comm )
H = *pH;
H->cap = cap;
H->m = m;
+ H->format = format;
H->start = (int*) smalloc( sizeof(int) * cap, "matrix_start", comm );
H->end = (int*) smalloc( sizeof(int) * cap, "matrix_end", comm );
H->entries = (sparse_matrix_entry*)
@@ -419,7 +425,8 @@ int Allocate_Matrix( sparse_matrix **pH, int cap, int m, MPI_Comm comm )
return SUCCESS;
}
-int Allocate_Matrix2( sparse_matrix **pH, int n, int cap, int m, MPI_Comm comm )
+int Allocate_Matrix2( sparse_matrix **pH, int n, int cap, int m,
+ int format, MPI_Comm comm )
{
sparse_matrix *H;
@@ -429,6 +436,7 @@ int Allocate_Matrix2( sparse_matrix **pH, int n, int cap, int m, MPI_Comm comm )
H->n = n;
H->cap = cap;
H->m = m;
+ H->format = format;
H->start = (int*) smalloc( sizeof(int) * cap, "matrix_start", comm );
H->end = (int*) smalloc( sizeof(int) * cap, "matrix_end", comm );
H->entries = (sparse_matrix_entry*)
@@ -450,8 +458,12 @@ void Deallocate_Matrix( sparse_matrix *H )
int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name,
MPI_Comm comm )
{
+ int format;
+
+ format = (*H)->format;
+
Deallocate_Matrix( *H );
- if ( !Allocate_Matrix( H, n, m, comm ) )
+ if ( !Allocate_Matrix( H, n, m, format, comm ) )
{
fprintf(stderr, "not enough space for %s matrix. terminating!\n", name);
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -469,7 +481,9 @@ int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name,
int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
MPI_Comm comm )
{
- int i, id, total_hbonds;
+ int i, id, total_hbonds, format;
+
+ format = hbonds->format;
total_hbonds = 0;
for ( i = 0; i < system->n; ++i )
@@ -483,7 +497,7 @@ int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
total_hbonds = (int)(MAX( total_hbonds * SAFER_ZONE, MIN_CAP * MIN_HBONDS ));
Delete_List( hbonds, comm );
- if ( !Make_List( system->Hcap, total_hbonds, TYP_HBOND, hbonds, comm ) )
+ if ( !Make_List( system->Hcap, total_hbonds, TYP_HBOND, format, hbonds, comm ) )
{
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -496,7 +510,9 @@ int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
int *total_bonds, int *est_3body, MPI_Comm comm )
{
- int i;
+ int i, format;
+
+ format = bonds->format;
*total_bonds = 0;
*est_3body = 0;
@@ -510,7 +526,7 @@ int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
*total_bonds = (int)(MAX( *total_bonds * SAFE_ZONE, MIN_CAP * MIN_BONDS ));
Delete_List( bonds, comm );
- if (!Make_List(system->total_cap, *total_bonds, TYP_BOND, bonds, comm))
+ if (!Make_List(system->total_cap, *total_bonds, TYP_BOND, format, bonds, comm))
{
fprintf( stderr, "not enough space for bonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -746,7 +762,7 @@ void ReAllocate( reax_system *system, control_params *control,
{
int i, j, k, p;
int num_bonds, est_3body, nflag, Nflag, Hflag, mpi_flag, ret, total_send;
- int renbr;
+ int renbr, format;
reallocate_data *realloc;
reax_list *far_nbrs;
sparse_matrix *H;
@@ -928,6 +944,9 @@ void ReAllocate( reax_system *system, control_params *control,
(int)(realloc->num_3body * sizeof(three_body_interaction_data) /
(1024 * 1024)) );
#endif
+
+ format = lists[THREE_BODIES]->format;
+
Delete_List( lists[THREE_BODIES], comm );
if ( num_bonds == -1 )
@@ -936,7 +955,7 @@ void ReAllocate( reax_system *system, control_params *control,
realloc->num_3body = (int)(MAX(realloc->num_3body * SAFE_ZONE, MIN_3BODIES));
if ( !Make_List( num_bonds, realloc->num_3body, TYP_THREE_BODY,
- lists[THREE_BODIES], comm ) )
+ format, lists[THREE_BODIES], comm ) )
{
fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
MPI_Abort( comm, CANNOT_INITIALIZE );
diff --git a/PuReMD/src/allocate.h b/PuReMD/src/allocate.h
index e4e7a86a..5458f29c 100644
--- a/PuReMD/src/allocate.h
+++ b/PuReMD/src/allocate.h
@@ -38,9 +38,9 @@ void Deallocate_Grid( grid* );
int Allocate_MPI_Buffers( mpi_datatypes*, int, neighbor_proc*, char* );
-int Allocate_Matrix( sparse_matrix**, int, int, MPI_Comm );
+int Allocate_Matrix( sparse_matrix**, int, int, int, MPI_Comm );
-int Allocate_Matrix2( sparse_matrix**, int, int, int, MPI_Comm );
+int Allocate_Matrix2( sparse_matrix**, int, int, int, int, MPI_Comm );
void Deallocate_Matrix( sparse_matrix * );
diff --git a/PuReMD/src/bond_orders.c b/PuReMD/src/bond_orders.c
index 286e15c0..b1e9b228 100644
--- a/PuReMD/src/bond_orders.c
+++ b/PuReMD/src/bond_orders.c
@@ -664,8 +664,8 @@ void Add_dBond_to_Forces( int i, int pj,
int BOp( storage *workspace, reax_list *bonds, real bo_cut,
int i, int btop_i, far_neighbor_data *nbr_pj,
- single_body_parameters *sbp_i, single_body_parameters *sbp_j,
- two_body_parameters *twbp )
+ int far_nbr_list_format, single_body_parameters *sbp_i,
+ single_body_parameters *sbp_j, two_body_parameters *twbp )
{
int j;
real r2, C12, C34, C56;
@@ -673,11 +673,9 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
real BO, BO_s, BO_pi, BO_pi2;
bond_data *ibond;
bond_order_data *bo_ij;
-#if defined(HALF_LIST)
int btop_j;
bond_data *jbond;
bond_order_data *bo_ji;
-#endif
j = nbr_pj->nbr;
r2 = SQR(nbr_pj->d);
@@ -710,46 +708,43 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
{
/****** bonds i-j and j-i ******/
ibond = &( bonds->bond_list[btop_i] );
-#if defined(HALF_LIST)
- btop_j = End_Index( j, bonds );
- jbond = &(bonds->bond_list[btop_j]);
-#endif
+ if ( far_nbr_list_format == HALF_LIST )
+ {
+ btop_j = End_Index( j, bonds );
+ jbond = &(bonds->bond_list[btop_j]);
+ }
-#if defined(HALF_LIST)
- ibond->nbr = j;
- ibond->d = nbr_pj->d;
- rvec_Copy( ibond->dvec, nbr_pj->dvec );
- ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
- ibond->dbond_index = btop_i;
- ibond->sym_index = btop_j;
- jbond->nbr = i;
- jbond->d = nbr_pj->d;
- rvec_Scale( jbond->dvec, -1.0, nbr_pj->dvec );
- ivec_Scale( jbond->rel_box, -1.0, nbr_pj->rel_box );
- jbond->dbond_index = btop_i;
- jbond->sym_index = btop_i;
-
- Set_End_Index( j, btop_j + 1, bonds );
-#else
ibond->nbr = j;
ibond->d = nbr_pj->d;
rvec_Copy( ibond->dvec, nbr_pj->dvec );
ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
ibond->dbond_index = btop_i;
-#endif
+ if ( far_nbr_list_format == HALF_LIST )
+ {
+ ibond->sym_index = btop_j;
+ jbond->nbr = i;
+ jbond->d = nbr_pj->d;
+ rvec_Scale( jbond->dvec, -1.0, nbr_pj->dvec );
+ ivec_Scale( jbond->rel_box, -1.0, nbr_pj->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;
-#if defined(HALF_LIST)
- 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;
-#endif
+ if ( far_nbr_list_format == HALF_LIST )
+ {
+ 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;
@@ -758,36 +753,31 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
/* 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 */
-#if defined(HALF_LIST)
- 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 );
-#else
- 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);
-#endif
+ 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 );
+ if ( far_nbr_list_format == HALF_LIST )
+ {
+ 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 */
-#if defined(HALF_LIST)
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 );
-#else
- 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 );
-#endif
+ if ( far_nbr_list_format == HALF_LIST )
+ {
+ rvec_Scale( bo_ji->dBOp, -1.0, bo_ij->dBOp );
+ }
rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
-#if defined(HALF_LIST)
- rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );
-#endif
+ if ( far_nbr_list_format == HALF_LIST )
+ {
+ rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );
+ }
bo_ij->BO_s -= bo_cut;
bo_ij->BO -= bo_cut;
@@ -795,14 +785,15 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
bo_ij->Cdbo = 0.0;
bo_ij->Cdbopi = 0.0;
bo_ij->Cdbopi2 = 0.0;
-#if defined(HALF_LIST)
- 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;
-#endif
+ if ( far_nbr_list_format == HALF_LIST )
+ {
+ 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;
}
diff --git a/PuReMD/src/bond_orders.h b/PuReMD/src/bond_orders.h
index 1975e20b..a692fccb 100644
--- a/PuReMD/src/bond_orders.h
+++ b/PuReMD/src/bond_orders.h
@@ -52,8 +52,9 @@ void Add_dDelta_to_Forces( reax_system *, reax_list**, int, real );
void Add_dBond_to_Forces( int, int, storage*, reax_list** );
void Add_dBond_to_Forces_NPT( int, int, simulation_data*,
storage*, reax_list** );
-int BOp(storage*, reax_list*, real, int, int, far_neighbor_data*,
- single_body_parameters*, single_body_parameters*, two_body_parameters*);
+int BOp( storage*, reax_list*, real, int, int, far_neighbor_data*,
+ int, single_body_parameters*, single_body_parameters*,
+ two_body_parameters* );
void BO( reax_system*, control_params*, simulation_data*,
storage*, reax_list**, output_controls* );
#endif
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 49e747a4..69ab40e7 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -336,12 +336,9 @@ void Init_Forces( reax_system *system, control_params *control,
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
-#if defined(HALF_LIST)
int jhb_top;
-#else
int start_j, end_j;
int btop_j;
-#endif
far_nbrs = lists[FAR_NBRS];
bonds = lists[BONDS];
@@ -369,13 +366,16 @@ void Init_Forces( reax_system *system, control_params *control,
type_i = atom_i->type;
start_i = Start_Index(i, far_nbrs);
end_i = End_Index(i, far_nbrs);
-#if defined(HALF_LIST)
- /* start at end because other atoms
- * can add to this atom's list (half-list) */
- btop_i = End_Index( i, bonds );
-#else
- btop_i = Start_Index( i, bonds );
-#endif
+ if ( far_nbrs->format == HALF_LIST )
+ {
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
+ btop_i = End_Index( i, bonds );
+ }
+ else if ( far_nbrs->format == FULL_LIST )
+ {
+ btop_i = Start_Index( i, bonds );
+ }
sbp_i = &(system->reax_param.sbp[type_i]);
if ( i < system->n )
@@ -403,13 +403,16 @@ void Init_Forces( reax_system *system, control_params *control,
ihb = sbp_i->p_hbond;
if ( ihb == 1 )
{
-#if defined(HALF_LIST)
- /* start at end because other atoms
- * can add to this atom's list (half-list) */
- ihb_top = End_Index( atom_i->Hindex, hbonds );
-#else
- ihb_top = Start_Index( atom_i->Hindex, hbonds );
-#endif
+ if ( far_nbrs->format == HALF_LIST )
+ {
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
+ ihb_top = End_Index( atom_i->Hindex, hbonds );
+ }
+ else if ( far_nbrs->format == FULL_LIST )
+ {
+ ihb_top = Start_Index( atom_i->Hindex, hbonds );
+ }
}
else
{
@@ -457,9 +460,9 @@ void Init_Forces( reax_system *system, control_params *control,
if ( local )
{
/* H matrix entry */
-#if defined(HALF_LIST)
- if ( j < system->n || atom_i->orig_id < atom_j->orig_id ) //tryQEq||1
-#endif
+ if ( (far_nbrs->format == HALF_LIST
+ && (j < system->n || atom_i->orig_id < atom_j->orig_id))
+ || far_nbrs->format == FULL_LIST )
{
H->entries[Htop].j = j;
if ( control->tabulate == 0 )
@@ -482,9 +485,9 @@ void Init_Forces( reax_system *system, control_params *control,
++ihb_top;
++num_hbonds;
}
-#if defined(HALF_LIST)
/* only add to list for local j (far nbrs is half-list) */
- else if ( j < system->n && ihb == 2 && jhb == 1 )
+ else if ( far_nbrs->format == HALF_LIST
+ && (j < system->n && ihb == 2 && jhb == 1) )
{
jhb_top = End_Index( atom_j->Hindex, hbonds );
hbonds->hbond_list[jhb_top].nbr = i;
@@ -493,7 +496,6 @@ void Init_Forces( reax_system *system, control_params *control,
Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
++num_hbonds;
}
-#endif
}
}
@@ -501,7 +503,8 @@ void Init_Forces( reax_system *system, control_params *control,
if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
nbr_pj->d <= control->bond_cut &&
BOp( workspace, bonds, control->bo_cut,
- i , btop_i, nbr_pj, sbp_i, sbp_j, twbp ) )
+ i, btop_i, nbr_pj, far_nbrs->format,
+ sbp_i, sbp_j, twbp ) )
{
num_bonds += 2;
++btop_i;
@@ -525,30 +528,31 @@ void Init_Forces( reax_system *system, control_params *control,
}
}
-#if !defined(HALF_LIST)
- /* set sym_index for bonds list (far_nbrs full list) */
- for ( i = 0; i < system->N; ++i )
+ if ( far_nbrs->format == FULL_LIST )
{
- start_i = Start_Index( i, bonds );
- end_i = End_Index( i, bonds );
-
- for ( btop_i = start_i; btop_i < end_i; ++btop_i )
+ /* set sym_index for bonds list (far_nbrs full list) */
+ for ( i = 0; i < system->N; ++i )
{
- j = bonds->bond_list[btop_i].nbr;
- start_j = Start_Index( j, bonds );
- end_j = End_Index( j, bonds );
+ start_i = Start_Index( i, bonds );
+ end_i = End_Index( i, bonds );
- for ( btop_j = start_j; btop_j < end_j; ++btop_j )
+ for ( btop_i = start_i; btop_i < end_i; ++btop_i )
{
- if ( bonds->bond_list[btop_j].nbr == i )
+ j = bonds->bond_list[btop_i].nbr;
+ start_j = Start_Index( j, bonds );
+ end_j = End_Index( j, bonds );
+
+ for ( btop_j = start_j; btop_j < end_j; ++btop_j )
{
- bonds->bond_list[btop_i].sym_index = btop_j;
- break;
+ if ( bonds->bond_list[btop_j].nbr == i )
+ {
+ bonds->bond_list[btop_i].sym_index = btop_j;
+ break;
+ }
}
}
}
}
-#endif
workspace->realloc.Htop = Htop;
workspace->realloc.num_bonds = num_bonds;
@@ -595,12 +599,9 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
-#if defined(HALF_LIST)
int jhb_top;
-#else
int start_j, end_j;
int btop_j;
-#endif
far_nbrs = lists[FAR_NBRS];
bonds = lists[BONDS];
@@ -625,13 +626,16 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
type_i = atom_i->type;
start_i = Start_Index(i, far_nbrs);
end_i = End_Index(i, far_nbrs);
-#if defined(HALF_LIST)
- /* start at end because other atoms
- * can add to this atom's list (half-list) */
- btop_i = End_Index( i, bonds );
-#else
- btop_i = Start_Index( i, bonds );
-#endif
+ if ( far_nbrs->format == HALF_LIST )
+ {
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
+ btop_i = End_Index( i, bonds );
+ }
+ else if ( far_nbrs->format == FULL_LIST )
+ {
+ btop_i = Start_Index( i, bonds );
+ }
sbp_i = &(system->reax_param.sbp[type_i]);
if ( i < system->n )
@@ -652,13 +656,16 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
ihb = sbp_i->p_hbond;
if ( ihb == 1 )
{
-#if defined(HALF_LIST)
- /* start at end because other atoms
- * can add to this atom's list (half-list) */
- ihb_top = End_Index( atom_i->Hindex, hbonds );
-#else
- ihb_top = Start_Index( atom_i->Hindex, hbonds );
-#endif
+ if ( far_nbrs->format == HALF_LIST )
+ {
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
+ ihb_top = End_Index( atom_i->Hindex, hbonds );
+ }
+ else if ( far_nbrs->format == FULL_LIST )
+ {
+ ihb_top = Start_Index( atom_i->Hindex, hbonds );
+ }
}
else
{
@@ -719,9 +726,9 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
++ihb_top;
++num_hbonds;
}
-#if defined(HALF_LIST)
/* only add to list for local j (far nbrs is half-list) */
- else if ( j < system->n && ihb == 2 && jhb == 1 )
+ else if ( far_nbrs->format == HALF_LIST
+ && (j < system->n && ihb == 2 && jhb == 1) )
{
jhb_top = End_Index( atom_j->Hindex, hbonds );
hbonds->hbond_list[jhb_top].nbr = i;
@@ -730,7 +737,6 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
++num_hbonds;
}
-#endif
}
}
@@ -739,7 +745,8 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
if ( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
nbr_pj->d <= control->bond_cut &&
BOp( workspace, bonds, control->bo_cut,
- i , btop_i, nbr_pj, sbp_i, sbp_j, twbp ) )
+ i, btop_i, nbr_pj, far_nbrs->format,
+ sbp_i, sbp_j, twbp ) )
{
num_bonds += 2;
++btop_i;
@@ -763,30 +770,31 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
}
-#if !defined(HALF_LIST)
- /* set sym_index for bonds list (far_nbrs full list) */
- for ( i = 0; i < system->N; ++i )
+ if ( far_nbrs->format == FULL_LIST )
{
- start_i = Start_Index( i, bonds );
- end_i = End_Index( i, bonds );
-
- for ( btop_i = start_i; btop_i < end_i; ++btop_i )
+ /* set sym_index for bonds list (far_nbrs full list) */
+ for ( i = 0; i < system->N; ++i )
{
- j = bonds->bond_list[btop_i].nbr;
- start_j = Start_Index( j, bonds );
- end_j = End_Index( j, bonds );
+ start_i = Start_Index( i, bonds );
+ end_i = End_Index( i, bonds );
- for ( btop_j = start_j; btop_j < end_j; ++btop_j )
+ for ( btop_i = start_i; btop_i < end_i; ++btop_i )
{
- if ( bonds->bond_list[btop_j].nbr == i )
+ j = bonds->bond_list[btop_i].nbr;
+ start_j = Start_Index( j, bonds );
+ end_j = End_Index( j, bonds );
+
+ for ( btop_j = start_j; btop_j < end_j; ++btop_j )
{
- bonds->bond_list[btop_i].sym_index = btop_j;
- break;
+ if ( bonds->bond_list[btop_j].nbr == i )
+ {
+ bonds->bond_list[btop_i].sym_index = btop_j;
+ break;
+ }
}
}
}
}
-#endif
workspace->realloc.num_bonds = num_bonds;
workspace->realloc.num_hbonds = num_hbonds;
@@ -823,10 +831,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
- reax_atom *atom_i;
-#if defined(HALF_LIST)
- reax_atom *atom_j;
-#endif
+ reax_atom *atom_i, *atom_j;
far_nbrs = lists[FAR_NBRS];
*Htop = 0;
@@ -860,9 +865,10 @@ void Estimate_Storages( reax_system *system, control_params *control,
{
nbr_pj = &( far_nbrs->far_nbr_list[pj] );
j = nbr_pj->nbr;
-#if defined(HALF_LIST)
- atom_j = &(system->my_atoms[j]);
-#endif
+ if ( far_nbrs->format == HALF_LIST )
+ {
+ atom_j = &(system->my_atoms[j]);
+ }
if (nbr_pj->d <= cutoff)
{
@@ -873,9 +879,9 @@ void Estimate_Storages( reax_system *system, control_params *control,
if ( local )
{
-#if defined(HALF_LIST)
- if ( j < system->n || atom_i->orig_id < atom_j->orig_id ) //tryQEq ||1
-#endif
+ if ( (far_nbrs->format == HALF_LIST
+ && (j < system->n || atom_i->orig_id < atom_j->orig_id))
+ || far_nbrs->format == FULL_LIST )
{
++(*Htop);
}
@@ -890,12 +896,11 @@ void Estimate_Storages( reax_system *system, control_params *control,
++hb_top[i];
}
/* only add to list for local j (far nbrs is half-list) */
-#if defined(HALF_LIST)
- else if ( j < system->n && ihb == 2 && jhb == 1 )
+ else if ( far_nbrs->format == HALF_LIST
+ && (j < system->n && ihb == 2 && jhb == 1) )
{
++hb_top[j];
}
-#endif
}
}
@@ -931,9 +936,10 @@ void Estimate_Storages( reax_system *system, control_params *control,
if ( BO >= control->bo_cut )
{
++bond_top[i];
-#if defined(HALF_LIST)
- ++bond_top[j];
-#endif
+ if ( far_nbrs->format == HALF_LIST )
+ {
+ ++bond_top[j];
+ }
}
}
}
diff --git a/PuReMD/src/init_md.c b/PuReMD/src/init_md.c
index ef991d27..4fac0f1d 100644
--- a/PuReMD/src/init_md.c
+++ b/PuReMD/src/init_md.c
@@ -568,7 +568,7 @@ int Init_Lists( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
mpi_datatypes *mpi_data, char *msg )
{
- int i, num_nbrs;
+ int i, num_nbrs, far_nbr_list_format, cm_format;
int total_hbonds, total_bonds, bond_cap, num_3body, cap_3body, Htop;
int *hb_top, *bond_top;
MPI_Comm comm;
@@ -579,10 +579,23 @@ int Init_Lists( reax_system *system, control_params *control,
#endif
comm = mpi_data->world;
+
+ if ( control->cm_solver_pre_comp_type == SAI_PC )
+ {
+ far_nbr_list_format = FULL_LIST;
+ cm_format = SYM_FULL_MATRIX;
+ }
+ else
+ {
+ far_nbr_list_format = HALF_LIST;
+ cm_format = SYM_HALF_MATRIX;
+ }
+
//for( i = 0; i < MAX_NBRS; ++i ) nrecv[i] = system->my_nbrs[i].est_recv;
//system->N = SendRecv( system, mpi_data, mpi_data->boundary_atom_type, nrecv,
// Sort_Boundary_Atoms, Unpack_Exchange_Message, 1 );
- num_nbrs = Estimate_NumNeighbors( system, lists );
+
+ num_nbrs = Estimate_NumNeighbors( system, lists, far_nbr_list_format );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: after est_nbrs - local_cap=%d, total_cap=%d\n",
@@ -590,11 +603,12 @@ int Init_Lists( reax_system *system, control_params *control,
#endif
if ( !Make_List( system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR,
- lists[FAR_NBRS], comm ) )
+ far_nbr_list_format, lists[FAR_NBRS], comm ) )
{
fprintf(stderr, "Problem in initializing far nbrs list. Terminating!\n");
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated far_nbrs: num_far=%d, space=%dMB\n",
system->my_rank, num_nbrs,
@@ -620,14 +634,18 @@ int Init_Lists( reax_system *system, control_params *control,
Estimate_Storages( system, control, lists,
&Htop, hb_top, bond_top, &num_3body, comm );
- Allocate_Matrix( &(workspace->H), system->local_cap, Htop, comm );
+ Allocate_Matrix( &(workspace->H), system->local_cap, Htop,
+ cm_format, comm );
workspace->L = NULL;
workspace->U = NULL;
//TODO: uncomment for SAI
- // Allocate_Matrix2( &(workspace->H_spar_patt), workspace->H->n, system->local_cap, workspace->H->m, comm );
- // Allocate_Matrix( &(workspace->H_spar_patt_full), workspace->H->n, 2 * workspace->H->m - workspace->H->n );
- // Allocate_Matrix2( &(workspace->H_app_inv), workspace->H->n, system->local_cap, workspace->H->m, comm );
+// Allocate_Matrix2( &(workspace->H_spar_patt), workspace->H->n, system->local_cap,
+// workspace->H->m, SYM_HALF_MATRIX, comm );
+// Allocate_Matrix( &(workspace->H_spar_patt_full), workspace->H->n,
+// 2 * workspace->H->m - workspace->H->n, SYM_FULL_MATRIX, comm );
+// Allocate_Matrix2( &(workspace->H_app_inv), workspace->H->n, system->local_cap,
+// workspace->H->m, SYM_FULL_MATRIX, comm );
workspace->H_spar_patt = NULL;
workspace->H_app_inv = NULL;
@@ -649,7 +667,7 @@ int Init_Lists( reax_system *system, control_params *control,
total_hbonds = MAX( total_hbonds * SAFER_ZONE, MIN_CAP * MIN_HBONDS );
if ( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
- lists[HBONDS], comm ) )
+ HALF_LIST, lists[HBONDS], comm ) )
{
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -673,7 +691,7 @@ int Init_Lists( reax_system *system, control_params *control,
bond_cap = MAX( total_bonds * SAFE_ZONE, MIN_CAP * MIN_BONDS );
if ( !Make_List( system->total_cap, bond_cap, TYP_BOND,
- lists[BONDS], comm ) )
+ HALF_LIST, lists[BONDS], comm ) )
{
fprintf( stderr, "not enough space for bonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -687,7 +705,7 @@ int Init_Lists( reax_system *system, control_params *control,
/* 3bodies list */
cap_3body = MAX( num_3body * SAFE_ZONE, MIN_3BODIES );
if ( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
- lists[THREE_BODIES], comm ) )
+ HALF_LIST, lists[THREE_BODIES], comm ) )
{
fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -700,7 +718,7 @@ int Init_Lists( reax_system *system, control_params *control,
#if defined(TEST_FORCES)
if ( !Make_List( system->total_cap, bond_cap * 8, TYP_DDELTA,
- lists[DDELTAS], comm ) )
+ HALF_LIST, lists[DDELTAS], comm ) )
{
fprintf( stderr, "Problem in initializing dDelta list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -709,7 +727,7 @@ int Init_Lists( reax_system *system, control_params *control,
system->my_rank, bond_cap * 30,
bond_cap * 8 * sizeof(dDelta_data) / (1024 * 1024) );
- if ( !Make_List( bond_cap, bond_cap * 50, TYP_DBO, lists[DBOS], comm ) )
+ if ( !Make_List( bond_cap, bond_cap * 50, TYP_DBO, HALF_LIST, lists[DBOS], comm ) )
{
fprintf( stderr, "Problem in initializing dBO list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -724,6 +742,8 @@ int Init_Lists( reax_system *system, control_params *control,
return SUCCESS;
}
+
+
#elif defined(LAMMPS_REAX)
int Init_Lists( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
@@ -736,6 +756,7 @@ int Init_Lists( reax_system *system, control_params *control,
MPI_Comm comm;
comm = mpi_data->world;
+
bond_top = (int*) calloc( system->total_cap, sizeof(int) );
hb_top = (int*) calloc( system->local_cap, sizeof(int) );
Estimate_Storages( system, control, lists,
@@ -753,7 +774,7 @@ int Init_Lists( reax_system *system, control_params *control,
total_hbonds = (int)(MAX( total_hbonds * SAFER_ZONE, MIN_CAP * MIN_HBONDS ));
if ( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
- lists[HBONDS], comm ) )
+ HALF_LIST, lists[HBONDS], comm ) )
{
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -777,7 +798,7 @@ int Init_Lists( reax_system *system, control_params *control,
bond_cap = (int)(MAX( total_bonds * SAFE_ZONE, MIN_CAP * MIN_BONDS ));
if ( !Make_List( system->total_cap, bond_cap, TYP_BOND,
- lists[BONDS], comm ) )
+ HALF_LIST, lists[BONDS], comm ) )
{
fprintf( stderr, "not enough space for bonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -791,7 +812,7 @@ int Init_Lists( reax_system *system, control_params *control,
/* 3bodies list */
cap_3body = (int)(MAX( num_3body * SAFE_ZONE, MIN_3BODIES ));
if ( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
- lists[THREE_BODIES], comm ) )
+ HALF_LIST, lists[THREE_BODIES], comm ) )
{
fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -804,7 +825,7 @@ int Init_Lists( reax_system *system, control_params *control,
#if defined(TEST_FORCES)
if ( !Make_List( system->total_cap, bond_cap * 8, TYP_DDELTA,
- lists[DDELTAS], comm ) )
+ HALF_LIST, lists[DDELTAS], comm ) )
{
fprintf( stderr, "Problem in initializing dDelta list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -813,7 +834,7 @@ int Init_Lists( reax_system *system, control_params *control,
system->my_rank, bond_cap * 30,
bond_cap * 8 * sizeof(dDelta_data) / (1024 * 1024) );
- if ( !Make_List( bond_cap, bond_cap * 50, TYP_DBO, lists[DBOS], comm ) )
+ if ( !Make_List( bond_cap, bond_cap * 50, TYP_DBO, HALF_LIST, lists[DBOS], comm ) )
{
fprintf( stderr, "Problem in initializing dBO list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
diff --git a/PuReMD/src/integrate.c b/PuReMD/src/integrate.c
index 3672651c..c3b2719c 100644
--- a/PuReMD/src/integrate.c
+++ b/PuReMD/src/integrate.c
@@ -51,6 +51,22 @@ void Velocity_Verlet_NVE( reax_system* system, control_params* control,
dt_sqr = SQR(dt);
steps = data->step - data->prev_steps;
renbr = (steps % control->reneighbor == 0);
+ if ( control->cm_solver_pre_comp_type == SAI_PC )
+ {
+ /* HACK: currently required that preconditioner (re)computation step
+ * and reneighbor step (i.e., (re)construct far nbr list)
+ * are the same value, so use reneighbor for now */
+ if ( renbr )
+ {
+ lists[FAR_NBRS]->format = FULL_LIST;
+ workspace->H->format = SYM_FULL_MATRIX;
+ }
+ else
+ {
+ lists[FAR_NBRS]->format = HALF_LIST;
+ workspace->H->format = SYM_HALF_MATRIX;
+ }
+ }
for ( i = 0; i < system->n; i++ )
{
@@ -114,6 +130,22 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system* system,
therm = &( data->therm );
steps = data->step - data->prev_steps;
renbr = (steps % control->reneighbor == 0);
+ if ( control->cm_solver_pre_comp_type == SAI_PC )
+ {
+ /* HACK: currently required that preconditioner (re)computation step
+ * and reneighbor step (i.e., (re)construct far nbr list)
+ * are the same value, so use reneighbor for now */
+ if ( renbr )
+ {
+ lists[FAR_NBRS]->format = FULL_LIST;
+ workspace->H->format = SYM_FULL_MATRIX;
+ }
+ else
+ {
+ lists[FAR_NBRS]->format = HALF_LIST;
+ workspace->H->format = SYM_HALF_MATRIX;
+ }
+ }
for ( i = 0; i < system->n; i++ )
{
@@ -209,6 +241,22 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
dt = control->dt;
steps = data->step - data->prev_steps;
renbr = (steps % control->reneighbor == 0);
+ if ( control->cm_solver_pre_comp_type == SAI_PC )
+ {
+ /* HACK: currently required that preconditioner (re)computation step
+ * and reneighbor step (i.e., (re)construct far nbr list)
+ * are the same value, so use reneighbor for now */
+ if ( renbr )
+ {
+ lists[FAR_NBRS]->format = FULL_LIST;
+ workspace->H->format = SYM_FULL_MATRIX;
+ }
+ else
+ {
+ lists[FAR_NBRS]->format = HALF_LIST;
+ workspace->H->format = SYM_HALF_MATRIX;
+ }
+ }
/* velocity verlet, 1st part */
for ( i = 0; i < system->n; i++ )
@@ -300,6 +348,22 @@ void Velocity_Verlet_Berendsen_NPT( reax_system* system,
dt = control->dt;
steps = data->step - data->prev_steps;
renbr = (steps % control->reneighbor == 0);
+ if ( control->cm_solver_pre_comp_type == SAI_PC )
+ {
+ /* HACK: currently required that preconditioner (re)computation step
+ * and reneighbor step (i.e., (re)construct far nbr list)
+ * are the same value, so use reneighbor for now */
+ if ( renbr )
+ {
+ lists[FAR_NBRS]->format = FULL_LIST;
+ workspace->H->format = SYM_FULL_MATRIX;
+ }
+ else
+ {
+ lists[FAR_NBRS]->format = HALF_LIST;
+ workspace->H->format = SYM_HALF_MATRIX;
+ }
+ }
/* velocity verlet, 1st part */
for ( i = 0; i < system->n; i++ )
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index 2d09f7da..a404bdb1 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -150,13 +150,15 @@ real setup_sparse_approx_inverse( reax_system *system, simulation_data *data, st
if ( *A_spar_patt == NULL )
{
- Allocate_Matrix2(A_spar_patt, A->n, system->local_cap, A->m, comm );
+ Allocate_Matrix2( A_spar_patt, A->n, system->local_cap, A->m,
+ A->format, comm );
}
else /*if ( (*A_spar_patt)->m < A->m )*/
{
Deallocate_Matrix( *A_spar_patt );
- Allocate_Matrix2( A_spar_patt, A->n, system->local_cap, A->m, comm );
+ Allocate_Matrix2( A_spar_patt, A->n, system->local_cap, A->m,
+ A->format, comm );
}
n_local = 0;
@@ -425,8 +427,10 @@ real setup_sparse_approx_inverse( reax_system *system, simulation_data *data, st
return Get_Timing_Info( start );
}
-real sparse_approx_inverse(reax_system *system, simulation_data *data, storage *workspace, mpi_datatypes *mpi_data,
- sparse_matrix *A, sparse_matrix *A_spar_patt, sparse_matrix **A_app_inv, int nprocs )
+real sparse_approx_inverse( reax_system *system, simulation_data *data,
+ storage *workspace, mpi_datatypes *mpi_data,
+ sparse_matrix *A, sparse_matrix *A_spar_patt,
+ sparse_matrix **A_app_inv, int nprocs )
{
int N, M, d_i, d_j;
int i, k, pj, j_temp;
@@ -434,13 +438,11 @@ real sparse_approx_inverse(reax_system *system, simulation_data *data, storage *
lapack_int m, n, nrhs, lda, ldb, info;
int *pos_x, *X;
real *e_j, *dense_matrix;
-
int cnt, scale;
reax_atom *atom;
int *row_nnz;
int **j_list;
real **val_list;
-
int d;
mpi_out_data *out_bufs;
MPI_Comm comm;
@@ -450,7 +452,6 @@ real sparse_approx_inverse(reax_system *system, simulation_data *data, storage *
const neighbor_proc *nbr1, *nbr2;
int *j_send, *j_recv1, *j_recv2;
real *val_send, *val_recv1, *val_recv2;
-
real start, t_start, t_comm;
real total_comm;
@@ -461,13 +462,15 @@ real sparse_approx_inverse(reax_system *system, simulation_data *data, storage *
if ( *A_app_inv == NULL)
{
- Allocate_Matrix2( A_app_inv, A_spar_patt->n, system->local_cap, A_spar_patt->m, comm );
+ Allocate_Matrix2( A_app_inv, A_spar_patt->n, system->local_cap, A_spar_patt->m,
+ SYM_FULL_MATRIX, comm );
}
else /* if ( (*A_app_inv)->m < A_spar_patt->m ) */
{
Deallocate_Matrix( *A_app_inv );
- Allocate_Matrix2( A_app_inv, A_spar_patt->n, system->local_cap, A_spar_patt->m, comm );
+ Allocate_Matrix2( A_app_inv, A_spar_patt->n, system->local_cap, A_spar_patt->m,
+ SYM_FULL_MATRIX, comm );
}
pos_x = NULL;
@@ -899,20 +902,17 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
b[i][0] = b[i][1] = 0;
}
- /* perform multiplication */
- for ( i = 0; i < A->n; ++i )
+ if ( A->format == SYM_HALF_MATRIX )
{
- si = A->start[i];
-#if defined(HALF_LIST)
- b[i][0] += A->entries[si].val * x[i][0];
- b[i][1] += A->entries[si].val * x[i][1];
-#endif
+ /* perform multiplication */
+ for ( i = 0; i < A->n; ++i )
+ {
+ si = A->start[i];
-#if defined(HALF_LIST)
- for ( k = si + 1; k < A->end[i]; ++k )
-#else
- for ( k = si; k < A->end[i]; ++k )
-#endif
+ b[i][0] += A->entries[si].val * x[i][0];
+ b[i][1] += A->entries[si].val * x[i][1];
+
+ for ( k = si + 1; k < A->end[i]; ++k )
{
j = A->entries[k].j;
H = A->entries[k].val;
@@ -920,14 +920,30 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
b[i][0] += H * x[j][0];
b[i][1] += H * x[j][1];
-#if defined(HALF_LIST)
// comment out for tryQEq
//if( j < A->n ) {
b[j][0] += H * x[i][0];
b[j][1] += H * x[i][1];
//}
-#endif
}
+ }
+ }
+ else if ( A->format == SYM_FULL_MATRIX )
+ {
+ /* perform multiplication */
+ for ( i = 0; i < A->n; ++i )
+ {
+ si = A->start[i];
+
+ for ( k = si; k < A->end[i]; ++k )
+ {
+ j = A->entries[k].j;
+ H = A->entries[k].val;
+
+ b[i][0] += H * x[j][0];
+ b[i][1] += H * x[j][1];
+ }
+ }
}
}
@@ -958,10 +974,10 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
Dist( system, mpi_data, x, mpi_data->mpi_rvec2, scale, rvec2_packer );
dual_Sparse_MatVec( H, x, workspace->q2, N );
-#if defined(HALF_LIST)
- // tryQEq
- Coll(system, mpi_data, workspace->q2, mpi_data->mpi_rvec2, scale, rvec2_unpacker);
-#endif
+ if ( H->format == SYM_HALF_MATRIX )
+ {
+ Coll(system, mpi_data, workspace->q2, mpi_data->mpi_rvec2, scale, rvec2_unpacker);
+ }
#if defined(CG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
@@ -1011,10 +1027,10 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
{
Dist(system, mpi_data, workspace->d2, mpi_data->mpi_rvec2, scale, rvec2_packer);
dual_Sparse_MatVec( H, workspace->d2, workspace->q2, N );
-#if defined(HALF_LIST)
- // tryQEq
- Coll(system, mpi_data, workspace->q2, mpi_data->mpi_rvec2, scale, rvec2_unpacker);
-#endif
+ if ( H->format == SYM_HALF_MATRIX )
+ {
+ Coll(system, mpi_data, workspace->q2, mpi_data->mpi_rvec2, scale, rvec2_unpacker);
+ }
#if defined(CG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
@@ -1124,66 +1140,50 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
{
- int i, j, k, si;
- real H;
+ int i, j, k, si;
+ real val;
for ( i = 0; i < N; ++i )
{
b[i] = 0;
}
- /* perform multiplication */
- for ( i = 0; i < A->n; ++i )
+ if ( A->format == SYM_HALF_MATRIX )
{
- si = A->start[i];
+ for ( i = 0; i < A->n; ++i )
+ {
+ si = A->start[i];
-#if defined(HALF_LIST)
- b[i] += A->entries[si].val * x[i];
-#endif
+ /* diagonal only contributes once */
+ b[i] += A->entries[si].val * x[i];
-#if defined(HALF_LIST)
- for ( k = si + 1; k < A->end[i]; ++k )
-#else
- for ( k = si; k < A->end[i]; ++k )
-#endif
+ for ( k = si + 1; k < A->end[i]; ++k )
{
j = A->entries[k].j;
- H = A->entries[k].val;
+ val = A->entries[k].val;
- b[i] += H * x[j];
-#if defined(HALF_LIST)
+ b[i] += val * x[j];
//if( j < A->n ) // comment out for tryQEq
- b[j] += H * x[i];
-#endif
+ b[j] += val * x[i];
}
+ }
}
-}
-
-
-/* sparse matrix-vector product Ax = b
- * where:
- * A: matrix, stored in CSR format
- * x: vector
- * b: vector (result)
-static void Sparse_MatVec_full( const sparse_matrix * const A,
- const real * const x, real * const b )
-{
- TODO: implement full SpMV in MPI
- int i, pj;
-
- Vector_MakeZero( b, A->n );
-
- #ifdef _OPENMP
- #pragma omp for schedule(static)
- #endif
+ else if ( A->format == SYM_FULL_MATRIX )
+ {
for ( i = 0; i < A->n; ++i )
{
- for ( pj = A->start[i]; pj < A->start[i + 1]; ++pj )
+ si = A->start[i];
+
+ for ( k = si; k < A->end[i]; ++k )
{
- b[i] += A->val[pj] * x[A->j[pj]];
+ j = A->entries[k].j;
+ val = A->entries[k].val;
+
+ b[i] += val * x[j];
}
}
-}*/
+ }
+}
int CG( reax_system *system, control_params *control, simulation_data *data,
@@ -1211,11 +1211,12 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
Sparse_MatVec( H, x, workspace->q, system->N );
t_spmv += Get_Timing_Info( t_start );
-#if defined(HALF_LIST)
- t_start = Get_Time( );
- Coll( system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker );
- t_comm += Get_Timing_Info( t_start );
-#endif
+ if ( H->format == SYM_HALF_MATRIX )
+ {
+ t_start = Get_Time( );
+ Coll( system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker );
+ t_comm += Get_Timing_Info( t_start );
+ }
t_start = Get_Time( );
Vector_Sum( workspace->r , 1., b, -1., workspace->q, system->n );
@@ -1258,11 +1259,12 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
Sparse_MatVec( H, workspace->d, workspace->q, system->N );
t_spmv += Get_Timing_Info( t_start );
-#if defined(HALF_LIST)
- t_start = Get_Time( );
- Coll(system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker);
- t_comm += Get_Timing_Info( t_start );
-#endif
+ if ( H->format == SYM_HALF_MATRIX )
+ {
+ t_start = Get_Time( );
+ Coll(system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker);
+ t_comm += Get_Timing_Info( t_start );
+ }
t_start = Get_Time( );
tmp = Parallel_Dot(workspace->d, workspace->q, system->n, mpi_data->world);
diff --git a/PuReMD/src/list.c b/PuReMD/src/list.c
index 4ccb03ed..5bfba249 100644
--- a/PuReMD/src/list.c
+++ b/PuReMD/src/list.c
@@ -30,7 +30,8 @@
/************* allocate list space ******************/
-int Make_List(int n, int num_intrs, int type, reax_list *l, MPI_Comm comm)
+int Make_List( int n, int num_intrs, int type, int format,
+ reax_list *l, MPI_Comm comm )
{
l->allocated = 1;
@@ -41,6 +42,7 @@ int Make_List(int n, int num_intrs, int type, reax_list *l, MPI_Comm comm)
l->end_index = (int*) smalloc( n * sizeof(int), "list:end_index", comm );
l->type = type;
+ l->format = format;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "list: n=%d num_intrs=%d type=%d\n", l->n, l->num_intrs, l->type );
diff --git a/PuReMD/src/list.h b/PuReMD/src/list.h
index da400b76..918256f5 100644
--- a/PuReMD/src/list.h
+++ b/PuReMD/src/list.h
@@ -25,7 +25,7 @@
#include "reax_types.h"
-int Make_List( int, int, int, reax_list*, MPI_Comm );
+int Make_List( int, int, int, int, reax_list*, MPI_Comm );
void Delete_List( reax_list*, MPI_Comm );
diff --git a/PuReMD/src/neighbors.c b/PuReMD/src/neighbors.c
index fa86a738..96800433 100644
--- a/PuReMD/src/neighbors.c
+++ b/PuReMD/src/neighbors.c
@@ -110,23 +110,18 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
itr = 0;
while ( (gcj = gci->nbrs[itr]) != NULL )
{
- if (
-#if defined(HALF_LIST)
- gci->str <= gcj->str &&
-#endif
- (DistSqr_to_Special_Point(gci->nbrs_cp[itr], atom1->x) <= cutoff) )
+ if ( ((far_nbrs->format == HALF_LIST && gci->str <= gcj->str)
+ || far_nbrs->format == FULL_LIST)
+ && (DistSqr_to_Special_Point(gci->nbrs_cp[itr], atom1->x) <= cutoff) )
{
/* pick up another atom from the neighbor cell */
for ( m = gcj->str; m < gcj->end; ++m )
{
-#if defined(HALF_LIST)
- /* prevent recounting same pairs within a gcell and
- * make half-list */
- if ( l < m )
-#else
- /* prevent recounting same pairs within a gcell */
- if ( l != m )
-#endif
+ /* HALF_LIST: prevent recounting same pairs within a gcell and
+ * make half-list
+ * FULL_LIST: prevent recounting same pairs within a gcell */
+ if ( (far_nbrs->format == HALF_LIST && l < m)
+ || (far_nbrs->format == FULL_LIST && l != m) )
{
atom2 = &(system->my_atoms[m]);
dvec[0] = atom2->x[0] - atom1->x[0];
@@ -181,7 +176,8 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
}
-int Estimate_NumNeighbors( reax_system *system, reax_list **lists )
+int Estimate_NumNeighbors( reax_system *system, reax_list **lists,
+ int far_nbr_list_format )
{
int i, j, k, l, m, itr, num_far; //, tmp, tested;
real d, cutoff;
@@ -214,23 +210,18 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists )
itr = 0;
while ( (gcj = gci->nbrs[itr]) != NULL )
{
- if (
-#if defined(HALF_LIST)
- gci->str <= gcj->str &&
-#endif
- (DistSqr_to_Special_Point(gci->nbrs_cp[itr], atom1->x) <= cutoff))
+ if ( ((far_nbr_list_format == HALF_LIST && gci->str <= gcj->str)
+ || far_nbr_list_format == FULL_LIST)
+ && (DistSqr_to_Special_Point(gci->nbrs_cp[itr], atom1->x) <= cutoff))
{
/* pick up another atom from the neighbor cell */
for ( m = gcj->str; m < gcj->end; ++m )
{
-#if defined(HALF_LIST)
- /* prevent recounting same pairs within a gcell and
- * make half-list */
- if ( l < m )
-#else
- /* prevent recounting same pairs within a gcell */
- if ( l != m )
-#endif
+ /* HALF_LIST: prevent recounting same pairs within a gcell and
+ * make half-list
+ * FULL_LIST: prevent recounting same pairs within a gcell */
+ if ( (far_nbr_list_format == HALF_LIST && l < m)
+ || (far_nbr_list_format == FULL_LIST && l != m) )
{
//fprintf( stderr, "\t\t\tatom2=%d\n", m );
atom2 = &(system->my_atoms[m]);
diff --git a/PuReMD/src/neighbors.h b/PuReMD/src/neighbors.h
index 0a1e3daf..818fc724 100644
--- a/PuReMD/src/neighbors.h
+++ b/PuReMD/src/neighbors.h
@@ -33,6 +33,6 @@
void Generate_Neighbor_Lists( reax_system*, simulation_data*, storage*,
reax_list** );
-int Estimate_NumNeighbors( reax_system*, reax_list** );
+int Estimate_NumNeighbors( reax_system*, reax_list**, int );
#endif
diff --git a/PuReMD/src/nonbonded.c b/PuReMD/src/nonbonded.c
index bb593b72..4f07103f 100644
--- a/PuReMD/src/nonbonded.c
+++ b/PuReMD/src/nonbonded.c
@@ -71,11 +71,9 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
j = nbr_pj->nbr;
orig_j = system->my_atoms[j].orig_id;
-#if defined(HALF_LIST)
- if ( nbr_pj->d <= control->nonb_cut && (j < natoms || orig_i < orig_j) )
-#else
- if ( nbr_pj->d <= control->nonb_cut && orig_i < orig_j )
-#endif
+ if ( nbr_pj->d <= control->nonb_cut
+ && ((far_nbrs->format == HALF_LIST && (j < natoms || orig_i < orig_j))
+ || (far_nbrs->format == FULL_LIST && orig_i < orig_j)) )
{
r_ij = nbr_pj->d;
twbp = &(system->reax_param.tbp[ system->my_atoms[i].type ]
@@ -253,11 +251,9 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
j = nbr_pj->nbr;
orig_j = system->my_atoms[j].orig_id;
-#if defined(HALF_LIST)
- if ( nbr_pj->d <= control->nonb_cut && (j < natoms || orig_i < orig_j) )
-#else
- if ( nbr_pj->d <= control->nonb_cut && orig_i < orig_j )
-#endif
+ if ( nbr_pj->d <= control->nonb_cut
+ && ((far_nbrs->format == HALF_LIST && (j < natoms || orig_i < orig_j))
+ || (far_nbrs->format == FULL_LIST && orig_i < orig_j)) )
{
j = nbr_pj->nbr;
type_j = system->my_atoms[j].type;
diff --git a/PuReMD/src/qEq.c b/PuReMD/src/qEq.c
index 3ec05786..15680d87 100644
--- a/PuReMD/src/qEq.c
+++ b/PuReMD/src/qEq.c
@@ -247,8 +247,8 @@ void Init_MatVec( reax_system *system, simulation_data *data,
if( workspace->L == NULL ) {
fillin = Estimate_LU_Fill( workspace->H, workspace->droptol );
- if( Allocate_Matrix( &(workspace->L), workspace->H->cap, fillin ) == 0 ||
- Allocate_Matrix( &(workspace->U), workspace->H->cap, fillin ) == 0 ) {
+ if( Allocate_Matrix( &(workspace->L), workspace->H->cap, fillin, FULL_MATRIX, comm ) == 0 ||
+ Allocate_Matrix( &(workspace->U), workspace->H->cap, fillin, FULL_MATRIX, comm ) == 0 ) {
fprintf( stderr, "not enough memory for LU matrices. terminating.\n" );
MPI_Abort( mpi_data->world, INSUFFICIENT_MEMORY );
}
diff --git a/PuReMD/src/reax_defs.h b/PuReMD/src/reax_defs.h
index a61ce245..dc2ac1b9 100644
--- a/PuReMD/src/reax_defs.h
+++ b/PuReMD/src/reax_defs.h
@@ -26,10 +26,10 @@
#define inline __inline__
#endif /*IBMC*/
-#define SUCCESS 1
-#define FAILURE 0
-#define TRUE 1
-#define FALSE 0
+#define SUCCESS (1)
+#define FAILURE (0)
+#define TRUE (1)
+#define FALSE (0)
#define SQR(x) ((x)*(x))
#define CUBE(x) ((x)*(x)*(x))
@@ -39,65 +39,65 @@
#define MIN(x,y) (((x) < (y)) ? (x) : (y))
#define MAX3(x,y,z) MAX( MAX(x,y), z)
-#define constPI 3.14159265
-#define C_ele 332.06371
-//#define K_B 503.398008 // kcal/mol/K
-#define K_B 0.831687 // amu A^2 / ps^2 / K
-#define F_CONV 1e6 / 48.88821291 / 48.88821291 // --> amu A / ps^2
-#define E_CONV 0.002391 // amu A^2 / ps^2 --> kcal/mol
-#define EV_to_KCALpMOL 14.400000 // ElectronVolt --> KCAL per MOLe
-#define KCALpMOL_to_EV 23.02 // 23.060549 //KCAL per MOLe --> ElectronVolt
-#define ECxA_to_DEBYE 4.803204 // elem. charge * Ang -> debye
-#define CAL_to_JOULES 4.184000 // CALories --> JOULES
-#define JOULES_to_CAL 1/4.184000 // JOULES --> CALories
-#define AMU_to_GRAM 1.6605e-24
-#define ANG_to_CM 1e-8
-#define AVOGNR 6.0221367e23
-#define P_CONV 1e-24 * AVOGNR * JOULES_to_CAL
-
-#define MAX_STR 1024
-#define MAX_LINE 1024
-#define MAX_TOKENS 1024
-#define MAX_TOKEN_LEN 1024
-
-#define MAX_ATOM_ID 100000
-#define MAX_RESTRICT 15
-#define MAX_MOLECULE_SIZE 20
-#define MAX_ATOM_TYPES 25
-
-#define NUM_INTRS 10
-#define ALMOST_ZERO 1e-10
-#define NEG_INF -1e10
-#define NO_BOND 1e-3 // 0.001
-#define HB_THRESHOLD 1e-2 // 0.01
-
-#define MIN_CAP 50
-#define MIN_NBRS 100
-#define MIN_HENTRIES 100
-#define MAX_BONDS 30
-#define MIN_BONDS 15
-#define MIN_HBONDS 25
-#define MIN_3BODIES 1000
-#define MIN_GCELL_POPL 50
-#define MIN_SEND 100
-#define SAFE_ZONE 1.2
-#define SAFER_ZONE 1.4
-#define DANGER_ZONE 0.90
-#define LOOSE_ZONE 0.75
-#define MAX_3BODY_PARAM 5
-#define MAX_4BODY_PARAM 5
-
-#define MAX_dV 1.01
-#define MIN_dV 0.99
-#define MAX_dT 4.00
-#define MIN_dT 0.00
-
-#define MASTER_NODE 0
-#define MAX_NBRS 6 //27
-#define MYSELF 13 // encoding of relative coordinate (0,0,0)
-
-#define MAX_ITR 10
-#define RESTART 30
+#define constPI (3.14159265)
+#define C_ele (332.06371)
+//#define K_B (503.398008) // kcal/mol/K
+#define K_B (0.831687) // amu A^2 / ps^2 / K
+#define F_CONV (1e6 / 48.88821291 / 48.88821291) // --> amu A / ps^2
+#define E_CONV (0.002391) // amu A^2 / ps^2 --> kcal/mol
+#define EV_to_KCALpMOL (14.400000) // ElectronVolt --> KCAL per MOLe
+#define KCALpMOL_to_EV (23.02) // 23.060549 //KCAL per MOLe --> ElectronVolt
+#define ECxA_to_DEBYE (4.803204) // elem. charge * Ang -> debye
+#define CAL_to_JOULES (4.184000) // CALories --> JOULES
+#define JOULES_to_CAL (1/4.184000) // JOULES --> CALories
+#define AMU_to_GRAM (1.6605e-24)
+#define ANG_to_CM (1e-8)
+#define AVOGNR (6.0221367e23)
+#define P_CONV (1e-24 * AVOGNR * JOULES_to_CAL)
+
+#define MAX_STR (1024)
+#define MAX_LINE (1024)
+#define MAX_TOKENS (1024)
+#define MAX_TOKEN_LEN (1024)
+
+#define MAX_ATOM_ID (100000)
+#define MAX_RESTRICT (15)
+#define MAX_MOLECULE_SIZE (20)
+#define MAX_ATOM_TYPES (25)
+
+#define NUM_INTRS (10)
+#define ALMOST_ZERO (1e-10)
+#define NEG_INF (-1e10)
+#define NO_BOND (1e-3) // 0.001
+#define HB_THRESHOLD (1e-2) // 0.01
+
+#define MIN_CAP (50)
+#define MIN_NBRS (100)
+#define MIN_HENTRIES (100)
+#define MAX_BONDS (30)
+#define MIN_BONDS (15)
+#define MIN_HBONDS (25)
+#define MIN_3BODIES (1000)
+#define MIN_GCELL_POPL (50)
+#define MIN_SEND (100)
+#define SAFE_ZONE (1.2)
+#define SAFER_ZONE (1.4)
+#define DANGER_ZONE (0.90)
+#define LOOSE_ZONE (0.75)
+#define MAX_3BODY_PARAM (5)
+#define MAX_4BODY_PARAM (5)
+
+#define MAX_dV (1.01)
+#define MIN_dV (0.99)
+#define MAX_dT (4.00)
+#define MIN_dT (0.00)
+
+#define MASTER_NODE (0)
+#define MAX_NBRS (6) //27
+#define MYSELF (13) // encoding of relative coordinate (0,0,0)
+
+#define MAX_ITR (10)
+#define RESTART (30)
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 56a8c1c9..25749323 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -41,7 +41,6 @@
#define PURE_REAX
//#define LAMMPS_REAX
-//#define HALF_LIST
//#define DEBUG
//#define DEBUG_FOCUS
//#define TEST_ENERGY
@@ -110,6 +109,26 @@ enum pre_app
JACOBI_ITER_PA = 3,
};
+/* interaction list (reax_list) storage format */
+enum reax_list_format
+{
+ /* store half of interactions, when i < j (atoms i and j) */
+ HALF_LIST = 0,
+ /* store all interactions */
+ FULL_LIST = 1,
+};
+
+/* sparse matrix (sparse_matrix) storage format */
+enum sparse_matrix_format
+{
+ /* store upper half of nonzeros in a symmetric matrix (a_{ij}, i >= j) */
+ SYM_HALF_MATRIX = 0,
+ /* store all nonzeros in a symmetric matrix */
+ SYM_FULL_MATRIX = 1,
+ /* store all nonzeros in a matrix */
+ FULL_MATRIX = 2,
+};
+
/********************** TYPE DEFINITIONS ********************/
typedef int ivec[3];
@@ -935,6 +954,8 @@ typedef struct
typedef struct
{
+ /* matrix storage format */
+ int format;
int cap, n, m;
int *start, *end;
sparse_matrix_entry *entries;
@@ -1066,6 +1087,9 @@ typedef struct
int type;
// list_type select;
+ /* list storage format (half or full) */
+ int format;
+
void *v;
three_body_interaction_data *three_body_list;
bond_data *bond_list;
--
GitLab