From 83e5e0c1ca9885f38654418388bf7e3fafe25f28 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Sat, 18 Mar 2017 15:55:45 -0400
Subject: [PATCH] PG-PuReMD: fix issue with hydrogen bond cutoff with systems
without hydrogen atoms. Remove debugging info for MPI-only version of
PG-PuReMD.
---
PG-PuReMD/src/cuda_forces.cu | 62 +++++++++++++++++++++-------------
PG-PuReMD/src/forces.c | 6 ++++
PG-PuReMD/src/init_md.c | 47 +++++++++++++++++---------
PG-PuReMD/src/linear_solvers.c | 8 +++++
PG-PuReMD/src/neighbors.c | 15 ++++++--
PG-PuReMD/src/qEq.c | 8 +++++
6 files changed, 104 insertions(+), 42 deletions(-)
diff --git a/PG-PuReMD/src/cuda_forces.cu b/PG-PuReMD/src/cuda_forces.cu
index cec6c0eb..acc3fa5b 100644
--- a/PG-PuReMD/src/cuda_forces.cu
+++ b/PG-PuReMD/src/cuda_forces.cu
@@ -30,7 +30,14 @@ extern "C" int Make_List( int, int, int, reax_list*);
extern "C" void Delete_List( reax_list*);
-CUDA_GLOBAL void ker_estimate_storages (reax_atom *my_atoms,
+
+CUDA_GLOBAL void k_disable_hydrogen_bonding( control_params *control )
+{
+ control->hbond_cut = 0.0;
+}
+
+
+CUDA_GLOBAL void k_estimate_storages (reax_atom *my_atoms,
single_body_parameters *sbp,
two_body_parameters *tbp,
control_params *control,
@@ -105,7 +112,7 @@ CUDA_GLOBAL void ker_estimate_storages (reax_atom *my_atoms,
&& (i > n)
)
{
- atomicAdd (&hb_top [i], 1);
+ atomicAdd( &hb_top[i], 1 );
}
if (i >= n)
@@ -145,14 +152,14 @@ CUDA_GLOBAL void ker_estimate_storages (reax_atom *my_atoms,
if( (ihb == 1) && (jhb == 2))
{
//++hb_top[i];
- atomicAdd (&hb_top[i], 1);
+ atomicAdd( &hb_top[i], 1 );
}
//else if( j < n && ihb == 2 && jhb == 1 )
//else if( ihb == 2 && jhb == 1 && j < n)
else if( ihb == 2 && jhb == 1 && j < n)
{
//++hb_top[j];
- atomicAdd (&hb_top[i], 1);
+ atomicAdd( &hb_top[i], 1 );
}
}
}
@@ -195,7 +202,7 @@ CUDA_GLOBAL void ker_estimate_storages (reax_atom *my_atoms,
}
-CUDA_GLOBAL void ker_init_system_atoms(reax_atom *my_atoms, int N,
+CUDA_GLOBAL void k_init_system_atoms(reax_atom *my_atoms, int N,
int *hb_top, int *bond_top)
{
int i;
@@ -242,7 +249,7 @@ void Cuda_Estimate_Storages(reax_system *system, control_params *control,
blocks = (int) CEIL((real)system->N / ST_BLOCK_SIZE);
- ker_estimate_storages <<< blocks, ST_BLOCK_SIZE>>>
+ k_estimate_storages <<< blocks, ST_BLOCK_SIZE>>>
(system->d_my_atoms, system->reax_param.d_sbp, system->reax_param.d_tbp,
(control_params *)control->d_control_params, *(*dev_lists + FAR_NBRS), system->reax_param.num_atom_types,
system->n, system->N, system->Hcap, system->total_cap,
@@ -284,7 +291,7 @@ void Cuda_Estimate_Storages(reax_system *system, control_params *control,
min_hbonds = hb_top[i];
}
- hbond_count += hb_top [i];
+ hbond_count += hb_top[i];
}
system->max_hbonds = max_hbonds * SAFER_ZONE;
@@ -296,7 +303,14 @@ void Cuda_Estimate_Storages(reax_system *system, control_params *control,
system->my_rank, min_bonds, max_bonds, min_hbonds, max_hbonds);
#endif
- ker_init_system_atoms <<<blocks, ST_BLOCK_SIZE>>>
+ // if number of hydrogen atoms is 0, disable hydrogen bond functionality
+ if ( hbond_count == 0 )
+ {
+ control->hbond_cut = 0.0;
+ k_disable_hydrogen_bonding <<<1,1>>> ( (control_params *)control->d_control_params );
+ }
+
+ k_init_system_atoms <<<blocks, ST_BLOCK_SIZE>>>
(system->d_my_atoms, system->N, d_hb_top, d_bond_top );
cudaThreadSynchronize();
@@ -345,7 +359,7 @@ CUDA_DEVICE real Compute_tabH( LR_lookup_table *t_LR, real r_ij, int ti, int tj,
}
-CUDA_GLOBAL void ker_estimate_sparse_matrix (reax_atom *my_atoms, control_params *control,
+CUDA_GLOBAL void k_estimate_sparse_matrix (reax_atom *my_atoms, control_params *control,
reax_list p_far_nbrs, int n, int N, int renbr, int *indices)
{
int i, j, pj;
@@ -459,7 +473,7 @@ int Cuda_Estimate_Sparse_Matrix (reax_system *system, control_params *control,
//TODO
//TODO
//TODO
- ker_estimate_sparse_matrix <<< blocks, DEF_BLOCK_SIZE >>>
+ k_estimate_sparse_matrix <<< blocks, DEF_BLOCK_SIZE >>>
(system->d_my_atoms, (control_params *)control->d_control_params,
*(*dev_lists + FAR_NBRS), system->n, system->N,
(((data->step-data->prev_steps) % control->reneighbor) == 0), indices);
@@ -483,7 +497,7 @@ int Cuda_Estimate_Sparse_Matrix (reax_system *system, control_params *control,
}
-CUDA_GLOBAL void ker_init_forces (reax_atom *my_atoms, single_body_parameters *sbp,
+CUDA_GLOBAL void k_init_forces (reax_atom *my_atoms, single_body_parameters *sbp,
two_body_parameters *tbp, storage workspace,
control_params *control,
reax_list far_nbrs, reax_list bonds, reax_list hbonds,
@@ -555,7 +569,7 @@ CUDA_GLOBAL void ker_init_forces (reax_atom *my_atoms, single_body_parameters *s
}
//CHANGE ORIGINAL
- if( control->hbond_cut > 0 ) {
+ if( control->hbond_cut > 0.0 ) {
ihb = sbp_i->p_hbond;
//CHANGE ORIGINAL
if( ihb == 1 || ihb == 2) {
@@ -756,7 +770,7 @@ CUDA_GLOBAL void ker_init_forces (reax_atom *my_atoms, single_body_parameters *s
}
-CUDA_GLOBAL void ker_init_bond_mark (int offset, int n, int *bond_mark)
+CUDA_GLOBAL void k_init_bond_mark (int offset, int n, int *bond_mark)
{
int i;
@@ -848,7 +862,7 @@ CUDA_GLOBAL void New_fix_sym_hbond_indices (reax_atom *my_atoms, reax_list hbond
////////////////////////
// HBOND ISSUE
-CUDA_GLOBAL void ker_update_bonds (reax_atom *my_atoms,
+CUDA_GLOBAL void k_update_bonds (reax_atom *my_atoms,
reax_list bonds,
int n)
{
@@ -860,7 +874,7 @@ CUDA_GLOBAL void ker_update_bonds (reax_atom *my_atoms,
}
-CUDA_GLOBAL void ker_update_hbonds (reax_atom *my_atoms,
+CUDA_GLOBAL void k_update_hbonds (reax_atom *my_atoms,
reax_list hbonds,
int n)
{
@@ -894,7 +908,7 @@ int Cuda_Validate_Lists (reax_system *system, storage *workspace, reax_list **li
blocks = system->n / DEF_BLOCK_SIZE +
((system->n % DEF_BLOCK_SIZE == 0) ? 0 : 1);
- ker_update_bonds <<< blocks, DEF_BLOCK_SIZE >>>
+ k_update_bonds <<< blocks, DEF_BLOCK_SIZE >>>
(system->d_my_atoms, *(*lists + BONDS),
system->n);
cudaThreadSynchronize ();
@@ -904,7 +918,7 @@ int Cuda_Validate_Lists (reax_system *system, storage *workspace, reax_list **li
// HBOND ISSUE
//FIX - 4 - Added this check for hydrogen bond issue
if ((control->hbond_cut > 0) && (system->numH > 0)){
- ker_update_hbonds <<< blocks, DEF_BLOCK_SIZE >>>
+ k_update_hbonds <<< blocks, DEF_BLOCK_SIZE >>>
(system->d_my_atoms, *(*lists + HBONDS),
system->n);
cudaThreadSynchronize ();
@@ -1084,7 +1098,7 @@ int Cuda_Validate_Lists (reax_system *system, storage *workspace, reax_list **li
}
-CUDA_GLOBAL void ker_init_bond_orders (reax_atom *my_atoms,
+CUDA_GLOBAL void k_init_bond_orders (reax_atom *my_atoms,
reax_list far_nbrs,
reax_list bonds,
real *total_bond_order,
@@ -1114,7 +1128,7 @@ CUDA_GLOBAL void ker_init_bond_orders (reax_atom *my_atoms,
}
-CUDA_GLOBAL void ker_bond_mark (reax_list p_bonds, storage p_workspace, int N)
+CUDA_GLOBAL void k_bond_mark (reax_list p_bonds, storage p_workspace, int N)
{
reax_list *bonds = &( p_bonds );
storage *workspace = &( p_workspace );
@@ -1153,7 +1167,7 @@ int Cuda_Init_Forces( reax_system *system, control_params *control,
blocks = (system->N - system->n) / DEF_BLOCK_SIZE +
(((system->N - system->n) % DEF_BLOCK_SIZE == 0) ? 0 : 1);
- ker_init_bond_mark <<< blocks, DEF_BLOCK_SIZE >>>
+ k_init_bond_mark <<< blocks, DEF_BLOCK_SIZE >>>
(system->n, (system->N - system->n), dev_workspace->bond_mark);
cudaThreadSynchronize ();
cudaCheckError ();
@@ -1165,14 +1179,14 @@ int Cuda_Init_Forces( reax_system *system, control_params *control,
(((system->N % DEF_BLOCK_SIZE) == 0) ? 0 : 1);
//fprintf (stderr, " Total atoms: %d, blocks: %d \n", system->N, init_blocks );
- // ker_init_bond_orders <<<init_blocks, DEF_BLOCK_SIZE >>>
+ // k_init_bond_orders <<<init_blocks, DEF_BLOCK_SIZE >>>
// ( system->d_my_atoms, *(*dev_lists + FAR_NBRS), *(*dev_lists + BONDS),
// dev_workspace->total_bond_order, system->N);
// cudaThreadSynchronize ();
// cudaCheckError ();
// fprintf (stderr, " DONE WITH VALIDATION \n");
- ker_init_forces <<<init_blocks, DEF_BLOCK_SIZE >>>
+ k_init_forces <<<init_blocks, DEF_BLOCK_SIZE >>>
(system->d_my_atoms, system->reax_param.d_sbp,
system->reax_param.d_tbp, *dev_workspace,
(control_params *)control->d_control_params,
@@ -1207,9 +1221,9 @@ int Cuda_Init_Forces( reax_system *system, control_params *control,
}
//update bond_mark
- //ker_bond_mark <<< init_blocks, DEF_BLOCK_SIZE>>>
+ //k_bond_mark <<< init_blocks, DEF_BLOCK_SIZE>>>
/*
- ker_bond_mark <<< 1, 1>>>
+ k_bond_mark <<< 1, 1>>>
( *(*dev_lists + BONDS), *dev_workspace, system->N);
cudaThreadSynchronize ();
cudaCheckError ();
diff --git a/PG-PuReMD/src/forces.c b/PG-PuReMD/src/forces.c
index 52092aac..ad98d3c3 100644
--- a/PG-PuReMD/src/forces.c
+++ b/PG-PuReMD/src/forces.c
@@ -378,8 +378,11 @@ int MPI_Not_GPU_Validate_Lists (reax_system *system, storage *workspace, reax_li
//update the current step max_sp_entries;
realloc->Htop = max_sp_entries;
+
+#if defined(DEBUG)
fprintf (stderr, "p:%d - MPI-Not-GPU Reallocate: Total H matrix entries: %d, cap: %d, used: %d \n",
system->my_rank, workspace->H.n, workspace->H.m, total_sp_entries);
+#endif
if (total_sp_entries >= workspace->H.m)
{
@@ -443,8 +446,11 @@ int MPI_Not_GPU_Validate_Lists (reax_system *system, storage *workspace, reax_li
{
if (end_index[i] - index[i] >= system->max_bonds)
{
+#if defined(DEBUG)
fprintf( stderr, "MPI-Not-GPU step%d-bondchk failed: i=%d start(i)=%d end(i)=%d max_bonds=%d\n",
step, i, index[i], end_index[i], system->max_bonds);
+#endif
+
return FAILURE;
}
if (end_index[i] - index[i] >= max_bonds)
diff --git a/PG-PuReMD/src/init_md.c b/PG-PuReMD/src/init_md.c
index f2df3ef7..80498ce8 100644
--- a/PG-PuReMD/src/init_md.c
+++ b/PG-PuReMD/src/init_md.c
@@ -154,24 +154,34 @@ int Init_System( reax_system *system, control_params *control,
Reorder_My_Atoms( system, workspace );
/* estimate N and total capacity */
- for ( i = 0; i < MAX_NBRS; ++i ) nrecv[i] = 0;
+ for ( i = 0; i < MAX_NBRS; ++i )
+ {
+ nrecv[i] = 0;
+ }
MPI_Barrier( MPI_COMM_WORLD );
system->N = SendRecv( system, mpi_data, mpi_data->boundary_atom_type, nrecv,
- Estimate_Boundary_Atoms, Unpack_Estimate_Message, 1 );
+ Estimate_Boundary_Atoms, Unpack_Estimate_Message, 1 );
system->total_cap = MAX( (int)(system->N * SAFE_ZONE), MIN_CAP );
Bin_Boundary_Atoms( system );
/* estimate numH and Hcap */
-
system->numH = 0;
- if ( control->hbond_cut > 0 )
+ if ( control->hbond_cut > 0.0 )
+ {
for ( i = 0; i < system->n; ++i )
{
atom = &(system->my_atoms[i]);
+
if ( system->reax_param.sbp[ atom->type ].p_hbond == 1 )
+ {
atom->Hindex = system->numH++;
- else atom->Hindex = -1;
+ }
+ else
+ {
+ atom->Hindex = -1;
+ }
}
+ }
//Tried fix
//system->Hcap = MAX( system->numH * SAFER_ZONE, MIN_CAP );
system->Hcap = MAX( system->n * SAFER_ZONE, MIN_CAP );
@@ -179,7 +189,7 @@ int Init_System( reax_system *system, control_params *control,
// Sudhir-style below
/*
system->numH = 0;
- if ( control->hbond_cut > 0 )
+ if ( control->hbond_cut > 0.0 )
for ( i = 0; i < system->n; ++i )
{
atom = &(system->my_atoms[i]);
@@ -253,7 +263,7 @@ int Cuda_Init_System( reax_system *system, control_params *control,
/* estimate numH and Hcap */
system->numH = 0;
- if ( control->hbond_cut > 0 )
+ if ( control->hbond_cut > 0.0 )
//TODO
//for( i = 0; i < system->n; ++i ) {
for ( i = 0; i < system->N; ++i )
@@ -778,7 +788,7 @@ int Init_Lists( reax_system *system, control_params *control,
(int)(Htop * sizeof(sparse_matrix_entry) / (1024 * 1024)) );
#endif
- if ( control->hbond_cut > 0 )
+ if ( control->hbond_cut > 0.0 )
{
// init H indexes
total_hbonds = 0;
@@ -880,7 +890,6 @@ int Cuda_Init_Lists( reax_system *system, control_params *control,
int total_hbonds, total_bonds, bond_cap, num_3body, cap_3body, Htop;
int *hb_top, *bond_top;
int nrecv[MAX_NBRS];
-
int *nbr_indices = (int *) host_scratch;
//num_nbrs = Estimate_NumNeighbors( system, lists );
@@ -890,19 +899,23 @@ int Cuda_Init_Lists( reax_system *system, control_params *control,
//fprintf (stderr, "atom: %d -- %d \n", i, nbr_indices[i]);
for (i = 0; i < system->N; i++)
- num_nbrs += nbr_indices [i] ;
+ {
+ num_nbrs += nbr_indices [i];
+ }
//fprintf (stderr, "DEVICE Total Neighbors: %d (%d)\n", num_nbrs, (int)(num_nbrs*SAFE_ZONE));
for (i = 0; i < system->N; i++)
+ {
nbr_indices[i] = MAX (nbr_indices [i] * SAFER_ZONE, MIN_NBRS);
+ }
num_nbrs = 0;
- num_nbrs += nbr_indices [0] ;
+ num_nbrs += nbr_indices[0];
for (i = 1; i < system->N; i++)
{
- num_nbrs += nbr_indices [i] ;
- nbr_indices [i] += nbr_indices [i - 1];
+ num_nbrs += nbr_indices[i];
+ nbr_indices[i] += nbr_indices[i - 1];
}
//fprintf (stderr, "DEVICE total neighbors entries: %d \n", nbr_indices [system->N - 1] );
@@ -929,7 +942,6 @@ int Cuda_Init_Lists( reax_system *system, control_params *control,
Cuda_Estimate_Storages( system, control, lists, system->local_cap, system->total_cap,
&Htop, hb_top, bond_top, &num_3body );
-
//TODO - CARVER FIX
//TODO - CARVER FIX
//TODO - CARVER FIX
@@ -965,7 +977,7 @@ int Cuda_Init_Lists( reax_system *system, control_params *control,
#endif
// FIX - 4 - Added addition check here for hydrogen Bonds
- if (( control->hbond_cut > 0 ) && (system->numH))
+ if (( control->hbond_cut > 0.0 ) && ( system->numH > 0 ))
{
/* init H indexes */
total_hbonds = 0;
@@ -978,7 +990,10 @@ int Cuda_Init_Lists( reax_system *system, control_params *control,
//TODO
hb_top [i] = MAX( hb_top[i] * 4, MIN_HBONDS * 4);
total_hbonds += hb_top[i];
- if (hb_top [i] > 0) count ++;
+ if (hb_top [i] > 0)
+ {
+ ++count;
+ }
}
total_hbonds = MAX( total_hbonds, MIN_CAP * MIN_HBONDS );
diff --git a/PG-PuReMD/src/linear_solvers.c b/PG-PuReMD/src/linear_solvers.c
index 3a8a5d3f..85e866ff 100644
--- a/PG-PuReMD/src/linear_solvers.c
+++ b/PG-PuReMD/src/linear_solvers.c
@@ -238,7 +238,11 @@ int dual_CG( reax_system *system, storage *workspace, sparse_matrix *H, rvec2
}
matvecs = CG( system, workspace, H, workspace->b_t, tol, workspace->t,
mpi_data, fout );
+
+#if defined(DEBUG)
fprintf (stderr, " CG1: iterations --> %d \n", matvecs );
+#endif
+
for ( j = 0; j < n; ++j )
{
workspace->x[j][1] = workspace->t[j];
@@ -252,7 +256,11 @@ int dual_CG( reax_system *system, storage *workspace, sparse_matrix *H, rvec2
}
matvecs = CG( system, workspace, H, workspace->b_s, tol, workspace->s,
mpi_data, fout );
+
+#if defined(DEBUG)
fprintf (stderr, " CG2: iterations --> %d \n", matvecs );
+#endif
+
for ( j = 0; j < system->n; ++j )
{
workspace->x[j][0] = workspace->s[j];
diff --git a/PG-PuReMD/src/neighbors.c b/PG-PuReMD/src/neighbors.c
index 01ed2482..ada23dc0 100644
--- a/PG-PuReMD/src/neighbors.c
+++ b/PG-PuReMD/src/neighbors.c
@@ -156,7 +156,9 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
}
}
+#if defined(DEBUG)
fprintf (stderr, " HOST NEIGHBOR COUNT: %d \n", num_far );
+#endif
workspace->realloc.num_far = num_far;
@@ -214,8 +216,15 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists )
for ( l = g->str[index_grid_3d (i, j, k, g)]; l < g->end[index_grid_3d (i, j, k, g)]; ++l )
{
atom1 = &(system->my_atoms[l]);
- if (l == 0) fprintf (stderr, "atom 0 has (%d %d %d) (%f %f %f) \n",
- i, j, k, atom1->x[0], atom1->x[1], atom1->x[2]);
+
+#if defined(DEBUG)
+ if (l == 0)
+ {
+ fprintf (stderr, "atom 0 has (%d %d %d) (%f %f %f) \n",
+ i, j, k, atom1->x[0], atom1->x[1], atom1->x[2]);
+ }
+#endif
+
//fprintf( stderr, "\tatom %d: ", l );
//tmp = num_far; tested = 0;
itr = 0;
@@ -251,7 +260,9 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists )
}
}
+#if defined(DEBUG)
fprintf (stderr, "Total number of host neighbors: %d \n", num_far);
+#endif
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: estimate nbrs done - num_far=%d\n",
diff --git a/PG-PuReMD/src/qEq.c b/PG-PuReMD/src/qEq.c
index 1e877887..578bf899 100644
--- a/PG-PuReMD/src/qEq.c
+++ b/PG-PuReMD/src/qEq.c
@@ -355,11 +355,19 @@ void Calculate_Charges( reax_system *system, storage *workspace,
my_sum[0] += workspace->x[i][0];
my_sum[1] += workspace->x[i][1];
}
+
+#if defined(DEBUG)
fprintf (stderr, "Host : my_sum[0]: %f and %f \n", my_sum[0], my_sum[1]);
+#endif
+
MPI_Allreduce( &my_sum, &all_sum, 2, MPI_DOUBLE, MPI_SUM, mpi_data->world );
u = all_sum[0] / all_sum[1];
+
+#if defined(DEBUG)
fprintf (stderr, "Host : u: %f \n", u);
+#endif
+
for ( i = 0; i < system->n; ++i )
{
atom = &( system->my_atoms[i] );
--
GitLab