diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index d2bf99717f4d0f4f9b5c1fda4db9026469699507..4935aaaebd878c2bc2a05fb848163d147745e2bf 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -219,6 +219,8 @@ void Reallocate_Part2( reax_system const * const system,
{
Reallocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max,
realloc->total_cm_entries );
+ Reallocate_Matrix( &workspace->H_sp, system->N_cm, system->N_cm_max,
+ realloc->total_cm_entries );
realloc->cm = FALSE;
}
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index d8c5b6aa22f8d48f69e6b3a025568692a006795c..df85a5e5cf6e154726a0d3deaaa3f59753ed4050 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -347,14 +347,16 @@ static inline real Init_Charge_Matrix_Entry( reax_system const * const system,
}
-static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const system,
+static int Init_Charge_Matrix_Remaining_Entries( reax_system const * const system,
control_params const * const control, reax_list const * const far_nbr_list,
sparse_matrix * const H, sparse_matrix * const H_sp,
int * const Htop, int * const H_sp_top )
{
- int i, j, pj, target, val_flag;
+ int i, j, pj, target, val_flag, flag_oom;
real d, xcut, bond_softness, * X_diag;
+ flag_oom = FALSE;
+
switch ( control->charge_method )
{
case QEQ_CM:
@@ -374,30 +376,48 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
{
H->j[*Htop] = i;
H->val[*Htop] = 1.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
H_sp->j[*H_sp_top] = i;
H_sp->val[*H_sp_top] = 1.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
}
- ++(*H_sp_top);
}
if ( *Htop < H->m )
{
H->j[*Htop] = system->N;
H->val[*Htop] = 0.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
H_sp->j[*H_sp_top] = system->N;
H_sp->val[*H_sp_top] = 0.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*H_sp_top);
}
else
{
@@ -414,15 +434,25 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
{
H->j[*Htop] = j - 1;
H->val[*Htop] = 1.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
H_sp->j[*H_sp_top] = j - 1;
H_sp->val[*H_sp_top] = 1.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
}
- ++(*H_sp_top);
}
/* explicit zeros on diagonals */
@@ -430,56 +460,85 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
{
H->j[*Htop] = system->N + i;
H->val[*Htop] = 0.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
H_sp->j[*H_sp_top] = system->N + i;
H_sp->val[*H_sp_top] = 0.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*H_sp_top);
}
for ( i = system->num_molec_charge_constraints;
i < system->num_molec_charge_constraints + system->num_custom_charge_constraints; ++i )
{
- H->start[system->N + i] = *Htop;
- H_sp->start[system->N + i] = *H_sp_top;
-
- for ( j = system->custom_charge_constraint_start[i - system->num_molec_charge_constraints];
- j <= system->custom_charge_constraint_start[i - system->num_molec_charge_constraints + 1]; ++j )
+ if ( flag_oom == FALSE )
{
- /* custom charge constraint on atoms */
+ H->start[system->N + i] = *Htop;
+ H_sp->start[system->N + i] = *H_sp_top;
+
+ for ( j = system->custom_charge_constraint_start[i - system->num_molec_charge_constraints];
+ j <= system->custom_charge_constraint_start[i - system->num_molec_charge_constraints + 1]; ++j )
+ {
+ /* custom charge constraint on atoms */
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->custom_charge_constraint_atom_index[j] - 1;
+ H->val[*Htop] = system->custom_charge_constraint_coeff[j];
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
+ }
+
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->custom_charge_constraint_atom_index[j] - 1;
+ H_sp->val[*H_sp_top] = system->custom_charge_constraint_coeff[j];
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
+ }
+ }
+
+ /* explicit zeros on diagonals */
if ( *Htop < H->m )
{
- H->j[*Htop] = system->custom_charge_constraint_atom_index[j] - 1;
- H->val[*Htop] = system->custom_charge_constraint_coeff[j];
+ H->j[*Htop] = system->N + i;
+ H->val[*Htop] = 0.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
- H_sp->j[*H_sp_top] = system->custom_charge_constraint_atom_index[j] - 1;
- H_sp->val[*H_sp_top] = system->custom_charge_constraint_coeff[j];
+ H_sp->j[*H_sp_top] = system->N + i;
+ H_sp->val[*H_sp_top] = 0.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*H_sp_top);
- }
-
- /* explicit zeros on diagonals */
- if ( *Htop < H->m )
- {
- H->j[*Htop] = system->N + i;
- H->val[*Htop] = 0.0;
- }
- ++(*Htop);
-
- if ( *H_sp_top < H_sp->m )
- {
- H_sp->j[*H_sp_top] = system->N + i;
- H_sp->val[*H_sp_top] = 0.0;
}
- ++(*H_sp_top);
}
}
break;
@@ -494,110 +553,139 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
for ( i = 0; i < system->N; ++i )
{
- H->start[system->N + i] = *Htop;
- H_sp->start[system->N + i] = *H_sp_top;
-
- /* constraint on ref. value for kinetic energy potential */
- if ( *Htop < H->m )
- {
- H->j[*Htop] = i;
- H->val[*Htop] = 1.0;
- }
- ++(*Htop);
-
- if ( *H_sp_top < H_sp->m )
+ if ( flag_oom == FALSE )
{
- H_sp->j[*H_sp_top] = i;
- H_sp->val[*H_sp_top] = 1.0;
- }
- ++(*H_sp_top);
+ H->start[system->N + i] = *Htop;
+ H_sp->start[system->N + i] = *H_sp_top;
- /* kinetic energy terms */
- for ( pj = Start_Index(i, far_nbr_list); pj < End_Index(i, far_nbr_list); ++pj )
- {
- /* exclude self-periodic images of atoms for
- * kinetic energy term because the effective
- * potential is the same on an atom and its periodic image */
- if ( far_nbr_list->far_nbr_list[pj].d <= control->nonb_cut )
+ /* constraint on ref. value for kinetic energy potential */
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = i;
+ H->val[*Htop] = 1.0;
+ ++(*Htop);
+ }
+ else
{
- j = far_nbr_list->far_nbr_list[pj].nbr;
+ flag_oom = TRUE;
+ }
- xcut = 0.5 * ( system->reax_param.sbp[ system->atoms[i].type ].b_s_acks2
- + system->reax_param.sbp[ system->atoms[j].type ].b_s_acks2 );
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = i;
+ H_sp->val[*H_sp_top] = 1.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ }
- if ( far_nbr_list->far_nbr_list[pj].d < xcut )
+ /* kinetic energy terms */
+ for ( pj = Start_Index(i, far_nbr_list); pj < End_Index(i, far_nbr_list); ++pj )
+ {
+ /* exclude self-periodic images of atoms for
+ * kinetic energy term because the effective
+ * potential is the same on an atom and its periodic image */
+ if ( far_nbr_list->far_nbr_list[pj].d <= control->nonb_cut )
{
- d = far_nbr_list->far_nbr_list[pj].d / xcut;
- bond_softness = system->reax_param.gp.l[34] * POW( d, 3.0 )
- * POW( 1.0 - d, 6.0 );
+ j = far_nbr_list->far_nbr_list[pj].nbr;
+
+ xcut = 0.5 * ( system->reax_param.sbp[ system->atoms[i].type ].b_s_acks2
+ + system->reax_param.sbp[ system->atoms[j].type ].b_s_acks2 );
- if ( bond_softness > 0.0 )
+ if ( far_nbr_list->far_nbr_list[pj].d < xcut )
{
- val_flag = FALSE;
+ d = far_nbr_list->far_nbr_list[pj].d / xcut;
+ bond_softness = system->reax_param.gp.l[34] * POW( d, 3.0 )
+ * POW( 1.0 - d, 6.0 );
- for ( target = H->start[system->N + i]; target < *Htop; ++target )
+ if ( bond_softness > 0.0 )
{
- if ( H->j[target] == system->N + j )
+ val_flag = FALSE;
+
+ for ( target = H->start[system->N + i]; target < *Htop; ++target )
{
- H->val[target] += bond_softness;
- val_flag = TRUE;
- break;
+ if ( H->j[target] == system->N + j )
+ {
+ H->val[target] += bond_softness;
+ val_flag = TRUE;
+ break;
+ }
}
- }
- if ( val_flag == FALSE )
- {
- if ( *Htop < H->m )
+ if ( val_flag == FALSE )
{
- H->j[*Htop] = system->N + j;
- H->val[*Htop] = bond_softness;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N + j;
+ H->val[*Htop] = bond_softness;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
+ }
}
- ++(*Htop);
- }
- val_flag = FALSE;
+ val_flag = FALSE;
- for ( target = H_sp->start[system->N + i]; target < *H_sp_top; ++target )
- {
- if ( H_sp->j[target] == system->N + j )
+ for ( target = H_sp->start[system->N + i]; target < *H_sp_top; ++target )
{
- H_sp->val[target] += bond_softness;
- val_flag = TRUE;
- break;
+ if ( H_sp->j[target] == system->N + j )
+ {
+ H_sp->val[target] += bond_softness;
+ val_flag = TRUE;
+ break;
+ }
}
- }
- if ( val_flag == FALSE )
- {
- if ( *H_sp_top < H_sp->m )
+ if ( val_flag == FALSE )
{
- H_sp->j[*H_sp_top] = system->N + j;
- H_sp->val[*H_sp_top] = bond_softness;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N + j;
+ H_sp->val[*H_sp_top] = bond_softness;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
+ }
}
- ++(*H_sp_top);
- }
- X_diag[i] -= bond_softness;
- X_diag[j] -= bond_softness;
+ X_diag[i] -= bond_softness;
+ X_diag[j] -= bond_softness;
+ }
}
}
}
- }
- /* placeholders for diagonal entries, to be replaced below */
- if ( *Htop < H->m )
- {
- H->j[*Htop] = system->N + i;
- H->val[*Htop] = 0.0;
- }
- ++(*Htop);
+ /* placeholders for diagonal entries, to be replaced below */
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N + i;
+ H->val[*Htop] = 0.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ }
- if ( *H_sp_top < H_sp->m )
- {
- H_sp->j[*H_sp_top] = system->N + i;
- H_sp->val[*H_sp_top] = 0.0;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N + i;
+ H_sp->val[*H_sp_top] = 0.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ }
}
- ++(*H_sp_top);
}
/* second to last row */
@@ -633,15 +721,25 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
{
H->j[*Htop] = system->N + i;
H->val[*Htop] = 1.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
H_sp->j[*H_sp_top] = system->N + i;
H_sp->val[*H_sp_top] = 1.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
}
- ++(*H_sp_top);
}
/* explicitly store zero on diagonal */
@@ -649,15 +747,23 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
{
H->j[*Htop] = system->N_cm - 2;
H->val[*Htop] = 0.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
H_sp->j[*H_sp_top] = system->N_cm - 2;
H_sp->val[*H_sp_top] = 0.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*H_sp_top);
/* last row */
H->start[system->N_cm - 1] = *Htop;
@@ -669,15 +775,25 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
{
H->j[*Htop] = i;
H->val[*Htop] = 1.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
H_sp->j[*H_sp_top] = i;
H_sp->val[*H_sp_top] = 1.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
}
- ++(*H_sp_top);
}
/* explicitly store zero on diagonal */
@@ -685,15 +801,23 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
{
H->j[*Htop] = system->N_cm - 1;
H->val[*Htop] = 0.0;
+ ++(*Htop);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*Htop);
if ( *H_sp_top < H_sp->m )
{
H_sp->j[*H_sp_top] = system->N_cm - 1;
H_sp->val[*H_sp_top] = 0.0;
+ ++(*H_sp_top);
+ }
+ else
+ {
+ flag_oom = TRUE;
}
- ++(*H_sp_top);
sfree( X_diag, __FILE__, __LINE__ );
break;
@@ -701,6 +825,8 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
default:
break;
}
+
+ return flag_oom;
}
@@ -746,7 +872,7 @@ static void Init_CM_Half( reax_system const * const system,
int i, j, pj, target;
int start_i, end_i;
int Htop, H_sp_top;
- int flag, flag_sp, val_flag;
+ int flag, flag_sp, val_flag, flag_oom;
real val;
sparse_matrix *H, *H_sp;
reax_list *far_nbrs;
@@ -756,67 +882,45 @@ static void Init_CM_Half( reax_system const * const system,
H_sp = &workspace->H_sp;
Htop = 0;
H_sp_top = 0;
+ flag_oom = FALSE;
for ( i = 0; i < far_nbrs->n; ++i )
{
- 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;
-
- for ( pj = start_i; pj < end_i; ++pj )
+ if ( flag_oom == FALSE )
{
- j = far_nbrs->far_nbr_list[pj].nbr;
- flag = FALSE;
- flag_sp = FALSE;
-
- if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
- {
- flag = TRUE;
-
- if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_sp_cut )
- {
- flag_sp = TRUE;
- }
- }
+ 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;
- if ( flag == TRUE )
+ for ( pj = start_i; pj < end_i; ++pj )
{
- val = Init_Charge_Matrix_Entry( system, control,
- workspace, i, j, far_nbrs->far_nbr_list[pj].d,
- OFF_DIAGONAL );
- val_flag = FALSE;
+ j = far_nbrs->far_nbr_list[pj].nbr;
+ flag = FALSE;
+ flag_sp = FALSE;
- for ( target = H->start[i]; target < Htop; ++target )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
{
- if ( H->j[target] == j )
- {
- H->val[target] += val;
- val_flag = TRUE;
- break;
- }
- }
+ flag = TRUE;
- if ( val_flag == FALSE )
- {
- if ( Htop < H->m )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_sp_cut )
{
- H->j[Htop] = j;
- H->val[Htop] = val;
+ flag_sp = TRUE;
}
- ++Htop;
}
- /* H_sp matrix entry */
- if ( flag_sp == TRUE )
+ if ( flag == TRUE )
{
+ val = Init_Charge_Matrix_Entry( system, control,
+ workspace, i, j, far_nbrs->far_nbr_list[pj].d,
+ OFF_DIAGONAL );
val_flag = FALSE;
- for ( target = H_sp->start[i]; target < H_sp_top; ++target )
+ for ( target = H->start[i]; target < Htop; ++target )
{
- if ( H_sp->j[target] == j )
+ if ( H->j[target] == j )
{
- H_sp->val[target] += val;
+ H->val[target] += val;
val_flag = TRUE;
break;
}
@@ -824,36 +928,88 @@ static void Init_CM_Half( reax_system const * const system,
if ( val_flag == FALSE )
{
- H_sp->j[H_sp_top] = j;
- H_sp->val[H_sp_top] = val;
- ++H_sp_top;
+ if ( Htop < H->m )
+ {
+ H->j[Htop] = j;
+ H->val[Htop] = val;
+ ++Htop;
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
+ }
+ }
+
+ /* H_sp matrix entry */
+ if ( flag_sp == TRUE )
+ {
+ val_flag = FALSE;
+
+ for ( target = H_sp->start[i]; target < H_sp_top; ++target )
+ {
+ if ( H_sp->j[target] == j )
+ {
+ H_sp->val[target] += val;
+ val_flag = TRUE;
+ break;
+ }
+ }
+
+ if ( val_flag == FALSE )
+ {
+ if ( H_sp_top < H_sp->m )
+ {
+ H_sp->j[H_sp_top] = j;
+ H_sp->val[H_sp_top] = val;
+ ++H_sp_top;
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
+ }
+ }
}
}
}
- }
- /* diagonal entry */
- if ( Htop < H->m )
- {
- H->j[Htop] = i;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
- workspace, i, i, 0.0, DIAGONAL );
- }
- ++Htop;
+ /* diagonal entry */
+ if ( Htop < H->m )
+ {
+ H->j[Htop] = i;
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
+ workspace, i, i, 0.0, DIAGONAL );
+ ++Htop;
+ }
+ else
+ {
+ flag_oom = TRUE;
+ }
- H_sp->j[H_sp_top] = i;
- H_sp->val[H_sp_top] = H->val[Htop - 1];
- ++H_sp_top;
+ if ( H_sp_top < H_sp->m )
+ {
+ H_sp->j[H_sp_top] = i;
+ H_sp->val[H_sp_top] = H->val[Htop - 1];
+ ++H_sp_top;
+ }
+ else
+ {
+ flag_oom = TRUE;
+ }
+ }
}
- Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
- H, H_sp, &Htop, &H_sp_top );
+ if ( flag_oom == FALSE )
+ {
+ flag_oom = Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
+ H, H_sp, &Htop, &H_sp_top );
+ }
H->start[system->N_cm] = Htop;
H_sp->start[system->N_cm] = H_sp_top;
-
- if ( Htop > H->m )
+ if ( flag_oom == TRUE )
{
workspace->realloc.cm = TRUE;
}
@@ -870,7 +1026,7 @@ static void Init_CM_Tab_Half( reax_system const * const system,
int i, j, pj, target;
int start_i, end_i;
int Htop, H_sp_top;
- int flag, flag_sp, val_flag;
+ int flag, flag_sp, val_flag, flag_oom;
real val;
sparse_matrix *H, *H_sp;
reax_list *far_nbrs;
@@ -880,67 +1036,45 @@ static void Init_CM_Tab_Half( reax_system const * const system,
H_sp = &workspace->H_sp;
Htop = 0;
H_sp_top = 0;
+ flag_oom = FALSE;
for ( i = 0; i < far_nbrs->n; ++i )
{
- 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;
-
- for ( pj = start_i; pj < end_i; ++pj )
+ if ( flag_oom == FALSE )
{
- j = far_nbrs->far_nbr_list[pj].nbr;
- flag = FALSE;
- flag_sp = FALSE;
-
- if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
- {
- flag = TRUE;
-
- if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_sp_cut )
- {
- flag_sp = TRUE;
- }
- }
+ 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;
- if ( flag == TRUE )
+ for ( pj = start_i; pj < end_i; ++pj )
{
- val = Init_Charge_Matrix_Entry_Tab( system, control,
- workspace->LR, i, j, far_nbrs->far_nbr_list[pj].d,
- OFF_DIAGONAL );
- val_flag = FALSE;
+ j = far_nbrs->far_nbr_list[pj].nbr;
+ flag = FALSE;
+ flag_sp = FALSE;
- for ( target = H->start[i]; target < Htop; ++target )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
{
- if ( H->j[target] == j )
- {
- H->val[target] += val;
- val_flag = TRUE;
- break;
- }
- }
+ flag = TRUE;
- if ( val_flag == FALSE )
- {
- if ( Htop < H->m )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_sp_cut )
{
- H->j[Htop] = j;
- H->val[Htop] = val;
+ flag_sp = TRUE;
}
- ++Htop;
}
- /* H_sp matrix entry */
- if ( flag_sp == TRUE )
+ if ( flag == TRUE )
{
+ val = Init_Charge_Matrix_Entry_Tab( system, control,
+ workspace->LR, i, j, far_nbrs->far_nbr_list[pj].d,
+ OFF_DIAGONAL );
val_flag = FALSE;
- for ( target = H_sp->start[i]; target < H_sp_top; ++target )
+ for ( target = H->start[i]; target < Htop; ++target )
{
- if ( H_sp->j[target] == j )
+ if ( H->j[target] == j )
{
- H_sp->val[target] += val;
+ H->val[target] += val;
val_flag = TRUE;
break;
}
@@ -948,36 +1082,88 @@ static void Init_CM_Tab_Half( reax_system const * const system,
if ( val_flag == FALSE )
{
- H_sp->j[H_sp_top] = j;
- H_sp->val[H_sp_top] = val;
- ++H_sp_top;
+ if ( Htop < H->m )
+ {
+ H->j[Htop] = j;
+ H->val[Htop] = val;
+ ++Htop;
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
+ }
+ }
+
+ /* H_sp matrix entry */
+ if ( flag_sp == TRUE )
+ {
+ val_flag = FALSE;
+
+ for ( target = H_sp->start[i]; target < H_sp_top; ++target )
+ {
+ if ( H_sp->j[target] == j )
+ {
+ H_sp->val[target] += val;
+ val_flag = TRUE;
+ break;
+ }
+ }
+
+ if ( val_flag == FALSE )
+ {
+ if ( H_sp_top < H_sp->m )
+ {
+ H_sp->j[H_sp_top] = j;
+ H_sp->val[H_sp_top] = val;
+ ++H_sp_top;
+ }
+ else
+ {
+ flag_oom = TRUE;
+ break;
+ }
+ }
}
}
}
- }
- /* diagonal entry */
- if ( Htop < H->m )
- {
- H->j[Htop] = i;
- H->val[Htop] = Init_Charge_Matrix_Entry_Tab( system, control,
- workspace->LR, i, i, 0.0, DIAGONAL );
- }
- ++Htop;
+ /* diagonal entry */
+ if ( Htop < H->m )
+ {
+ H->j[Htop] = i;
+ H->val[Htop] = Init_Charge_Matrix_Entry_Tab( system, control,
+ workspace->LR, i, i, 0.0, DIAGONAL );
+ ++Htop;
+ }
+ else
+ {
+ flag_oom = TRUE;
+ }
- H_sp->j[H_sp_top] = i;
- H_sp->val[H_sp_top] = H->val[Htop - 1];
- ++H_sp_top;
+ if ( H_sp_top < H_sp->m )
+ {
+ H_sp->j[H_sp_top] = i;
+ H_sp->val[H_sp_top] = H->val[Htop - 1];
+ ++H_sp_top;
+ }
+ else
+ {
+ flag_oom = TRUE;
+ }
+ }
}
- Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
- H, H_sp, &Htop, &H_sp_top );
-
- H->start[system->N_cm] = Htop;
- H_sp->start[system->N_cm] = H_sp_top;
+ if ( flag_oom == FALSE )
+ {
+ flag_oom = Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
+ H, H_sp, &Htop, &H_sp_top );
+ H->start[system->N_cm] = Htop;
+ H_sp->start[system->N_cm] = H_sp_top;
+ }
- if ( Htop > H->m )
+ if ( flag_oom == TRUE )
{
workspace->realloc.cm = TRUE;
}
@@ -995,7 +1181,7 @@ static void Init_Bond_Full( reax_system const * const system,
int type_i, type_j;
int btop_i, btop_j;
int ihb, jhb, ihb_top, jhb_top;
- int num_bonds, num_hbonds;
+ int num_bonds, num_hbonds, flag_oom_bonds, flag_oom_hbonds;
real r_ij, r2;
real C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
@@ -1014,6 +1200,8 @@ static void Init_Bond_Full( reax_system const * const system,
num_hbonds = 0;
btop_i = 0;
btop_j = 0;
+ flag_oom_bonds = FALSE;
+ flag_oom_hbonds = FALSE;
for ( i = 0; i < far_nbrs->n; ++i )
{
@@ -1052,7 +1240,7 @@ static void Init_Bond_Full( reax_system const * const system,
|| system->atoms[j].qmmm_mask == TRUE )
{
#endif
- if ( nbr_pj->d <= control->nonb_cut )
+ if ( nbr_pj->d <= control->nonb_cut )
{
type_j = system->atoms[j].type;
sbp_j = &system->reax_param.sbp[type_j];
@@ -1083,8 +1271,12 @@ static void Init_Bond_Full( reax_system const * const system,
hbonds->hbond_list[ihb_top].scl = 1;
hbonds->hbond_list[ihb_top].ptr = nbr_pj;
++ihb_top;
+ ++num_hbonds;
+ }
+ else
+ {
+ flag_oom_hbonds = TRUE;
}
- ++num_hbonds;
}
else if ( ihb == H_BONDING_ATOM && jhb == H_ATOM )
{
@@ -1095,8 +1287,12 @@ static void Init_Bond_Full( reax_system const * const system,
hbonds->hbond_list[jhb_top].scl = -1;
hbonds->hbond_list[jhb_top].ptr = nbr_pj;
Set_End_Index( j, jhb_top + 1, hbonds );
+ ++num_hbonds;
+ }
+ else
+ {
+ flag_oom_hbonds = TRUE;
}
- ++num_hbonds;
}
}
@@ -1215,8 +1411,12 @@ static void Init_Bond_Full( reax_system const * const system,
bo_ji->Cdbopi2 = 0.0;
Set_End_Index( j, btop_j + 1, bonds );
+ num_bonds += 2;
+ }
+ else
+ {
+ flag_oom_bonds = TRUE;
}
- num_bonds += 2;
}
}
}
@@ -1236,27 +1436,15 @@ static void Init_Bond_Full( reax_system const * const system,
}
}
- for ( i = 0; i < bonds->n; ++i )
+ if ( flag_oom_bonds == TRUE )
{
- if ( Num_Entries( i, bonds ) > system->bonds[i] )
- {
- workspace->realloc.bonds = TRUE;
- break;
- }
-
+ workspace->realloc.bonds = TRUE;
}
- if ( control->hbond_cut > 0.0 && workspace->num_H > 0 )
+ if ( control->hbond_cut > 0.0 && workspace->num_H > 0
+ && flag_oom_hbonds == TRUE )
{
- for ( i = 0; i < hbonds->n; ++i )
- {
- if ( Num_Entries( i, hbonds ) > system->hbonds[i] )
- {
- workspace->realloc.hbonds = TRUE;
- break;
- }
-
- }
+ workspace->realloc.hbonds = TRUE;
}
}
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index 49accae1128f0f1c2a4875fb9baf1ade9bc364fe..be82ce7b0f39b8b96ec21e0a0b8fa1feadda9c8b 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -709,8 +709,8 @@ static void reax_atom_Copy( reax_atom * const dest, reax_atom * const src )
static void Copy_Storage( reax_system const * const system,static_storage * const workspace,
control_params const * const control, int top, int old_id, int old_type,
- int * const num_H, real ** const v, real ** const s, real ** const t,
- int * const orig_id, rvec * const f_old )
+ real ** const v, real ** const s, real ** const t, int * const orig_id,
+ rvec * const f_old )
{
int i;
@@ -772,7 +772,7 @@ static void Assign_New_Storage( static_storage *workspace,
void Reorder_Atoms( reax_system * const system, static_storage * const workspace,
control_params const * const control )
{
- int i, j, k, l, top, old_id, num_H;
+ int i, j, k, l, top, old_id;
reax_atom *old_atom, *new_atoms;
grid *g;
int *orig_id;
@@ -780,7 +780,6 @@ void Reorder_Atoms( reax_system * const system, static_storage * const workspace
real **s, **t;
rvec *f_old;
- num_H = 0;
top = 0;
g = &system->g;
@@ -818,7 +817,7 @@ void Reorder_Atoms( reax_system * const system, static_storage * const workspace
reax_atom_Copy( &new_atoms[top], old_atom );
Copy_Storage( system, workspace, control, top, old_id, old_atom->type,
- &num_H, v, s, t, orig_id, f_old );
+ v, s, t, orig_id, f_old );
++top;
}
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 29db67758b46efbf658bdab1d297c0adf0282584..2a0e9dbbc0262ebbd0ee3ec3cf14dd908094e6b3 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -854,7 +854,13 @@ static void Init_Lists( reax_system * const system,
MAX( workspace->realloc.total_far_nbrs, lists[FAR_NBRS]->total_intrs ),
TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
- Generate_Neighbor_Lists( system, control, data, workspace, lists );
+ ret = Generate_Neighbor_Lists( system, control, data, workspace, lists );
+
+ if ( ret != SUCCESS )
+ {
+ fprintf( stderr, "[ERROR] Unrecoverable memory allocation issue (Generate_Neighbor_Lists). Terminating...\n" );
+ exit( INVALID_INPUT );
+ }
}
if ( realloc == TRUE )
@@ -917,7 +923,7 @@ static void Init_Lists( reax_system * const system,
{
if ( system->reax_param.sbp[ system->atoms[i].type ].p_hbond == H_ATOM )
{
- workspace->num_H++;
+ ++(workspace->num_H);
}
}
}
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index c75db365b6b1e4c45c9c0fac1afd2caffaaa40ce..3fe9919161455678ff39429276d882ba85385cef 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -220,7 +220,7 @@ int Generate_Neighbor_Lists( reax_system * const system,
int i, j, k, l, m, itr;
int x, y, z;
int atom1, atom2, max;
- int num_far, count, ret;
+ int num_far, count, flag_oom, ret;
int *nbr_atoms;
ivec *nbrs;
rvec *nbrs_cp;
@@ -232,6 +232,7 @@ int Generate_Neighbor_Lists( reax_system * const system,
t_start = Get_Time( );
num_far = 0;
far_nbrs = lists[FAR_NBRS];
+ flag_oom = FALSE;
ret = SUCCESS;
Choose_Neighbor_Finder( system, control, &Find_Far_Neighbors );
@@ -292,6 +293,11 @@ int Generate_Neighbor_Lists( reax_system * const system,
far_nbrs->total_intrs - num_far );
num_far += count;
+
+ if ( num_far >= far_nbrs->total_intrs )
+ {
+ goto OUT_OF_MEMORY;
+ }
}
}
}
@@ -311,9 +317,10 @@ int Generate_Neighbor_Lists( reax_system * const system,
ivec_MakeZero( system->atoms[i].rel_map );
}
- if ( num_far > far_nbrs->total_intrs )
+ if ( flag_oom == TRUE )
{
- ret = FAILURE;
+ OUT_OF_MEMORY:
+ ret = FAILURE;
}
#if defined(DEBUG_FOCUS)