From 90f7d97bdc0047bb3179cc0023ad354da80f4eff Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Thu, 19 Oct 2017 22:19:49 -0400
Subject: [PATCH] sPuReMD: add safe memory allocation functions. Fix issue with
certain parameters for shielded van der Waals interactions not being
correctly initialized. Other general refactoring.
---
sPuReMD/src/allocate.c | 34 ++-
sPuReMD/src/analyze.c | 8 +-
sPuReMD/src/box.c | 32 +--
sPuReMD/src/charges.c | 44 ++--
sPuReMD/src/control.c | 6 +-
sPuReMD/src/ffield.c | 37 ++-
sPuReMD/src/forces.c | 125 ++--------
sPuReMD/src/geo_tools.c | 70 +++---
sPuReMD/src/grid.c | 119 +++++-----
sPuReMD/src/init_md.c | 347 +++++++++++++++++-----------
sPuReMD/src/lin_alg.c | 197 ++++++++++++----
sPuReMD/src/list.c | 25 +-
sPuReMD/src/lookup.c | 198 ++++++++++------
sPuReMD/src/lookup.h | 2 +
sPuReMD/src/mytypes.h | 1 -
sPuReMD/src/neighbors.c | 16 +-
sPuReMD/src/reset_utils.c | 6 +
sPuReMD/src/system_props.h | 4 +-
sPuReMD/src/testmd.c | 17 +-
sPuReMD/src/tool_box.c | 149 +++++++++---
sPuReMD/src/tool_box.h | 47 +++-
sPuReMD/src/two_body_interactions.c | 7 +-
22 files changed, 910 insertions(+), 581 deletions(-)
diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 6da598bd..5a1cbfae 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "allocate.h"
+
#include "list.h"
#include "tool_box.h"
@@ -101,10 +102,10 @@ int Allocate_Matrix( sparse_matrix **pH, int n, int m )
/* deallocate memory for matrix in CSR format */
void Deallocate_Matrix( sparse_matrix *H )
{
- free(H->start);
- free(H->j);
- free(H->val);
- free(H);
+ sfree( H->start, "Deallocate_Matrix::H->start" );
+ sfree( H->j, "Deallocate_Matrix::H->j" );
+ sfree( H->val, "Deallocate_Matrix::H->val" );
+ sfree( H, "Deallocate_Matrix::H" );
}
@@ -130,7 +131,7 @@ int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name )
int Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
- list *hbonds )
+ list *hbonds )
{
int i, num_hbonds;
@@ -179,7 +180,8 @@ int Reallocate_HBonds_List( int n, int num_h, int *h_index, list *hbonds )
#if defined(DEBUG_FOCUS)
fprintf( stderr, "reallocating hbonds\n" );
#endif
- hb_top = calloc( n, sizeof(int) );
+
+ hb_top = (int *) calloc( n, sizeof(int) );
for ( i = 0; i < n; ++i )
{
if ( h_index[i] >= 0 )
@@ -192,7 +194,7 @@ int Reallocate_HBonds_List( int n, int num_h, int *h_index, list *hbonds )
Allocate_HBond_List( n, num_h, h_index, hb_top, hbonds );
- free( hb_top );
+ sfree( hb_top, "Reallocate_HBonds_List::hb_top" );
return SUCCESS;
}
@@ -241,8 +243,10 @@ int Reallocate_Bonds_List( int n, list *bonds, int *num_bonds, int *est_3body )
#if defined(DEBUG_FOCUS)
fprintf( stderr, "reallocating bonds\n" );
#endif
- bond_top = calloc( n, sizeof(int) );
+
+ bond_top = (int *) calloc( n, sizeof(int) );
*est_3body = 0;
+
for ( i = 0; i < n; ++i )
{
*est_3body += SQR( Num_Entries( i, bonds ) );
@@ -254,14 +258,14 @@ int Reallocate_Bonds_List( int n, list *bonds, int *num_bonds, int *est_3body )
Allocate_Bond_List( n, bond_top, bonds );
*num_bonds = bond_top[n - 1];
- free( bond_top );
+ sfree( bond_top, "Reallocate_Bonds_List::bond_top" );
return SUCCESS;
}
void Reallocate( reax_system *system, static_storage *workspace, list **lists,
- int nbr_flag )
+ int nbr_flag )
{
int i, j, k;
int num_bonds, est_3body;
@@ -334,15 +338,21 @@ void Reallocate( reax_system *system, static_storage *workspace, list **lists,
#if defined(DEBUG_FOCUS)
fprintf(stderr, "reallocating gcell: g->max_atoms: %d\n", g->max_atoms);
#endif
+
for ( i = 0; i < g->ncell[0]; i++ )
+ {
for ( j = 0; j < g->ncell[1]; j++ )
+ {
for ( k = 0; k < g->ncell[2]; k++ )
{
// reallocate g->atoms
- free( g->atoms[i][j][k] );
+ sfree( g->atoms[i][j][k], "Reallocate::g->atoms[i][j][k]" );
g->atoms[i][j][k] = (int*)
- calloc(workspace->realloc.gcell_atoms, sizeof(int));
+ calloc(workspace->realloc.gcell_atoms, sizeof(int));
}
+ }
+ }
+
realloc->gcell_atoms = -1;
}
}
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index e78f533d..4bb9cf4b 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -217,8 +217,8 @@ void Print_Molecule( reax_system *system, molecule *m, int mode, char *s )
void Analyze_Molecules( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- list **lists, FILE *fout )
+ simulation_data *data, static_storage *workspace,
+ list **lists, FILE *fout )
{
int atom, i, j, k, l, flag;
int *mark = workspace->mark;
@@ -304,8 +304,8 @@ void Analyze_Molecules( reax_system *system, control_params *control,
Copy_Bond_List( system, control, lists );
- //free( mark );
- //free( old_mark );
+ //sfree( mark, "Analyze_Molecules::mark" );
+ //sfree( old_mark, "Analyze_Molecules::old_mark" );
fprintf( fout, "\n" );
fflush( fout );
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index 3cddaf2d..f3646ae3 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -162,10 +162,10 @@ void Setup_Box( real a, real b, real c, real alpha, real beta, real gamma,
exit( INVALID_INPUT );
}
- c_alpha = COS(DEG2RAD(alpha));
- c_beta = COS(DEG2RAD(beta));
- c_gamma = COS(DEG2RAD(gamma));
- s_gamma = SIN(DEG2RAD(gamma));
+ c_alpha = COS( DEG2RAD(alpha) );
+ c_beta = COS( DEG2RAD(beta) );
+ c_gamma = COS( DEG2RAD(gamma) );
+ s_gamma = SIN( DEG2RAD(gamma) );
zi = (c_alpha - c_beta * c_gamma) / s_gamma;
box->box[0][0] = a;
@@ -177,6 +177,7 @@ void Setup_Box( real a, real b, real c, real alpha, real beta, real gamma,
box->box[2][0] = c * c_beta;
box->box[2][1] = c * zi;
box->box[2][2] = c * SQRT(1.0 - SQR(c_beta) - SQR(zi));
+
#if defined(DEBUG)
fprintf( stderr, "box is %8.2f x %8.2f x %8.2f\n",
box->box[0][0], box->box[1][1], box->box[2][2] );
@@ -215,7 +216,7 @@ void Update_Box_Isotropic( simulation_box *box, real mu )
box->box[2][2] *= mu;
box->volume = box->box[0][0] * box->box[1][1] * box->box[2][2];
- Make_Consistent(box/*, periodic*/);
+ Make_Consistent( box );
}
@@ -231,7 +232,7 @@ void Update_Box_SemiIsotropic( simulation_box *box, rvec mu )
box->box[2][2] *= mu[2];
box->volume = box->box[0][0] * box->box[1][1] * box->box[2][2];
- Make_Consistent(box);
+ Make_Consistent( box );
}
@@ -240,18 +241,20 @@ void Inc_on_T3( rvec x, rvec dx, simulation_box *box )
int i;
real tmp;
- for (i = 0; i < 3; i++)
+ for ( i = 0; i < 3; i++ )
{
tmp = x[i] + dx[i];
+
if ( tmp <= -box->box_norms[i] || tmp >= box->box_norms[i] )
{
tmp = FMOD( tmp, box->box_norms[i] );
}
- if ( tmp < 0 )
+ if ( tmp < 0.0 )
{
tmp += box->box_norms[i];
}
+
x[i] = tmp;
}
}
@@ -349,15 +352,16 @@ real Metric_Product( rvec x1, rvec x2, simulation_box* box )
int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
real cutoff, far_neighbor_data *data )
{
- real norm_sqr, d, tmp;
+ real norm_sqr, d, tmp, ret;
int i;
- norm_sqr = 0;
+ norm_sqr = 0.0;
+ ret = FALSE;
for ( i = 0; i < 3; i++ )
{
d = x2[i] - x1[i];
- tmp = SQR(d);
+ tmp = SQR( d );
if ( tmp >= SQR( box->box_norms[i] / 2.0 ) )
{
@@ -373,7 +377,7 @@ int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
}
data->dvec[i] = d;
- norm_sqr += SQR(d);
+ norm_sqr += SQR( d );
}
else
{
@@ -386,10 +390,10 @@ int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
if ( norm_sqr <= SQR( cutoff ) )
{
data->d = SQRT( norm_sqr );
- return 1;
+ ret = TRUE;
}
- return 0;
+ return ret;
}
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 07641a69..63ed940d 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -28,7 +28,7 @@
#include "tool_box.h"
#include "vector.h"
#if defined(HAVE_SUPERLU_MT)
-#include "slu_mt_ddefs.h"
+ #include "slu_mt_ddefs.h"
#endif
@@ -142,7 +142,7 @@ static void Sort_Matrix_Rows( sparse_matrix * const A )
}
}
- free( temp );
+ sfree( temp, "Sort_Matrix_Rows::temp" );
}
}
@@ -287,7 +287,7 @@ static int Estimate_LU_Fill( const sparse_matrix * const A, const real * const d
#if defined(HAVE_SUPERLU_MT)
static real SuperLU_Factorize( const sparse_matrix * const A,
- sparse_matrix * const L, sparse_matrix * const U )
+ sparse_matrix * const L, sparse_matrix * const U )
{
unsigned int i, pj, count, *Ltop, *Utop, r;
sparse_matrix *A_t;
@@ -554,9 +554,9 @@ static real SuperLU_Factorize( const sparse_matrix * const A,
------------------------------------------------------------*/
pxgstrf_finalize( &superlumt_options, &AC_S );
Deallocate_Matrix( A_t );
- free( xa );
- free( asub );
- free( a );
+ sfree( xa, "SuperLU_Factorize::xa" );
+ sfree( asub, "SuperLU_Factorize::asub" );
+ sfree( a, "SuperLU_Factorize::a" );
SUPERLU_FREE( perm_r );
SUPERLU_FREE( perm_c );
SUPERLU_FREE( ((NCformat *)A_S.Store)->rowind );
@@ -574,8 +574,8 @@ static real SuperLU_Factorize( const sparse_matrix * const A,
}
StatFree(&Gstat);
- free( Utop );
- free( Ltop );
+ sfree( Utop, "SuperLU_Factorize::Utop" );
+ sfree( Ltop, "SuperLU_Factorize::Ltop" );
//TODO: return iters
return 0.;
@@ -750,9 +750,9 @@ static real ICHOLT( const sparse_matrix * const A, const real * const droptol,
// fprintf( stderr, "nnz(U): %d, max: %d\n", Utop[U->n], U->n * 50 );
- free( tmp_val );
- free( tmp_j );
- free( Utop );
+ sfree( tmp_val, "ICHOLT::tmp_val" );
+ sfree( tmp_j, "ICHOLT::tmp_j" );
+ sfree( Utop, "ICHOLT::Utop" );
return Get_Timing_Info( start );
}
@@ -766,7 +766,7 @@ static real ICHOLT( const sparse_matrix * const A, const real * const droptol,
* SIAM J. Sci. Comp. */
#if defined(TESTING)
static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
- sparse_matrix * const U_t, sparse_matrix * const U )
+ sparse_matrix * const U_t, sparse_matrix * const U )
{
unsigned int i, j, k, pj, x = 0, y = 0, ei_x, ei_y;
real *D, *D_inv, sum, start;
@@ -925,9 +925,9 @@ static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
#endif
Deallocate_Matrix( DAD );
- free(D_inv);
- free(D);
- free(Utop);
+ sfree( D_inv, "ICHOL_PAR::D_inv" );
+ sfree( D, "ICHOL_PAR::D" );
+ sfree( Utop, "ICHOL_PAR::Utop" );
return Get_Timing_Info( start );
}
@@ -1141,8 +1141,8 @@ static real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
#endif
Deallocate_Matrix( DAD );
- free( D_inv );
- free( D );
+ sfree( D_inv, "ILU_PAR::D_inv" );
+ sfree( D, "ILU_PAR::D_inv" );
return Get_Timing_Info( start );
}
@@ -1407,8 +1407,8 @@ static real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
Deallocate_Matrix( U_temp );
Deallocate_Matrix( L_temp );
Deallocate_Matrix( DAD );
- free( D_inv );
- free( D );
+ sfree( D_inv, "ILUT_PAR::D_inv" );
+ sfree( D, "ILUT_PAR::D_inv" );
return Get_Timing_Info( start );
}
@@ -2540,6 +2540,12 @@ static void Calculate_Charges_EE( const reax_system * const system,
for ( i = 0; i < system->N; ++i )
{
system->atoms[i].q = workspace->s[0][i];
+
+#if defined(DEBUG_FOCUS)
+ printf( "atom %4d: %f\n", i, system->atoms[i].q );
+ printf( " x[0]: %10.5f, x[1]: %10.5f, x[2]: %10.5f\n",
+ system->atoms[i].x[0], system->atoms[i].x[1], system->atoms[i].x[2] );
+#endif
}
}
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index ad2d3d21..ee07a2ce 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -580,10 +580,10 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
/* free memory allocations at the top */
for ( i = 0; i < MAX_TOKENS; i++ )
{
- free( tmp[i] );
+ sfree( tmp[i], "Read_Control_File::tmp[i]" );
}
- free( tmp );
- free( s );
+ sfree( tmp, "Read_Control_File::tmp" );
+ sfree( s, "Read_Control_File::s" );
#if defined(DEBUG_FOCUS)
fprintf( stderr,
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index a39aacf1..226faa40 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -457,6 +457,29 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
POW(reax->sbp[j].gamma *
reax->sbp[i].gamma, -1.5);
+ reax->tbp[i][j].acore =
+ SQRT( reax->sbp[i].acore2 *
+ reax->sbp[j].acore2 );
+
+ reax->tbp[j][i].acore =
+ SQRT( reax->sbp[j].acore2 *
+ reax->sbp[i].acore2 );
+
+ reax->tbp[i][j].ecore =
+ SQRT( reax->sbp[i].ecore2 *
+ reax->sbp[j].ecore2 );
+
+ reax->tbp[j][i].ecore =
+ SQRT( reax->sbp[j].ecore2 *
+ reax->sbp[i].ecore2 );
+
+ reax->tbp[i][j].rcore =
+ SQRT( reax->sbp[i].rcore2 *
+ reax->sbp[j].rcore2 );
+
+ reax->tbp[j][i].rcore =
+ SQRT( reax->sbp[j].rcore2 *
+ reax->sbp[i].rcore2 );
}
/* next line is number of 2-body offdiagonal combinations and some comments */
@@ -719,10 +742,10 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
/* deallocate helper storage */
for ( i = 0; i < MAX_TOKENS; i++ )
{
- free( tmp[i] );
+ sfree( tmp[i], "Read_Force_Field::tmp[i]" );
}
- free( tmp );
- free( s );
+ sfree( tmp, "Read_Force_Field::tmp" );
+ sfree( s, "Read_Force_Field::s" );
/* deallocate tor_flag */
for ( i = 0; i < reax->num_atom_types; i++ )
@@ -731,16 +754,16 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
{
for ( k = 0; k < reax->num_atom_types; k++ )
{
- free( tor_flag[i][j][k] );
+ sfree( tor_flag[i][j][k], "Read_Force_Field::tor_flag[i][j][k]" );
}
- free( tor_flag[i][j] );
+ sfree( tor_flag[i][j], "Read_Force_Field::tor_flag[i][j]" );
}
- free( tor_flag[i] );
+ sfree( tor_flag[i], "Read_Force_Field::tor_flag[i]" );
}
- free( tor_flag );
+ sfree( tor_flag, "Read_Force_Field::tor_flag" );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "force field read\n" );
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index b5c1a830..09b92b34 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -23,14 +23,14 @@
#include "box.h"
#include "bond_orders.h"
-#include "single_body_interactions.h"
-#include "two_body_interactions.h"
-#include "three_body_interactions.h"
+#include "charges.h"
#include "four_body_interactions.h"
#include "list.h"
-#include "print_utils.h"
#include "system_props.h"
-#include "charges.h"
+#include "single_body_interactions.h"
+#include "three_body_interactions.h"
+#include "tool_box.h"
+#include "two_body_interactions.h"
#include "vector.h"
@@ -162,6 +162,7 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
Tabulated_vdW_Coulomb_Energy( system, control, data, workspace,
lists, out_control );
}
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "nonb forces - " );
#endif
@@ -338,6 +339,9 @@ static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
switch ( control->charge_method )
{
case QEQ_CM:
+ //TODO: tabulate other portions of matrices for EE, ACKS2?
+ case EE_CM:
+ case ACKS2_CM:
switch ( pos )
{
case OFF_DIAGONAL:
@@ -347,48 +351,23 @@ static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
/* cubic spline interpolation */
r = (int)(r_ij * t->inv_dx);
- if ( r == 0 ) ++r;
+ if ( r == 0 )
+ {
+ ++r;
+ }
base = (real)(r + 1) * t->dx;
dif = r_ij - base;
val = ((t->ele[r].d * dif + t->ele[r].c) * dif + t->ele[r].b) * dif +
- t->ele[r].a;
+ t->ele[r].a;
val *= EV_to_KCALpMOL / C_ele;
ret = ((i == j) ? 0.5 : 1.0) * val;
break;
- case DIAGONAL:
- ret = system->reaxprm.sbp[system->atoms[i].type].eta;
- break;
- default:
- fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
- }
- break;
- case EE_CM:
- //TODO
- switch ( pos )
- {
- case OFF_DIAGONAL:
- break;
case DIAGONAL:
+ ret = system->reaxprm.sbp[system->atoms[i].type].eta;
break;
- default:
- fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
- }
- break;
- case ACKS2_CM:
- //TODO
- switch ( pos )
- {
- case OFF_DIAGONAL:
- break;
- case DIAGONAL:
- break;
default:
fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
exit( INVALID_INPUT );
@@ -417,6 +396,8 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
switch ( control->charge_method )
{
case QEQ_CM:
+ case EE_CM:
+ case ACKS2_CM:
switch ( pos )
{
case OFF_DIAGONAL:
@@ -447,76 +428,6 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
}
break;
- case EE_CM:
- switch ( pos )
- {
- case OFF_DIAGONAL:
- if ( r_ij < control->r_cut && r_ij > 0.001 )
- {
- Tap = control->Tap7 * r_ij + control->Tap6;
- Tap = Tap * r_ij + control->Tap5;
- Tap = Tap * r_ij + control->Tap4;
- Tap = Tap * r_ij + control->Tap3;
- Tap = Tap * r_ij + control->Tap2;
- Tap = Tap * r_ij + control->Tap1;
- Tap = Tap * r_ij + control->Tap0;
-
- gamij = SQRT( system->reaxprm.sbp[system->atoms[i].type].gamma
- * system->reaxprm.sbp[system->atoms[j].type].gamma );
- /* shielding */
- dr3gamij_1 = POW( r_ij, 3.0 ) + 1.0 / POW( gamij, 3.0 );
- dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
-
- ret = Tap * EV_to_KCALpMOL / dr3gamij_3;
- }
- break;
-
- case DIAGONAL:
- ret = system->reaxprm.sbp[system->atoms[i].type].eta;
- break;
-
- default:
- fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
- }
- break;
-
- case ACKS2_CM:
- switch ( pos )
- {
- case OFF_DIAGONAL:
- if ( r_ij < control->r_cut && r_ij > 0.001 )
- {
- Tap = control->Tap7 * r_ij + control->Tap6;
- Tap = Tap * r_ij + control->Tap5;
- Tap = Tap * r_ij + control->Tap4;
- Tap = Tap * r_ij + control->Tap3;
- Tap = Tap * r_ij + control->Tap2;
- Tap = Tap * r_ij + control->Tap1;
- Tap = Tap * r_ij + control->Tap0;
-
- gamij = SQRT( system->reaxprm.sbp[system->atoms[i].type].gamma
- * system->reaxprm.sbp[system->atoms[j].type].gamma );
- /* shielding */
- dr3gamij_1 = POW( r_ij, 3.0 ) + 1.0 / POW( gamij, 3.0 );
- dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
-
- ret = Tap * EV_to_KCALpMOL / dr3gamij_3;
- }
- break;
-
- case DIAGONAL:
- ret = system->reaxprm.sbp[system->atoms[i].type].eta;
- break;
-
- default:
- fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
- }
- break;
-
default:
fprintf( stderr, "Invalid charge method. Terminating...\n" );
exit( INVALID_INPUT );
@@ -688,7 +599,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
H_sp->val[*H_sp_top] = 0.0;
*H_sp_top = *H_sp_top + 1;
- free( X_diag );
+ sfree( X_diag, "Init_Charge_Matrix_Remaining_Entries::X_diag" );
break;
default:
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 442f1e0c..dbcb9ab8 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -19,13 +19,10 @@
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
-#include <ctype.h>
-
#include "geo_tools.h"
+
#include "allocate.h"
#include "box.h"
-#include "list.h"
-#include "restart.h"
#include "tool_box.h"
#include "vector.h"
@@ -59,7 +56,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
char Read_Geo( char* geo_file, reax_system* system, control_params *control,
- simulation_data *data, static_storage *workspace )
+ simulation_data *data, static_storage *workspace )
{
FILE *geo;
@@ -128,22 +125,28 @@ char Read_Geo( char* geo_file, reax_system* system, control_params *control,
int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
{
+ int ret;
char *cryst;
- char line[MAX_LINE + 1];
- char descriptor[9];
- char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
- char s_group[12], s_zValue[12];
+ char line[MAX_LINE + 1];
+ char descriptor[9];
+ char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
+ char s_group[12], s_zValue[12];
+
+ /* set the pointer to the beginning of the file */
+ fseek( geo, 0, SEEK_SET );
/* initialize variables */
- fseek( geo, 0, SEEK_SET ); // set the pointer to the beginning of the file
+ ret = FAILURE;
switch ( geo_format )
{
case PDB:
cryst = "CRYST1";
break;
+
default:
cryst = "BOX";
+ break;
}
/* locate the cryst line in the geo file, read it and
@@ -153,25 +156,28 @@ int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
if ( strncmp( line, cryst, 6 ) == 0 )
{
if ( geo_format == PDB )
+ {
sscanf( line, PDB_CRYST1_FORMAT,
&descriptor[0],
&s_a[0], &s_b[0], &s_c[0],
&s_alpha[0], &s_beta[0], &s_gamma[0],
&s_group[0], &s_zValue[0] );
- /* compute full volume tensor from the angles */
- Setup_Box( atof(s_a), atof(s_b), atof(s_c),
- atof(s_alpha), atof(s_beta), atof(s_gamma),
- &(system->box) );
- return SUCCESS;
+ /* compute full volume tensor from the angles */
+ Setup_Box( atof(s_a), atof(s_b), atof(s_c),
+ atof(s_alpha), atof(s_beta), atof(s_gamma),
+ &(system->box) );
+
+ ret = SUCCESS;
+ }
}
}
if ( ferror( geo ) )
{
- return FAILURE;
+ ret = FAILURE;
}
- return FAILURE;
+ return ret;
}
@@ -200,8 +206,10 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
s_y[8] = 0;
strncpy( s_z, line + 46, 8 );
s_z[8] = 0;
+
Make_Point( strtod( s_x, &endptr ), strtod( s_y, &endptr ),
- strtod( s_z, &endptr ), &x );
+ strtod( s_z, &endptr ), &x );
+
Fit_to_Periodic_Box( &(system->box), &x );
}
}
@@ -372,20 +380,6 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
atom->q = 0;
top++;
- // fprintf( stderr, "p%d: %6d%2d x:%8.3f%8.3f%8.3f"
- // "q:%8.3f occ:%s temp:%s seg:%s elmnt:%s\n",
- // system->my_rank,
- // c, system->my_atoms[top].type,
- // system->my_atoms[top].x[0],
- // system->my_atoms[top].x[1],
- // system->my_atoms[top].x[2],
- // system->my_atoms[top].q, occupancy, temp_factor,
- // seg_id, element );
-
-// fprintf( stderr, "atom( %8.3f %8.3f %8.3f )\n",
-// atom->x[0], atom->x[1],
-// atom->x[2] );
-
c++;
}
@@ -542,15 +536,15 @@ char Write_PDB( reax_system* system, list* bonds, simulation_data *data,
}
*/
- free( buffer );
- free( line );
+ sfree( buffer, "Write_PDB::buffer" );
+ sfree( line, "Write_PDB::line" );
return SUCCESS;
}
char Read_BGF( char* bgf_file, reax_system* system, control_params *control,
- simulation_data *data, static_storage *workspace )
+ simulation_data *data, static_storage *workspace )
{
FILE *bgf;
char **tokens;
@@ -738,9 +732,9 @@ char Read_BGF( char* bgf_file, reax_system* system, control_params *control,
&s_gamma[0] );
/* Compute full volume tensor from the angles */
- Setup_Box( atof(s_a), atof(s_b), atof(s_c),
- atof(s_alpha), atof(s_beta), atof(s_gamma),
- &(system->box) );
+ Setup_Box( atof(s_a), atof(s_b), atof(s_c),
+ atof(s_alpha), atof(s_beta), atof(s_gamma),
+ &(system->box) );
}
else if (!strncmp( tokens[0], "CONECT", 6 ))
{
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index d48644f6..3b7a45ce 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -20,7 +20,9 @@
----------------------------------------------------------------------*/
#include "grid.h"
+
#include "reset_utils.h"
+#include "tool_box.h"
#include "vector.h"
@@ -39,9 +41,10 @@ int Estimate_GCell_Population( reax_system* system )
j = (int)(system->atoms[l].x[1] * g->inv_len[1]);
k = (int)(system->atoms[l].x[2] * g->inv_len[2]);
g->top[i][j][k]++;
- // fprintf( stderr, "\tatom%-6d (%8.3f%8.3f%8.3f) --> (%3d%3d%3d)\n",
- // l, system->atoms[l].x[0], system->atoms[l].x[1], system->atoms[l].x[2],
- // i, j, k );
+
+// fprintf( stderr, "\tatom%-6d (%8.3f%8.3f%8.3f) --> (%3d%3d%3d)\n",
+// l, system->atoms[l].x[0], system->atoms[l].x[1], system->atoms[l].x[2],
+// i, j, k );
}
max_atoms = 0;
@@ -59,7 +62,7 @@ int Estimate_GCell_Population( reax_system* system )
}
}
- return MAX(max_atoms * SAFE_ZONE, MIN_GCELL_POPL);
+ return MAX( max_atoms * SAFE_ZONE, MIN_GCELL_POPL );
}
@@ -107,8 +110,8 @@ void Allocate_Space_for_Grid( reax_system *system )
g->mark[i][j][k] = 0;
g->start[i][j][k] = 0;
g->end[i][j][k] = 0;
- g->nbrs[i][j][k] = (ivec*) calloc( g->max_nbrs, sizeof( ivec ) );
- g->nbrs_cp[i][j][k] = (rvec*) calloc( g->max_nbrs, sizeof( rvec ) );
+ g->nbrs[i][j][k] = (ivec*) malloc( g->max_nbrs * sizeof( ivec ) );
+ g->nbrs_cp[i][j][k] = (rvec*) malloc( g->max_nbrs * sizeof( rvec ) );
for ( l = 0; l < g->max_nbrs; ++l )
{
@@ -116,22 +119,28 @@ void Allocate_Space_for_Grid( reax_system *system )
g->nbrs[i][j][k][l][1] = -1;
g->nbrs[i][j][k][l][2] = -1;
- g->nbrs_cp[i][j][k][l][0] = -1;
- g->nbrs_cp[i][j][k][l][1] = -1;
- g->nbrs_cp[i][j][k][l][2] = -1;
+ g->nbrs_cp[i][j][k][l][0] = -1.0;
+ g->nbrs_cp[i][j][k][l][1] = -1.0;
+ g->nbrs_cp[i][j][k][l][2] = -1.0;
}
}
}
}
g->max_atoms = Estimate_GCell_Population( system );
+
for ( i = 0; i < g->ncell[0]; i++ )
{
for ( j = 0; j < g->ncell[1]; j++ )
{
for ( k = 0; k < g->ncell[2]; k++ )
{
- g->atoms[i][j][k] = (int*) calloc( g->max_atoms, sizeof(int) );
+ g->atoms[i][j][k] = (int*) malloc( g->max_atoms * sizeof( int ) );
+
+ for ( l = 0; l < g->max_atoms; ++l )
+ {
+ g->atoms[i][j][k][l] = -1;
+ }
}
}
}
@@ -150,36 +159,36 @@ void Deallocate_Grid_Space( grid *g )
{
for ( k = 0; k < g->ncell[2]; k++ )
{
- free( g->atoms[i][j][k] );
- free( g->nbrs[i][j][k] );
- free( g->nbrs_cp[i][j][k] );
+ sfree( g->atoms[i][j][k], "Deallocate_Grid_Space::g->atoms[i][j][k]" );
+ sfree( g->nbrs[i][j][k], "Deallocate_Grid_Space::g->nbrs[i][j][k]" );
+ sfree( g->nbrs_cp[i][j][k], "Deallocate_Grid_Space::g->nbrs_cp[i][j][k]" );
}
- free( g->atoms[i][j] );
- free( g->top[i][j] );
- free( g->start[i][j] );
- free( g->end[i][j] );
- free( g->mark[i][j] );
- free( g->nbrs[i][j] );
- free( g->nbrs_cp[i][j] );
+ sfree( g->atoms[i][j], "Deallocate_Grid_Space::g->atoms[i][j]" );
+ sfree( g->top[i][j], "Deallocate_Grid_Space::g->top[i][j]" );
+ sfree( g->mark[i][j], "Deallocate_Grid_Space::g->mark[i][j]" );
+ sfree( g->start[i][j], "Deallocate_Grid_Space::g->start[i][j]" );
+ sfree( g->end[i][j], "Deallocate_Grid_Space::g->end[i][j]" );
+ sfree( g->nbrs[i][j], "Deallocate_Grid_Space::g->nbrs[i][j]" );
+ sfree( g->nbrs_cp[i][j], "Deallocate_Grid_Space::g->nbrs_cp[i][j]" );
}
- free( g->atoms[i] );
- free( g->top[i] );
- free( g->start[i] );
- free( g->end[i] );
- free( g->mark[i] );
- free( g->nbrs[i] );
- free( g->nbrs_cp[i] );
+ sfree( g->atoms[i], "Deallocate_Grid_Space::g->atoms[i]" );
+ sfree( g->top[i], "Deallocate_Grid_Space::g->top[i]" );
+ sfree( g->mark[i], "Deallocate_Grid_Space::g->mark[i]" );
+ sfree( g->start[i], "Deallocate_Grid_Space::g->start[i]" );
+ sfree( g->end[i], "Deallocate_Grid_Space::g->end[i]" );
+ sfree( g->nbrs[i], "Deallocate_Grid_Space::g->nbrs[i]" );
+ sfree( g->nbrs_cp[i], "Deallocate_Grid_Space::g->nbrs_cp[i]" );
}
- free( g->atoms );
- free( g->top );
- free( g->start );
- free( g->end );
- free( g->mark );
- free( g->nbrs );
- free( g->nbrs_cp );
+ sfree( g->atoms, "Deallocate_Grid_Space::g->atoms" );
+ sfree( g->top, "Deallocate_Grid_Space::g->top" );
+ sfree( g->mark, "Deallocate_Grid_Space::g->mark" );
+ sfree( g->start, "Deallocate_Grid_Space::g->start" );
+ sfree( g->end, "Deallocate_Grid_Space::g->end" );
+ sfree( g->nbrs, "Deallocate_Grid_Space::g->nbrs" );
+ sfree( g->nbrs_cp, "Deallocate_Grid_Space::g->nbrs_cp" );
}
@@ -425,7 +434,8 @@ void Update_Grid( reax_system* system )
}
}
}
- else /* at least one of ncell has changed */
+ /* at least one of ncell has changed */
+ else
{
Deallocate_Grid_Space( g );
/* update number of grid cells */
@@ -437,6 +447,7 @@ void Update_Grid( reax_system* system )
Allocate_Space_for_Grid( system );
Find_Neighbor_GridCells( g );
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "updated grid: " );
fprintf( stderr, "ncell[%d %d %d] ",
@@ -463,6 +474,7 @@ void Bin_Atoms( reax_system* system, static_storage *workspace )
k = (int)(system->atoms[l].x[2] * g->inv_len[2]);
g->atoms[i][j][k][g->top[i][j][k]] = l;
g->top[i][j][k]++;
+
// fprintf( stderr, "\tatom%-6d (%8.3f%8.3f%8.3f) --> (%3d%3d%3d)\n",
// l, system->atoms[l].x[0], system->atoms[l].x[1], system->atoms[l].x[2],
// i, j, k );
@@ -486,7 +498,8 @@ void Bin_Atoms( reax_system* system, static_storage *workspace )
/* check if current gcell->max_atoms is safe */
if ( max_atoms >= g->max_atoms * SAFE_ZONE )
{
- workspace->realloc.gcell_atoms = MAX(max_atoms * SAFE_ZONE, MIN_GCELL_POPL);
+ workspace->realloc.gcell_atoms = MAX( max_atoms * SAFE_ZONE,
+ MIN_GCELL_POPL );
}
}
@@ -547,33 +560,29 @@ void Free_Storage( static_storage *workspace, control_params * control )
for ( i = 0; i < control->cm_solver_restart + 1; ++i )
{
- free( workspace->v[i] );
+ sfree( workspace->v[i], "Free_Storage::workspace->v[i]" );
}
- free( workspace->v );
+ sfree( workspace->v, "Free_Storage::workspace->v" );
for ( i = 0; i < 3; ++i )
{
- free( workspace->s[i] );
- free( workspace->t[i] );
+ sfree( workspace->s[i], "Free_Storage::workspace->s[i]" );
+ sfree( workspace->t[i], "Free_Storage::workspace->t[i]" );
}
- free( workspace->s );
- free( workspace->t );
+ sfree( workspace->s, "Free_Storage::workspace->s" );
+ sfree( workspace->t, "Free_Storage::workspace->t" );
- free( workspace->orig_id );
+ sfree( workspace->orig_id, "Free_Storage::workspace->orig_id" );
}
void Assign_New_Storage( static_storage *workspace,
- real **v, real **s, real **t,
- int *orig_id, rvec *f_old )
+ real **v, real **s, real **t, int *orig_id, rvec *f_old )
{
workspace->v = v;
-
workspace->s = s;
workspace->t = t;
-
workspace->orig_id = orig_id;
-
workspace->f_old = f_old;
}
@@ -581,14 +590,14 @@ void Assign_New_Storage( static_storage *workspace,
void Cluster_Atoms( reax_system *system, static_storage *workspace,
control_params *control )
{
- int i, j, k, l, top, old_id, num_H;
+ int i, j, k, l, top, old_id, num_H;
reax_atom *old_atom;
- grid *g;
+ grid *g;
reax_atom *new_atoms;
- int *orig_id ;
- real **v;
- real **s, **t;
- rvec *f_old;
+ int *orig_id ;
+ real **v;
+ real **s, **t;
+ rvec *f_old;
num_H = 0;
g = &( system->g );
@@ -638,7 +647,7 @@ void Cluster_Atoms( reax_system *system, static_storage *workspace,
}
- free( system->atoms );
+ sfree( system->atoms, "Cluster_Atoms::system->atoms" );
Free_Storage( workspace, control );
system->atoms = new_atoms;
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 26121468..44e86371 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -20,19 +20,18 @@
----------------------------------------------------------------------*/
#include "init_md.h"
+
#include "allocate.h"
#include "box.h"
#include "forces.h"
#include "grid.h"
#include "integrate.h"
-#include "lin_alg.h"
#include "neighbors.h"
#include "list.h"
#include "lookup.h"
-#include "print_utils.h"
#include "reset_utils.h"
#include "system_props.h"
-#include "traj.h"
+#include "tool_box.h"
#include "vector.h"
@@ -639,8 +638,8 @@ void Init_Lists( reax_system *system, control_params *control,
}
#endif
- free( hb_top );
- free( bond_top );
+ sfree( hb_top, "Init_Lists::hb_top" );
+ sfree( bond_top, "Init_Lists::bond_top" );
}
@@ -936,7 +935,7 @@ void Finalize_System( reax_system *system, control_params *control,
Finalize_Grid( system );
- free( reax->gp.l );
+ sfree( reax->gp.l, "Finalize_System::reax->gp.l" );
for ( i = 0; i < reax->num_atom_types; i++ )
{
@@ -944,27 +943,27 @@ void Finalize_System( reax_system *system, control_params *control,
{
for ( k = 0; k < reax->num_atom_types; k++ )
{
- free( reax->fbp[i][j][k] );
+ sfree( reax->fbp[i][j][k], "Finalize_System::reax->fbp[i][j][k]" );
}
- free( reax->thbp[i][j] );
- free( reax->hbp[i][j] );
- free( reax->fbp[i][j] );
+ sfree( reax->thbp[i][j], "Finalize_System::reax->thbp[i][j]" );
+ sfree( reax->hbp[i][j], "Finalize_System::reax->hbp[i][j]" );
+ sfree( reax->fbp[i][j], "Finalize_System::reax->fbp[i][j]" );
}
- free( reax->tbp[i] );
- free( reax->thbp[i] );
- free( reax->hbp[i] );
- free( reax->fbp[i] );
+ sfree( reax->tbp[i], "Finalize_System::reax->tbp[i]" );
+ sfree( reax->thbp[i], "Finalize_System::reax->thbp[i]" );
+ sfree( reax->hbp[i], "Finalize_System::reax->hbp[i]" );
+ sfree( reax->fbp[i], "Finalize_System::reax->fbp[i]" );
}
- free( reax->sbp );
- free( reax->tbp );
- free( reax->thbp );
- free( reax->hbp );
- free( reax->fbp );
+ sfree( reax->sbp, "Finalize_System::reax->sbp" );
+ sfree( reax->tbp, "Finalize_System::reax->tbp" );
+ sfree( reax->thbp, "Finalize_System::reax->thbp" );
+ sfree( reax->hbp, "Finalize_System::reax->hbp" );
+ sfree( reax->fbp, "Finalize_System::reax->fbp" );
- free( system->atoms );
+ sfree( system->atoms, "Finalize_System::system->atoms" );
}
@@ -979,23 +978,23 @@ void Finalize_Workspace( reax_system *system, control_params *control,
{
int i;
- free( workspace->hbond_index );
- free( workspace->total_bond_order );
- free( workspace->Deltap );
- free( workspace->Deltap_boc );
- free( workspace->dDeltap_self );
- free( workspace->Delta );
- free( workspace->Delta_lp );
- free( workspace->Delta_lp_temp );
- free( workspace->dDelta_lp );
- free( workspace->dDelta_lp_temp );
- free( workspace->Delta_e );
- free( workspace->Delta_boc );
- free( workspace->nlp );
- free( workspace->nlp_temp );
- free( workspace->Clp );
- free( workspace->CdDelta );
- free( workspace->vlpex );
+ sfree( workspace->hbond_index, "Finalize_Workspace::workspace->hbond_index" );
+ sfree( workspace->total_bond_order, "Finalize_Workspace::workspace->total_bond_order" );
+ sfree( workspace->Deltap, "Finalize_Workspace::workspace->Deltap" );
+ sfree( workspace->Deltap_boc, "Finalize_Workspace::workspace->Deltap_boc" );
+ sfree( workspace->dDeltap_self, "Finalize_Workspace::workspace->dDeltap_self" );
+ sfree( workspace->Delta, "Finalize_Workspace::workspace->Delta" );
+ sfree( workspace->Delta_lp, "Finalize_Workspace::workspace->Delta_lp" );
+ sfree( workspace->Delta_lp_temp, "Finalize_Workspace::workspace->Delta_lp_temp" );
+ sfree( workspace->dDelta_lp, "Finalize_Workspace::workspace->dDelta_lp" );
+ sfree( workspace->dDelta_lp_temp, "Finalize_Workspace::workspace->dDelta_lp_temp" );
+ sfree( workspace->Delta_e, "Finalize_Workspace::workspace->Delta_e" );
+ sfree( workspace->Delta_boc, "Finalize_Workspace::workspace->Delta_boc" );
+ sfree( workspace->nlp, "Finalize_Workspace::workspace->nlp" );
+ sfree( workspace->nlp_temp, "Finalize_Workspace::workspace->nlp_temp" );
+ sfree( workspace->Clp, "Finalize_Workspace::workspace->Clp" );
+ sfree( workspace->CdDelta, "Finalize_Workspace::workspace->CdDelta" );
+ sfree( workspace->vlpex, "Finalize_Workspace::workspace->vlpex" );
Deallocate_Matrix( workspace->H );
Deallocate_Matrix( workspace->H_sp );
@@ -1009,26 +1008,29 @@ void Finalize_Workspace( reax_system *system, control_params *control,
for ( i = 0; i < 5; ++i )
{
- free( workspace->s[i] );
- free( workspace->t[i] );
+ sfree( workspace->s[i], "Finalize_Workspace::workspace->s[i]" );
+ sfree( workspace->t[i], "Finalize_Workspace::workspace->t[i]" );
}
- free( workspace->Hdia_inv );
+ if ( control->cm_solver_pre_comp_type == DIAG_PC )
+ {
+ sfree( workspace->Hdia_inv, "Finalize_Workspace::workspace->Hdia_inv" );
+ }
if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
control->cm_solver_pre_comp_type == ILUT_PAR_PC )
{
- free( workspace->droptol );
+ sfree( workspace->droptol, "Finalize_Workspace::workspace->droptol" );
}
//TODO: check if unused
- //free( workspace->w );
+ //sfree( workspace->w, "Finalize_Workspace::workspace->w" );
//TODO: check if unused
- free( workspace->b );
- free( workspace->b_s );
- free( workspace->b_t );
- free( workspace->b_prc );
- free( workspace->b_prm );
- free( workspace->s );
- free( workspace->t );
+ sfree( workspace->b, "Finalize_Workspace::workspace->b" );
+ sfree( workspace->b_s, "Finalize_Workspace::workspace->b_s" );
+ sfree( workspace->b_t, "Finalize_Workspace::workspace->b_t" );
+ sfree( workspace->b_prc, "Finalize_Workspace::workspace->b_prc" );
+ sfree( workspace->b_prm, "Finalize_Workspace::workspace->b_prm" );
+ sfree( workspace->s, "Finalize_Workspace::workspace->s" );
+ sfree( workspace->t, "Finalize_Workspace::workspace->t" );
switch ( control->cm_solver_type )
{
@@ -1037,38 +1039,38 @@ void Finalize_Workspace( reax_system *system, control_params *control,
case GMRES_H_S:
for ( i = 0; i < control->cm_solver_restart + 1; ++i )
{
- free( workspace->h[i] );
- free( workspace->rn[i] );
- free( workspace->v[i] );
+ sfree( workspace->h[i], "Finalize_Workspace::workspace->h[i]" );
+ sfree( workspace->rn[i], "Finalize_Workspace::workspace->rn[i]" );
+ sfree( workspace->v[i], "Finalize_Workspace::workspace->v[i]" );
}
- free( workspace->y );
- free( workspace->z );
- free( workspace->g );
- free( workspace->h );
- free( workspace->hs );
- free( workspace->hc );
- free( workspace->rn );
- free( workspace->v );
-
- free( workspace->r );
- free( workspace->d );
- free( workspace->q );
- free( workspace->p );
+ sfree( workspace->y, "Finalize_Workspace::workspace->y" );
+ sfree( workspace->z, "Finalize_Workspace::workspace->z" );
+ sfree( workspace->g, "Finalize_Workspace::workspace->g" );
+ sfree( workspace->h, "Finalize_Workspace::workspace->h" );
+ sfree( workspace->hs, "Finalize_Workspace::workspace->hs" );
+ sfree( workspace->hc, "Finalize_Workspace::workspace->hc" );
+ sfree( workspace->rn, "Finalize_Workspace::workspace->rn" );
+ sfree( workspace->v, "Finalize_Workspace::workspace->v" );
+
+ sfree( workspace->r, "Finalize_Workspace::workspace->r" );
+ sfree( workspace->d, "Finalize_Workspace::workspace->d" );
+ sfree( workspace->q, "Finalize_Workspace::workspace->q" );
+ sfree( workspace->p, "Finalize_Workspace::workspace->p" );
break;
/* CG storage */
case CG_S:
- free( workspace->r );
- free( workspace->d );
- free( workspace->q );
- free( workspace->p );
+ sfree( workspace->r, "Finalize_Workspace::workspace->r" );
+ sfree( workspace->d, "Finalize_Workspace::workspace->d" );
+ sfree( workspace->q, "Finalize_Workspace::workspace->q" );
+ sfree( workspace->p, "Finalize_Workspace::workspace->p" );
break;
case SDM_S:
- free( workspace->r );
- free( workspace->d );
- free( workspace->q );
+ sfree( workspace->r, "Finalize_Workspace::workspace->r" );
+ sfree( workspace->d, "Finalize_Workspace::workspace->d" );
+ sfree( workspace->q, "Finalize_Workspace::workspace->q" );
break;
default:
@@ -1078,63 +1080,56 @@ void Finalize_Workspace( reax_system *system, control_params *control,
}
/* integrator storage */
- free( workspace->a );
- free( workspace->f_old );
- free( workspace->v_const );
+ sfree( workspace->a, "Finalize_Workspace::workspace->a" );
+ sfree( workspace->f_old, "Finalize_Workspace::workspace->f_old" );
+ sfree( workspace->v_const, "Finalize_Workspace::workspace->v_const" );
#ifdef _OPENMP
- free( workspace->f_local );
+ sfree( workspace->f_local, "Finalize_Workspace::workspace->f_local" );
#endif
/* storage for analysis */
if ( control->molec_anal || control->diffusion_coef )
{
- free( workspace->mark );
- free( workspace->old_mark );
- }
- else
- {
- free( workspace->mark );
+ sfree( workspace->mark, "Finalize_Workspace::workspace->mark" );
+ sfree( workspace->old_mark, "Finalize_Workspace::workspace->old_mark" );
}
if ( control->diffusion_coef )
{
- free( workspace->x_old );
- }
- else
- {
- free( workspace->x_old );
+ sfree( workspace->x_old, "Finalize_Workspace::workspace->x_old" );
}
- free( workspace->orig_id );
+ sfree( workspace->orig_id, "Finalize_Workspace::workspace->orig_id" );
/* space for keeping restriction info, if any */
if ( control->restrict_bonds )
{
for ( i = 0; i < system->N; ++i )
{
- free( workspace->restricted_list[i] );
+ sfree( workspace->restricted_list[i],
+ "Finalize_Workspace::workspace->restricted_list[i]" );
}
- free( workspace->restricted );
- free( workspace->restricted_list );
+ sfree( workspace->restricted, "Finalize_Workspace::workspace->restricted" );
+ sfree( workspace->restricted_list, "Finalize_Workspace::workspace->restricted_list" );
}
#ifdef TEST_FORCES
- free( workspace->dDelta );
- free( workspace->f_ele );
- free( workspace->f_vdw );
- free( workspace->f_bo );
- free( workspace->f_be );
- free( workspace->f_lp );
- free( workspace->f_ov );
- free( workspace->f_un );
- free( workspace->f_ang );
- free( workspace->f_coa );
- free( workspace->f_pen );
- free( workspace->f_hb );
- free( workspace->f_tor );
- free( workspace->f_con );
+ sfree( workspace->dDelta, "Finalize_Workspace::workspace->dDelta" );
+ sfree( workspace->f_ele, "Finalize_Workspace::workspace->f_ele" );
+ sfree( workspace->f_vdw, "Finalize_Workspace::workspace->f_vdw" );
+ sfree( workspace->f_bo, "Finalize_Workspace::workspace->f_bo" );
+ sfree( workspace->f_be, "Finalize_Workspace::workspace->f_be" );
+ sfree( workspace->f_lp, "Finalize_Workspace::workspace->f_lp" );
+ sfree( workspace->f_ov, "Finalize_Workspace::workspace->f_ov" );
+ sfree( workspace->f_un, "Finalize_Workspace::workspace->f_un" );
+ sfree( workspace->f_ang, "Finalize_Workspace::workspace->f_ang" );
+ sfree( workspace->f_coa, "Finalize_Workspace::workspace->f_coa" );
+ sfree( workspace->f_pen, "Finalize_Workspace::workspace->f_pen" );
+ sfree( workspace->f_hb, "Finalize_Workspace::workspace->f_hb" );
+ sfree( workspace->f_tor, "Finalize_Workspace::workspace->f_tor" );
+ sfree( workspace->f_con, "Finalize_Workspace::workspace->f_con" );
#endif
}
@@ -1159,19 +1154,31 @@ void Finalize_Out_Controls( reax_system *system, control_params *control,
/* close trajectory file */
if ( out_control->write_steps > 0 )
{
- fclose( out_control->trj );
+ if ( out_control->trj )
+ {
+ fclose( out_control->trj );
+ }
}
if ( out_control->energy_update_freq > 0 )
{
/* close out file */
- fclose( out_control->out );
+ if ( out_control->out )
+ {
+ fclose( out_control->out );
+ }
/* close potentials file */
- fclose( out_control->pot );
+ if ( out_control->pot )
+ {
+ fclose( out_control->pot );
+ }
/* close log file */
- fclose( out_control->log );
+ if ( out_control->log )
+ {
+ fclose( out_control->log );
+ }
}
/* close pressure file */
@@ -1179,7 +1186,10 @@ void Finalize_Out_Controls( reax_system *system, control_params *control,
control->ensemble == iNPT ||
control->ensemble == sNPT )
{
- fclose( out_control->prs );
+ if ( out_control->prs )
+ {
+ fclose( out_control->prs );
+ }
}
/* close molecular analysis file */
@@ -1203,76 +1213,145 @@ void Finalize_Out_Controls( reax_system *system, control_params *control,
#ifdef TEST_ENERGY
/* close bond energy file */
- fclose( out_control->ebond );
+ if ( out_control->ebond )
+ {
+ fclose( out_control->ebond );
+ }
/* close lone-pair energy file */
- fclose( out_control->elp );
+ if ( out_control->help )
+ {
+ fclose( out_control->elp );
+ }
/* close overcoordination energy file */
- fclose( out_control->eov );
+ if ( out_control->eov )
+ {
+ fclose( out_control->eov );
+ }
/* close undercoordination energy file */
- fclose( out_control->eun );
+ if ( out_control->eun )
+ {
+ fclose( out_control->eun );
+ }
/* close angle energy file */
- fclose( out_control->eval );
+ if ( out_control->eval )
+ {
+ fclose( out_control->eval );
+ }
/* close penalty energy file */
- fclose( out_control->epen );
+ if ( out_control->epen )
+ {
+ fclose( out_control->epen );
+ }
/* close coalition energy file */
- fclose( out_control->ecoa );
+ if ( out_control->ecoa )
+ {
+ fclose( out_control->ecoa );
+ }
/* close hydrogen bond energy file */
- fclose( out_control->ehb );
+ if ( out_control->ehb )
+ {
+ fclose( out_control->ehb );
+ }
/* close torsion energy file */
- fclose( out_control->etor );
+ if ( out_control->etor )
+ {
+ fclose( out_control->etor );
+ }
/* close conjugation energy file */
- fclose( out_control->econ );
+ if ( out_control->econ )
+ {
+ fclose( out_control->econ );
+ }
/* close vdWaals energy file */
- fclose( out_control->evdw );
+ if ( out_control->evdw )
+ {
+ fclose( out_control->evdw );
+ }
/* close coulomb energy file */
- fclose( out_control->ecou );
+ if ( out_control->ecou )
+ {
+ fclose( out_control->ecou );
+ }
#endif
#ifdef TEST_FORCES
/* close bond orders file */
- fclose( out_control->fbo );
+ if ( out_control->fbo )
+ {
+ fclose( out_control->fbo );
+ }
/* close bond orders derivatives file */
- fclose( out_control->fdbo );
+ if ( out_control->fdbo )
+ {
+ fclose( out_control->fdbo );
+ }
/* close bond forces file */
- fclose( out_control->fbond );
+ if ( out_control->fbond )
+ {
+ fclose( out_control->fbond );
+ }
/* close lone-pair forces file */
- fclose( out_control->flp );
+ if ( out_control->flp )
+ {
+ fclose( out_control->flp );
+ }
/* close overcoordination forces file */
- fclose( out_control->fatom );
+ if ( out_control->fatom )
+ {
+ fclose( out_control->fatom );
+ }
/* close angle forces file */
- fclose( out_control->f3body );
+ if ( out_control->f3body )
+ {
+ fclose( out_control->f3body );
+ }
/* close hydrogen bond forces file */
- fclose( out_control->fhb );
+ if ( out_control->fhb )
+ {
+ fclose( out_control->fhb );
+ }
/* close torsion forces file */
- fclose( out_control->f4body );
+ if ( out_control->f4body )
+ {
+ fclose( out_control->f4body );
+ }
/* close nonbonded forces file */
- fclose( out_control->fnonb );
+ if ( out_control->fnonb )
+ {
+ fclose( out_control->fnonb );
+ }
/* close total force file */
- fclose( out_control->ftot );
+ if ( out_control->ftot )
+ {
+ fclose( out_control->ftot );
+ }
/* close coulomb forces file */
- fclose( out_control->ftot2 );
+ if ( out_control->ftot2 )
+ {
+ fclose( out_control->ftot2 );
+ }
#endif
}
@@ -1283,7 +1362,7 @@ void Finalize( reax_system *system, control_params *control,
{
if ( control->tabulate )
{
-// Finalize_LR_Lookup_Table( system, control );
+ Finalize_LR_Lookup_Table( system, control );
}
Finalize_Out_Controls( system, control, workspace, out_control );
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 208791b8..6a48b042 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -22,8 +22,6 @@
#include "lin_alg.h"
#include "allocate.h"
-#include "list.h"
-#include "print_utils.h"
#include "tool_box.h"
#include "vector.h"
@@ -81,11 +79,10 @@ static void Sparse_MatVec( const sparse_matrix * const A,
Vector_MakeZero( b, n );
#ifdef _OPENMP
- tid = omp_get_thread_num();
+ tid = omp_get_thread_num( );
#pragma omp master
{
-
/* keep b_local for program duration to avoid allocate/free
* overhead per Sparse_MatVec call*/
if ( b_local == NULL )
@@ -101,8 +98,8 @@ static void Sparse_MatVec( const sparse_matrix * const A,
Vector_MakeZero( (real * const)b_local, omp_get_num_threads() * n );
-#endif
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < n; ++i )
{
si = A->start[i];
@@ -128,6 +125,7 @@ static void Sparse_MatVec( const sparse_matrix * const A,
b[i] += A->val[k] * x[i];
#endif
}
+
#ifdef _OPENMP
#pragma omp for schedule(static)
for ( i = 0; i < n; ++i )
@@ -138,7 +136,6 @@ static void Sparse_MatVec( const sparse_matrix * const A,
}
}
#endif
-
}
@@ -186,7 +183,7 @@ void Transpose( const sparse_matrix * const A, sparse_matrix const *A_t )
}
}
- free( A_t_top );
+ sfree( A_t_top, "Transpose::A_t_top" );
}
@@ -226,7 +223,9 @@ static void diag_pre_app( const real * const Hdia_inv, const real * const y,
{
unsigned int i;
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < N; ++i )
{
x[i] = y[i] * Hdia_inv[i];
@@ -251,7 +250,9 @@ void tri_solve( const sparse_matrix * const LU, const real * const y,
int i, pj, j, si, ei;
real val;
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
if ( tri == LOWER )
{
@@ -286,6 +287,10 @@ void tri_solve( const sparse_matrix * const LU, const real * const y,
}
}
}
+
+#ifdef _OPENMP
+ #pragma omp barrier
+#endif
}
@@ -307,7 +312,9 @@ void tri_solve_level_sched( const sparse_matrix * const LU,
{
int i, j, pj, local_row, local_level;
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
if ( tri == LOWER )
{
@@ -407,14 +414,18 @@ void tri_solve_level_sched( const sparse_matrix * const LU,
}
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
/* perform substitutions by level */
if ( tri == LOWER )
{
for ( i = 0; i < levels; ++i )
{
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( j = level_rows_cnt[i]; j < level_rows_cnt[i + 1]; ++j )
{
local_row = level_rows[j];
@@ -432,7 +443,9 @@ void tri_solve_level_sched( const sparse_matrix * const LU,
{
for ( i = 0; i < levels; ++i )
{
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( j = level_rows_cnt[i]; j < level_rows_cnt[i + 1]; ++j )
{
local_row = level_rows[j];
@@ -447,7 +460,9 @@ void tri_solve_level_sched( const sparse_matrix * const LU,
}
}
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
/* save level info for re-use if performing repeated triangular solves via preconditioning */
if ( tri == LOWER )
@@ -466,7 +481,9 @@ void tri_solve_level_sched( const sparse_matrix * const LU,
}
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
}
@@ -527,7 +544,9 @@ static void compute_H_full( const sparse_matrix * const H )
*/
void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
{
+#ifdef _OPENMP
#pragma omp parallel
+#endif
{
#define MAX_COLOR (500)
int i, pj, v;
@@ -542,7 +561,9 @@ void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
num_thread = 1;
#endif
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
memset( color, 0, sizeof(unsigned int) * A->n );
recolor_cnt = A->n;
@@ -552,7 +573,9 @@ void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
* for which coloring is to be used for */
if ( tri == LOWER )
{
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < A->n; ++i )
{
to_color[i] = i;
@@ -560,7 +583,9 @@ void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
}
else
{
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < A->n; ++i )
{
to_color[i] = A->n - 1 - i;
@@ -574,14 +599,18 @@ void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
exit( INSUFFICIENT_MEMORY );
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
while ( recolor_cnt > 0 )
{
memset( fb_color, -1, sizeof(int) * MAX_COLOR );
/* color vertices */
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < recolor_cnt; ++i )
{
v = to_color[i];
@@ -607,16 +636,22 @@ void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
temp = recolor_cnt;
recolor_cnt_local = 0;
+#ifdef _OPENMP
#pragma omp barrier
+#endif
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
recolor_cnt = 0;
}
-
+
+#ifdef _OPENMP
#pragma omp barrier
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < temp; ++i )
{
v = to_color[i];
@@ -634,21 +669,25 @@ void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
}
/* count thread-local conflicts and compute offsets for copying into shared buffer */
- conflict_cnt[tid + 1] = recolor_cnt_local;
+ conflict_cnt[tid + 1] = recolor_cnt_local;
+#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
+#endif
{
conflict_cnt[0] = 0;
for ( i = 1; i < num_thread + 1; ++i )
{
- conflict_cnt[i] += conflict_cnt[i - 1];
+ conflict_cnt[i] += conflict_cnt[i - 1];
}
recolor_cnt = conflict_cnt[num_thread];
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
/* copy thread-local conflicts into shared buffer */
for ( i = 0; i < recolor_cnt_local; ++i )
@@ -657,22 +696,38 @@ void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
color[conflict_local[i]] = 0;
}
+#ifdef _OPENMP
#pragma omp barrier
- #pragma omp master
+ #pragma omp single
+#endif
{
temp_ptr = to_color;
to_color = conflict;
conflict = temp_ptr;
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
}
- free( conflict_local );
- free( fb_color );
+ sfree( conflict_local, "graph_coloring::conflict_local" );
+ sfree( fb_color, "graph_coloring::fb_color" );
+//#if defined(DEBUG)
+//#ifdef _OPENMP
+// #pragma omp master
+//#endif
+// {
+// for ( i = 0; i < A->n; ++i )
+// printf("Vertex: %5d, Color: %5d\n", i, color[i] );
+// }
+//#endif
+
+#ifdef _OPENMP
#pragma omp barrier
+#endif
}
}
@@ -729,7 +784,9 @@ static void permute_vector( real * const x, const unsigned int n, const int inve
{
unsigned int i;
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
if ( x_p == NULL )
{
@@ -750,20 +807,26 @@ static void permute_vector( real * const x, const unsigned int n, const int inve
}
}
+#ifdef _OPENMP
#pragma omp barrier
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < n; ++i )
{
x_p[i] = x[mapping[i]];
}
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
memcpy( x, x_p, sizeof(real) * n );
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
}
@@ -974,7 +1037,9 @@ void jacobi_iter( const sparse_matrix * const R, const real * const Dinv,
iter = 0;
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
if ( Dinv_b == NULL )
{
@@ -1002,12 +1067,16 @@ void jacobi_iter( const sparse_matrix * const R, const real * const Dinv,
}
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
Vector_MakeZero( rp, R->n );
/* precompute and cache, as invariant in loop below */
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < R->n; ++i )
{
Dinv_b[i] = Dinv[i] * b[i];
@@ -1016,7 +1085,9 @@ void jacobi_iter( const sparse_matrix * const R, const real * const Dinv,
do
{
// x_{k+1} = G*x_{k} + Dinv*b;
+#ifdef _OPENMP
#pragma omp for schedule(guided)
+#endif
for ( i = 0; i < R->n; ++i )
{
if (tri == LOWER)
@@ -1042,14 +1113,18 @@ void jacobi_iter( const sparse_matrix * const R, const real * const Dinv,
rp2[i] += Dinv_b[i];
}
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
rp3 = rp;
rp = rp2;
rp2 = rp3;
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
++iter;
}
@@ -1130,12 +1205,16 @@ static void apply_preconditioner( const static_storage * const workspace, const
case ICHOLT_PC:
case ILU_PAR_PC:
case ILUT_PAR_PC:
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
memcpy( y_p, y, sizeof(real) * workspace->H->n );
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
permute_vector( y_p, workspace->H->n, FALSE, LOWER );
tri_solve_level_sched( workspace->L, y_p, x, workspace->L->n, LOWER, fresh_pre );
@@ -1158,7 +1237,9 @@ static void apply_preconditioner( const static_storage * const workspace, const
case ICHOLT_PC:
case ILU_PAR_PC:
case ILUT_PAR_PC:
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
if ( Dinv_L == NULL )
{
@@ -1170,12 +1251,16 @@ static void apply_preconditioner( const static_storage * const workspace, const
}
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
/* construct D^{-1}_L */
if ( fresh_pre == TRUE )
{
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < workspace->L->n; ++i )
{
si = workspace->L->start[i + 1] - 1;
@@ -1185,7 +1270,9 @@ static void apply_preconditioner( const static_storage * const workspace, const
jacobi_iter( workspace->L, Dinv_L, y, x, LOWER, control->cm_solver_pre_app_jacobi_iters );
- #pragma omp master
+#ifdef _OPENMP
+ #pragma omp single
+#endif
{
if ( Dinv_U == NULL )
{
@@ -1197,12 +1284,16 @@ static void apply_preconditioner( const static_storage * const workspace, const
}
}
+#ifdef _OPENMP
#pragma omp barrier
+#endif
/* construct D^{-1}_U */
if ( fresh_pre == TRUE )
{
+#ifdef _OPENMP
#pragma omp for schedule(static)
+#endif
for ( i = 0; i < workspace->U->n; ++i )
{
si = workspace->U->start[i];
@@ -1238,9 +1329,14 @@ int GMRES( const static_storage * const workspace, const control_params * const
N = H->n;
+#ifdef _OPENMP
#pragma omp parallel default(none) private(i, j, k, itr, bnorm, ret_temp) \
shared(N, cc, tmp1, tmp2, temp, time_start, g_itr, g_j, stderr)
+#endif
{
+ j = 0;
+ itr = 0;
+
#pragma omp master
{
time_start = Get_Time( );
@@ -1253,7 +1349,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
if ( control->cm_solver_pre_comp_type == DIAG_PC )
{
- /* apply preconditioner to RHS */
+ /* apply preconditioner to residual */
#pragma omp master
{
time_start = Get_Time( );
@@ -1340,6 +1436,9 @@ int GMRES( const static_storage * const workspace, const control_params * const
{
workspace->g[0] = ret_temp;
}
+
+ #pragma omp barrier
+
Vector_Scale( workspace->v[0], 1. / workspace->g[0], workspace->v[0], N );
#pragma omp master
{
@@ -1355,6 +1454,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
time_start = Get_Time( );
}
Sparse_MatVec( H, workspace->v[j], workspace->v[j + 1] );
+
#pragma omp master
{
data->timing.cm_solver_spmv += Get_Timing_Info( time_start );
@@ -1365,6 +1465,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
time_start = Get_Time( );
}
apply_preconditioner( workspace, control, workspace->v[j + 1], workspace->v[j + 1], FALSE );
+
#pragma omp master
{
data->timing.cm_solver_pre_app += Get_Timing_Info( time_start );
@@ -1380,7 +1481,11 @@ int GMRES( const static_storage * const workspace, const control_params * const
for ( i = 0; i <= j; i++ )
{
workspace->h[i][j] = Dot( workspace->v[i], workspace->v[j + 1], N );
+
+ #pragma omp barrier
+
Vector_Add( workspace->v[j + 1], -workspace->h[i][j], workspace->v[i], N );
+
}
#pragma omp master
{
@@ -1394,9 +1499,12 @@ int GMRES( const static_storage * const workspace, const control_params * const
#pragma omp master
{
time_start = Get_Time( );
+ }
+ #pragma omp single
+ {
for ( i = 0; i < j - 1; i++ )
{
- workspace->h[i][j] = 0;
+ workspace->h[i][j] = 0.0;
}
}
@@ -1407,6 +1515,9 @@ int GMRES( const static_storage * const workspace, const control_params * const
{
workspace->h[i][j] = ret_temp;
}
+
+ #pragma omp barrier
+
Vector_Add( workspace->v[j + 1], -workspace->h[i][j], workspace->v[i], N );
}
#pragma omp master
@@ -1424,15 +1535,16 @@ int GMRES( const static_storage * const workspace, const control_params * const
{
workspace->h[j + 1][j] = ret_temp;
}
+
+ #pragma omp barrier
+
Vector_Scale( workspace->v[j + 1],
- 1. / workspace->h[j + 1][j], workspace->v[j + 1], N );
+ 1.0 / workspace->h[j + 1][j], workspace->v[j + 1], N );
+
#pragma omp master
{
data->timing.cm_solver_vector_ops += Get_Timing_Info( time_start );
}
-#if defined(DEBUG)
- fprintf( stderr, "%d-%d: orthogonalization completed.\n", itr, j );
-#endif
#pragma omp master
{
@@ -1487,22 +1599,18 @@ int GMRES( const static_storage * const workspace, const control_params * const
tmp2 = -workspace->hs[j] * workspace->g[j];
workspace->g[j] = tmp1;
workspace->g[j + 1] = tmp2;
+
data->timing.cm_solver_orthog += Get_Timing_Info( time_start );
}
#pragma omp barrier
-
- //fprintf( stderr, "h: " );
- //for( i = 0; i <= j+1; ++i )
- //fprintf( stderr, "%.6f ", workspace->h[i][j] );
- //fprintf( stderr, "\n" );
- //fprintf( stderr, "res: %.15e\n", workspace->g[j+1] );
}
/* solve Hy = g: H is now upper-triangular, do back-substitution */
#pragma omp master
{
time_start = Get_Time( );
+
for ( i = j - 1; i >= 0; i-- )
{
temp = workspace->g[i];
@@ -1513,18 +1621,22 @@ int GMRES( const static_storage * const workspace, const control_params * const
workspace->y[i] = temp / workspace->h[i][i];
}
+
data->timing.cm_solver_tri_solve += Get_Timing_Info( time_start );
/* update x = x_0 + Vy */
time_start = Get_Time( );
}
+
Vector_MakeZero( workspace->p, N );
+
for ( i = 0; i < j; i++ )
{
Vector_Add( workspace->p, workspace->y[i], workspace->v[i], N );
}
Vector_Add( x, 1., workspace->p, N );
+
#pragma omp master
{
data->timing.cm_solver_vector_ops += Get_Timing_Info( time_start );
@@ -1559,11 +1671,12 @@ int GMRES_HouseHolder( const static_storage * const workspace,
const sparse_matrix * const H, const real * const b, real tol,
real * const x, const int fresh_pre )
{
- int i, j, k, itr, N;
+ int i, j, k, itr, N;
real cc, tmp1, tmp2, temp, bnorm;
real v[10000], z[control->cm_solver_restart + 2][10000], w[control->cm_solver_restart + 2];
real u[control->cm_solver_restart + 2][10000];
+ j = 0;
N = H->n;
bnorm = Norm( b, N );
@@ -1597,7 +1710,7 @@ int GMRES_HouseHolder( const static_storage * const workspace,
// fprintf( stderr, "\n\n%12.6f\n", w[0] );
/* GMRES inner-loop */
- for ( j = 0; j < control->cm_solver_restart && fabs( w[j] ) / bnorm > tol; j++ )
+ for ( j = 0; j < control->cm_solver_restart && FABS( w[j] ) / bnorm > tol; j++ )
{
/* compute v_j */
Vector_Scale( z[j], -2 * u[j][j], u[j], N );
@@ -1653,7 +1766,7 @@ int GMRES_HouseHolder( const static_storage * const workspace,
}
/* apply the new Givens rotation to H and right-hand side */
- if ( fabs(v[j + 1]) >= ALMOST_ZERO )
+ if ( FABS(v[j + 1]) >= ALMOST_ZERO )
{
cc = SQRT( SQR( v[j] ) + SQR( v[j + 1] ) );
workspace->hc[j] = v[j] / cc;
@@ -1730,7 +1843,7 @@ int GMRES_HouseHolder( const static_storage * const workspace,
}
/* stopping condition */
- if ( fabs( w[j] ) / bnorm <= tol )
+ if ( FABS( w[j] ) / bnorm <= tol )
{
break;
}
@@ -1906,7 +2019,7 @@ real condest( const sparse_matrix * const L, const sparse_matrix * const U )
}
- free( e );
+ sfree( e, "condest::e" );
return c;
}
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index e02b22fb..d31077c4 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -21,6 +21,8 @@
#include "list.h"
+#include "tool_box.h"
+
char Make_List( int n, int num_intrs, int type, list* l )
{
@@ -104,11 +106,11 @@ void Delete_List( int type, list* l )
{
if ( l->index != NULL )
{
- free( l->index );
+ sfree( l->index, "Delete_List::l->index" );
}
if ( l->end_index != NULL )
{
- free( l->end_index );
+ sfree( l->end_index, "Delete_List::l->end_index" );
}
switch ( type )
@@ -116,54 +118,55 @@ void Delete_List( int type, list* l )
case TYP_VOID:
if ( l->select.v != NULL )
{
- free( l->select.v );
+ sfree( l->select.v, "Delete_List::l->select.v" );
}
break;
case TYP_THREE_BODY:
if ( l->select.three_body_list != NULL )
{
- free( l->select.three_body_list );
+ sfree( l->select.three_body_list, "Delete_List::l->select.three_body_list" );
}
break;
case TYP_BOND:
if ( l->select.bond_list != NULL )
{
- free( l->select.bond_list );
+ sfree( l->select.bond_list, "Delete_List::l->select.bond_list" );
}
break;
case TYP_DBO:
if ( l->select.dbo_list != NULL )
{
- free( l->select.dbo_list );
+ sfree( l->select.dbo_list, "Delete_List::l->select.dbo_list" );
}
break;
case TYP_DDELTA:
if ( l->select.dDelta_list != NULL )
{
- free( l->select.dDelta_list );
+ sfree( l->select.dDelta_list, "Delete_List::l->select.dDelta_list" );
}
break;
case TYP_FAR_NEIGHBOR:
if ( l->select.far_nbr_list != NULL )
{
- free( l->select.far_nbr_list );
+ sfree( l->select.far_nbr_list, "Delete_List::l->select.far_nbr_list" );
}
break;
case TYP_NEAR_NEIGHBOR:
if ( l->select.near_nbr_list != NULL )
{
- free( l->select.near_nbr_list );
+ sfree( l->select.near_nbr_list, "Delete_List::l->select.near_nbr_list" );
}
break;
case TYP_HBOND:
if ( l->select.hbond_list != NULL )
{
- free( l->select.hbond_list );
+ sfree( l->select.hbond_list, "Delete_List::l->select.hbond_list" );
}
break;
default:
- // Report fatal error
+ fprintf( stderr, "[ERROR] unknown list type. Terminating...\n" );
+ exit( INVALID_INPUT );
break;
}
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index 973ba5fa..fef08501 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -21,6 +21,7 @@
#include "lookup.h"
+#include "tool_box.h"
#include "two_body_interactions.h"
@@ -89,49 +90,58 @@ void Natural_Cubic_Spline( const real *h, const real *f,
v = (real*) malloc( n * sizeof(real) );
/* build the linear system */
- a[0] = a[1] = a[n - 1] = 0;
+ a[0] = 0.0;
+ a[1] = 0.0;
+ a[n - 1] = 0.0;
for ( i = 2; i < n - 1; ++i )
{
a[i] = h[i - 1];
}
- b[0] = b[n - 1] = 0;
+ b[0] = 0.0;
+ b[n - 1] = 0.0;
for ( i = 1; i < n - 1; ++i )
{
- b[i] = 2 * (h[i - 1] + h[i]);
+ b[i] = 2.0 * (h[i - 1] + h[i]);
}
- c[0] = c[n - 2] = c[n - 1] = 0;
+ c[0] = 0.0;
+ c[n - 2] = 0.0;
+ c[n - 1] = 0.0;
for ( i = 1; i < n - 2; ++i )
{
c[i] = h[i];
}
- d[0] = d[n - 1] = 0;
+ d[0] = 0.0;
+ d[n - 1] = 0.0;
for ( i = 1; i < n - 1; ++i )
{
- d[i] = 6 * ((f[i + 1] - f[i]) / h[i] - (f[i] - f[i - 1]) / h[i - 1]);
+ d[i] = 6.0 * ((f[i + 1] - f[i])
+ / h[i] - (f[i] - f[i - 1]) / h[i - 1]);
}
/*fprintf( stderr, "i a b c d\n" );
for( i = 0; i < n; ++i )
fprintf( stderr, "%d %f %f %f %f\n", i, a[i], b[i], c[i], d[i] );*/
- v[0] = 0;
- v[n - 1] = 0;
- Tridiagonal_Solve( &(a[1]), &(b[1]), &(c[1]), &(d[1]), &(v[1]), n - 2 );
+ v[0] = 0.0;
+ v[n - 1] = 0.0;
+ Tridiagonal_Solve( a + 1, b + 1, c + 1, d + 1, v + 1, n - 2 );
for ( i = 1; i < n; ++i )
{
- coef[i - 1].d = (v[i] - v[i - 1]) / (6 * h[i - 1]);
- coef[i - 1].c = v[i] / 2;
- coef[i - 1].b = (f[i] - f[i - 1]) / h[i - 1] + h[i - 1] * (2 * v[i] + v[i - 1]) / 6;
+ coef[i - 1].d = (v[i] - v[i - 1]) / (6.0 * h[i - 1]);
+ coef[i - 1].c = v[i] / 2.0;
+ coef[i - 1].b = (f[i] - f[i - 1]) / h[i - 1] + h[i - 1]
+ * (2.0 * v[i] + v[i - 1]) / 6.0;
coef[i - 1].a = f[i];
}
- /*fprintf( stderr, "i v coef\n" );
- for( i = 0; i < n; ++i )
- fprintf( stderr, "%d %f %f %f %f %f\n",
- i, v[i], coef[i].a, coef[i].b, coef[i].c, coef[i].d ); */
+ sfree( a, "Natural_Cubic_Spline::a" );
+ sfree( b, "Natural_Cubic_Spline::b" );
+ sfree( c, "Natural_Cubic_Spline::c" );
+ sfree( d, "Natural_Cubic_Spline::d" );
+ sfree( v, "Natural_Cubic_Spline::v" );
}
@@ -150,49 +160,51 @@ void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
v = (real*) malloc( n * sizeof(real) );
/* build the linear system */
- a[0] = 0;
+ a[0] = 0.0;
for ( i = 1; i < n; ++i )
{
a[i] = h[i - 1];
}
- b[0] = 2 * h[0];
+ b[0] = 2.0 * h[0];
for ( i = 1; i < n; ++i )
{
- b[i] = 2 * (h[i - 1] + h[i]);
+ b[i] = 2.0 * (h[i - 1] + h[i]);
}
- c[n - 1] = 0;
+ c[n - 1] = 0.0;
for ( i = 0; i < n - 1; ++i )
{
c[i] = h[i];
}
- d[0] = 6 * (f[1] - f[0]) / h[0] - 6 * v0;
- d[n - 1] = 6 * vlast - 6 * (f[n - 1] - f[n - 2] / h[n - 2]);
+ d[0] = 6.0 * (f[1] - f[0]) / h[0] - 6.0 * v0;
+ d[n - 1] = 6.0 * vlast - 6.0 * (f[n - 1] - f[n - 2] / h[n - 2]);
for ( i = 1; i < n - 1; ++i )
{
- d[i] = 6 * ((f[i + 1] - f[i]) / h[i] - (f[i] - f[i - 1]) / h[i - 1]);
+ d[i] = 6.0 * ((f[i + 1] - f[i])
+ / h[i] - (f[i] - f[i - 1]) / h[i - 1]);
}
/*fprintf( stderr, "i a b c d\n" );
for( i = 0; i < n; ++i )
fprintf( stderr, "%d %f %f %f %f\n", i, a[i], b[i], c[i], d[i] );*/
- Tridiagonal_Solve( &(a[0]), &(b[0]), &(c[0]), &(d[0]), &(v[0]), n );
- // Tridiagonal_Solve( &(a[1]), &(b[1]), &(c[1]), &(d[1]), &(v[1]), n-2 );
+
+ Tridiagonal_Solve( a, b, c, d, v, n );
for ( i = 1; i < n; ++i )
{
- coef[i - 1].d = (v[i] - v[i - 1]) / (6 * h[i - 1]);
- coef[i - 1].c = v[i] / 2;
- coef[i - 1].b = (f[i] - f[i - 1]) / h[i - 1] + h[i - 1] * (2 * v[i] + v[i - 1]) / 6;
+ coef[i - 1].d = (v[i] - v[i - 1]) / (6.0 * h[i - 1]);
+ coef[i - 1].c = v[i] / 2.0;
+ coef[i - 1].b = (f[i] - f[i - 1]) / h[i - 1] + h[i - 1] * (2.0 * v[i] + v[i - 1]) / 6.0;
coef[i - 1].a = f[i];
}
- /*fprintf( stderr, "i v coef\n" );
- for( i = 0; i < n; ++i )
- fprintf( stderr, "%d %f %f %f %f %f\n",
- i, v[i], coef[i].a, coef[i].b, coef[i].c, coef[i].d ); */
+ sfree( a, "Complete_Cubic_Spline::a" );
+ sfree( b, "Complete_Cubic_Spline::b" );
+ sfree( c, "Complete_Cubic_Spline::c" );
+ sfree( d, "Complete_Cubic_Spline::d" );
+ sfree( v, "Complete_Cubic_Spline::v" );
}
@@ -248,12 +260,12 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control )
num_atom_types = system->reaxprm.num_atom_types;
dr = control->r_cut / control->tabulate;
- h = (real*) malloc( (control->tabulate + 1) * sizeof(real) );
- fh = (real*) malloc( (control->tabulate + 1) * sizeof(real) );
- fvdw = (real*) malloc( (control->tabulate + 1) * sizeof(real) );
- fCEvd = (real*) malloc( (control->tabulate + 1) * sizeof(real) );
- fele = (real*) malloc( (control->tabulate + 1) * sizeof(real) );
- fCEclmb = (real*) malloc( (control->tabulate + 1) * sizeof(real) );
+ h = (real*) calloc( (control->tabulate + 1), sizeof(real) );
+ fh = (real*) calloc( (control->tabulate + 1), sizeof(real) );
+ fvdw = (real*) calloc( (control->tabulate + 1), sizeof(real) );
+ fCEvd = (real*) calloc( (control->tabulate + 1), sizeof(real) );
+ fele = (real*) calloc( (control->tabulate + 1), sizeof(real) );
+ fCEclmb = (real*) calloc( (control->tabulate + 1), sizeof(real) );
/* allocate Long-Range LookUp Table space based on
number of atom types in the ffield file */
@@ -291,17 +303,17 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control )
LR[i][j].dx = dr;
LR[i][j].inv_dx = control->tabulate / control->r_cut;
LR[i][j].y = (LR_data*)
- malloc(LR[i][j].n * sizeof(LR_data));
+ malloc( LR[i][j].n * sizeof(LR_data) );
LR[i][j].H = (cubic_spline_coef*)
- malloc(LR[i][j].n * sizeof(cubic_spline_coef));
+ malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
LR[i][j].vdW = (cubic_spline_coef*)
- malloc(LR[i][j].n * sizeof(cubic_spline_coef));
+ malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
LR[i][j].CEvd = (cubic_spline_coef*)
- malloc(LR[i][j].n * sizeof(cubic_spline_coef));
+ malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
LR[i][j].ele = (cubic_spline_coef*)
- malloc(LR[i][j].n * sizeof(cubic_spline_coef));
+ malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
LR[i][j].CEclmb = (cubic_spline_coef*)
- malloc(LR[i][j].n * sizeof(cubic_spline_coef));
+ malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
for ( r = 1; r <= control->tabulate; ++r )
{
@@ -325,31 +337,39 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control )
}
}
- /*fprintf( stderr, "%-6s %-6s %-6s\n", "r", "h", "fh" );
- for( r = 1; r <= control->tabulate; ++r )
- fprintf( stderr, "%f %f %f\n", r * dr, h[r], fh[r] ); */
Natural_Cubic_Spline( &h[1], &fh[1],
- &(LR[i][j].H[1]), control->tabulate + 1 );
+ &(LR[i][j].H[1]), control->tabulate + 1 );
+
+// fprintf( stderr, "%-6s %-6s %-6s\n", "r", "h", "fh" );
+// for( r = 1; r <= control->tabulate; ++r )
+// fprintf( stderr, "%f %f %f\n", r * dr, h[r], fh[r] );
- /*fprintf( stderr, "%-6s %-6s %-6s\n", "r", "h", "fvdw" );
- for( r = 1; r <= control->tabulate; ++r )
- fprintf( stderr, "%f %f %f\n", r * dr, h[r], fvdw[r] );
- fprintf( stderr, "v0_vdw: %f, vlast_vdw: %f\n", v0_vdw, vlast_vdw );
- */
Complete_Cubic_Spline( &h[1], &fvdw[1], v0_vdw, vlast_vdw,
- &(LR[i][j].vdW[1]), control->tabulate + 1 );
+ &(LR[i][j].vdW[1]), control->tabulate + 1 );
+
+// fprintf( stderr, "%-6s %-6s %-6s\n", "r", "h", "fvdw" );
+// for( r = 1; r <= control->tabulate; ++r )
+// fprintf( stderr, "%f %f %f\n", r * dr, h[r], fvdw[r] );
+// fprintf( stderr, "v0_vdw: %f, vlast_vdw: %f\n", v0_vdw, vlast_vdw );
+
Natural_Cubic_Spline( &h[1], &fCEvd[1],
- &(LR[i][j].CEvd[1]), control->tabulate + 1 );
+ &(LR[i][j].CEvd[1]), control->tabulate + 1 );
+
+// fprintf( stderr, "%-6s %-6s %-6s\n", "r", "h", "fele" );
+// for( r = 1; r <= control->tabulate; ++r )
+// fprintf( stderr, "%f %f %f\n", r * dr, h[r], fele[r] );
+// fprintf( stderr, "v0_ele: %f, vlast_ele: %f\n", v0_ele, vlast_ele );
- /*fprintf( stderr, "%-6s %-6s %-6s\n", "r", "h", "fele" );
- for( r = 1; r <= control->tabulate; ++r )
- fprintf( stderr, "%f %f %f\n", r * dr, h[r], fele[r] );
- fprintf( stderr, "v0_ele: %f, vlast_ele: %f\n", v0_ele, vlast_ele );
- */
Complete_Cubic_Spline( &h[1], &fele[1], v0_ele, vlast_ele,
- &(LR[i][j].ele[1]), control->tabulate + 1 );
+ &(LR[i][j].ele[1]), control->tabulate + 1 );
+
+// fprintf( stderr, "%-6s %-6s %-6s\n", "r", "h", "fele" );
+// for( r = 1; r <= control->tabulate; ++r )
+// fprintf( stderr, "%f %f %f\n", r * dr, h[r], fele[r] );
+// fprintf( stderr, "v0_ele: %f, vlast_ele: %f\n", v0_ele, vlast_ele );
+
Natural_Cubic_Spline( &h[1], &fCEclmb[1],
- &(LR[i][j].CEclmb[1]), control->tabulate + 1 );
+ &(LR[i][j].CEclmb[1]), control->tabulate + 1 );
}
}
}
@@ -403,12 +423,54 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control )
fprintf( stderr, "eele_maxerr: %24.15e\n", eele_maxerr );
*******/
- free(h);
- free(fh);
- free(fvdw);
- free(fCEvd);
- free(fele);
- free(fCEclmb);
+ sfree( h, "Make_LR_Lookup_Table::h" );
+ sfree( fh, "Make_LR_Lookup_Table::fh" );
+ sfree( fvdw, "Make_LR_Lookup_Table::fvdw" );
+ sfree( fCEvd, "Make_LR_Lookup_Table::fCEvd" );
+ sfree( fele, "Make_LR_Lookup_Table::fele" );
+ sfree( fCEclmb, "Make_LR_Lookup_Table::fCEclmb" );
+}
+
+
+void Finalize_LR_Lookup_Table( reax_system *system, control_params *control )
+{
+ int i, j;
+ int num_atom_types;
+ int existing_types[MAX_ATOM_TYPES];
+
+ num_atom_types = system->reaxprm.num_atom_types;
+
+ for ( i = 0; i < MAX_ATOM_TYPES; ++i )
+ {
+ existing_types[i] = 0;
+ }
+ for ( i = 0; i < system->N; ++i )
+ {
+ existing_types[ system->atoms[i].type ] = 1;
+ }
+
+ for ( i = 0; i < num_atom_types; ++i )
+ {
+ if ( existing_types[i] )
+ {
+ for ( j = i; j < num_atom_types; ++j )
+ {
+ if ( existing_types[j] )
+ {
+ sfree( LR[i][j].y, "Finalize_LR_Lookup_Table::LR[i][j].y" );
+ sfree( LR[i][j].H, "Finalize_LR_Lookup_Table::LR[i][j].H" );
+ sfree( LR[i][j].vdW, "Finalize_LR_Lookup_Table::LR[i][j].vdW" );
+ sfree( LR[i][j].CEvd, "Finalize_LR_Lookup_Table::LR[i][j].CEvd" );
+ sfree( LR[i][j].ele, "Finalize_LR_Lookup_Table::LR[i][j].ele" );
+ sfree( LR[i][j].CEclmb, "Finalize_LR_Lookup_Table::LR[i][j].CEclmb" );
+ }
+ }
+ }
+
+ sfree( LR[i], "Finalize_LR_Lookup_Table::LR[i]" );
+ }
+
+ sfree( LR, "Finalize_LR_Lookup_Table::LR" );
}
diff --git a/sPuReMD/src/lookup.h b/sPuReMD/src/lookup.h
index bb1ba468..29720cfb 100644
--- a/sPuReMD/src/lookup.h
+++ b/sPuReMD/src/lookup.h
@@ -37,5 +37,7 @@ real Lookup( real, lookup_table* );
void Make_LR_Lookup_Table( reax_system*, control_params* );
+void Finalize_LR_Lookup_Table( reax_system *, control_params * );
+
#endif
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 91f411e2..a350a654 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -891,7 +891,6 @@ typedef struct
sparse_matrix *L;
sparse_matrix *U;
real *droptol;
- real *w;
real *Hdia_inv;
real *b;
real *b_s;
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index 1963bd14..97280473 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -116,10 +116,10 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
{
nbr_data = &(far_nbrs->select.far_nbr_list[num_far]);
//fprintf (stderr, " %f %f %f \n", nbr_data->dvec[0], nbr_data->dvec[1], nbr_data->dvec[2]);
- if (Are_Far_Neighbors(system->atoms[atom1].x,
- system->atoms[atom2].x,
- &(system->box), control->vlist_cut,
- nbr_data))
+ if ( Are_Far_Neighbors(system->atoms[atom1].x,
+ system->atoms[atom2].x,
+ &(system->box), control->vlist_cut,
+ nbr_data) )
{
nbr_data->nbr = atom2;
++num_far;
@@ -231,10 +231,10 @@ int Estimate_NumNeighbors( reax_system *system, control_params *control,
//if( nbrs[itr+1][0] >= 0 || atom1 > atom2 ) {
if ( atom1 > atom2 )
{
- if (Are_Far_Neighbors(system->atoms[atom1].x,
- system->atoms[atom2].x,
- &(system->box), control->vlist_cut,
- &nbr_data))
+ if ( Are_Far_Neighbors(system->atoms[atom1].x,
+ system->atoms[atom2].x,
+ &(system->box), control->vlist_cut,
+ &nbr_data) )
++num_far;
}
}
diff --git a/sPuReMD/src/reset_utils.c b/sPuReMD/src/reset_utils.c
index 0a364e69..bf5330dd 100644
--- a/sPuReMD/src/reset_utils.c
+++ b/sPuReMD/src/reset_utils.c
@@ -168,9 +168,15 @@ void Reset_Grid( grid *g )
int i, j, k;
for ( i = 0; i < g->ncell[0]; i++ )
+ {
for ( j = 0; j < g->ncell[1]; j++ )
+ {
for ( k = 0; k < g->ncell[2]; k++ )
+ {
g->top[i][j][k] = 0;
+ }
+ }
+ }
}
diff --git a/sPuReMD/src/system_props.h b/sPuReMD/src/system_props.h
index 6a5ed1d9..bf8f36e0 100644
--- a/sPuReMD/src/system_props.h
+++ b/sPuReMD/src/system_props.h
@@ -24,9 +24,6 @@
#include <mytypes.h>
-real Get_Time( );
-
-real Get_Timing_Info( real );
void Temperature_Control( control_params*, simulation_data*, output_controls* );
@@ -40,4 +37,5 @@ void Compute_Pressure( reax_system*, simulation_data*, static_storage* );
void Compute_Pressure_Isotropic( reax_system*, control_params*, simulation_data*, output_controls* );
+
#endif
diff --git a/sPuReMD/src/testmd.c b/sPuReMD/src/testmd.c
index d1159ef1..620b8a01 100644
--- a/sPuReMD/src/testmd.c
+++ b/sPuReMD/src/testmd.c
@@ -20,6 +20,7 @@
----------------------------------------------------------------------*/
#include "mytypes.h"
+
#include "analyze.h"
#include "control.h"
#include "ffield.h"
@@ -31,15 +32,13 @@
#include "reset_utils.h"
#include "restart.h"
#include "system_props.h"
+#include "tool_box.h"
#include "vector.h"
-static void Post_Evolve( reax_system * const system,
- control_params * const control,
- simulation_data * const data,
- static_storage * const workspace,
- list ** const lists,
- output_controls * const out_control )
+static void Post_Evolve( reax_system * const system, control_params * const control,
+ simulation_data * const data, static_storage * const workspace,
+ list ** const lists, output_controls * const out_control )
{
int i;
rvec diff, cross;
@@ -145,13 +144,13 @@ void static Read_System( char * const geo_file,
}
-static void usage(char* argv[])
+static void usage( char* argv[] )
{
fprintf(stderr, "usage: ./%s geometry ffield control\n", argv[0]);
}
-int main(int argc, char* argv[])
+int main( int argc, char* argv[] )
{
reax_system system;
control_params control;
@@ -222,7 +221,7 @@ int main(int argc, char* argv[])
Finalize( &system, &control, &data, &workspace, &lists,
&out_control );
- free( lists );
+ sfree( lists, "main::lists" );
return SUCCESS;
}
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 245a4b76..87ef382a 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -389,12 +389,12 @@ void Deallocate_Tokenizer_Space( char **line, char **backup, char ***tokens )
for ( i = 0; i < MAX_TOKENS; i++ )
{
- free( (*tokens)[i] );
+ sfree( (*tokens)[i], "Deallocate_Tokenizer_Space::tokens[i]" );
}
- free( *line );
- free( *backup );
- free( *tokens );
+ sfree( *line, "Deallocate_Tokenizer_Space::line" );
+ sfree( *backup, "Deallocate_Tokenizer_Space::backup" );
+ sfree( *tokens, "Deallocate_Tokenizer_Space::tokens" );
}
@@ -418,74 +418,161 @@ int Tokenize( char* s, char*** tok )
/***************** taken from lammps ************************/
-/* safe malloc */
-void *smalloc( long n, char *name )
+/* Safe wrapper around libc malloc
+ *
+ * n: num. of bytes to allocated
+ * name: message with details about pointer, used for warnings/errors
+ *
+ * returns: ptr to allocated memory
+ * */
+void *smalloc( size_t n, const char *name )
{
void *ptr;
- if ( n <= 0 )
+ if ( n == 0 )
{
- fprintf( stderr, "WARNING: trying to allocate %ld bytes for array %s. ",
- n, name );
- fprintf( stderr, "returning NULL.\n" );
- return NULL;
+ fprintf( stderr, "[ERROR] failed to allocate %zu bytes for array %s.\n",
+ n, name );
+ exit( INSUFFICIENT_MEMORY );
}
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] requesting memory for %s\n", name );
+ fflush( stderr );
+#endif
+
ptr = malloc( n );
+
if ( ptr == NULL )
{
- fprintf( stderr, "ERROR: failed to allocate %ld bytes for array %s",
- n, name );
+ fprintf( stderr, "[ERROR] failed to allocate %zu bytes for array %s.\n",
+ n, name );
exit( INSUFFICIENT_MEMORY );
}
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] address: %p [SMALLOC]\n", (void *) ptr );
+ fflush( stderr );
+#endif
+
return ptr;
}
-/* safe calloc */
-void *scalloc( int n, int size, char *name )
+/* Safe wrapper around libc realloc
+ *
+ * n: num. of bytes to reallocated
+ * name: message with details about pointer, used for warnings/errors
+ *
+ * returns: ptr to reallocated memory
+ * */
+void* srealloc( void *ptr, size_t n, const char *name )
{
- void *ptr;
+ void *new_ptr;
- if ( n <= 0 )
+ if ( n == 0 )
{
- fprintf( stderr, "WARNING: trying to allocate %d elements for array %s. ",
- n, name );
- fprintf( stderr, "returning NULL.\n" );
- return NULL;
+ fprintf( stderr, "[ERROR] failed to reallocate %zu bytes for array %s.\n",
+ n, name );
+ exit( INSUFFICIENT_MEMORY );
}
- if ( size <= 0 )
+ if ( ptr == NULL )
{
- fprintf( stderr, "WARNING: elements size for array %s is %d. ",
- name, size );
- fprintf( stderr, "returning NULL.\n" );
- return NULL;
+ fprintf( stderr, "[INFO] trying to allocate %zu NEW bytes for array %s.\n",
+ n, name );
}
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] requesting memory for %s\n", name );
+ fflush( stderr );
+#endif
+
+ new_ptr = realloc( ptr, n );
+
+ /* technically, ptr may still be allocated and valid,
+ * but we needed more memory, so abort */
+ if ( new_ptr == NULL )
+ {
+ fprintf( stderr, "[ERROR] failed to reallocate %zu bytes for array %s.\n",
+ n, name );
+ exit( INSUFFICIENT_MEMORY );
+ }
+
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] address: %p [SREALLOC]\n", (void *) new_ptr );
+ fflush( stderr );
+#endif
+
+ return new_ptr;
+}
+
+
+/* Safe wrapper around libc calloc
+ *
+ * n: num. of elements to allocated (each of size bytes)
+ * size: num. of bytes per element
+ * name: message with details about pointer, used for warnings/errors
+ *
+ * returns: ptr to allocated memory, all bits initialized to zeros
+ * */
+void* scalloc( size_t n, size_t size, const char *name )
+{
+ void *ptr;
+
+ if ( n == 0 )
+ {
+ fprintf( stderr, "[ERROR] failed to allocate %zu bytes for array %s.\n",
+ n * size, name );
+ exit( INSUFFICIENT_MEMORY );
+ }
+
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] requesting memory for %s\n", name );
+ fflush( stderr );
+#endif
+
ptr = calloc( n, size );
+
if ( ptr == NULL )
{
- fprintf( stderr, "ERROR: failed to allocate %d bytes for array %s",
- n * size, name );
+ fprintf( stderr, "[ERROR] failed to allocate %zu bytes for array %s.\n",
+ n * size, name );
exit( INSUFFICIENT_MEMORY );
}
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] address: %p [SCALLOC]\n", (void *) ptr );
+ fflush( stderr );
+#endif
+
return ptr;
}
/* safe free */
-void sfree( void *ptr, char *name )
+/* Safe wrapper around libc free
+ *
+ * ptr: pointer to dynamically allocated memory which will be deallocated
+ * name: message with details about pointer, used for warnings/errors
+ * */
+void sfree( void *ptr, const char *name )
{
if ( ptr == NULL )
{
- fprintf( stderr, "WARNING: trying to free the already NULL pointer %s!\n",
- name );
+ fprintf( stderr, "[WARNING] trying to free the already NULL pointer %s!\n",
+ name );
return;
}
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] trying to free pointer %s\n", name );
+ fflush( stderr );
+ fprintf( stderr, "[INFO] address: %p [SFREE]\n", (void *) ptr );
+ fflush( stderr );
+#endif
+
free( ptr );
+
ptr = NULL;
}
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index 38e6ba0c..ba103f9a 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -24,45 +24,68 @@
#include "mytypes.h"
+
/* from box.h */
void Transform( rvec, simulation_box*, char, rvec );
+
void Transform_to_UnitBox( rvec, simulation_box*, char, rvec );
+
void Fit_to_Periodic_Box( simulation_box*, rvec* );
+
//void Box_Touch_Point( simulation_box*, ivec, rvec );
-//int is_Inside_Box( simulation_box*, rvec );
-//int iown_midpoint( simulation_box*, rvec, rvec );
+
+//int is_Inside_Box( simulation_box*, rvec );
+
+//int iown_midpoint( simulation_box*, rvec, rvec );
/* from grid.h */
-/*
-void GridCell_Closest_Point( grid_cell*, grid_cell*, ivec, ivec, rvec );
-void GridCell_to_Box_Points( grid_cell*, ivec, rvec, rvec );
-real DistSqr_between_Special_Points( rvec, rvec );
-real DistSqr_to_Special_Point( rvec, rvec );
-int Relative_Coord_Encoding( ivec );
-*/
+//void GridCell_Closest_Point( grid_cell*, grid_cell*, ivec, ivec, rvec );
+
+//void GridCell_to_Box_Points( grid_cell*, ivec, rvec, rvec );
+
+//real DistSqr_between_Special_Points( rvec, rvec );
+
+//real DistSqr_to_Special_Point( rvec, rvec );
+
+//int Relative_Coord_Encoding( ivec );
/* from geo_tools.h */
void Make_Point( real, real, real, rvec* );
+
int is_Valid_Serial( static_storage*, int );
+
int Check_Input_Range( int, int, int, char* );
+
void Trim_Spaces( char* );
/* from system_props.h */
real Get_Time( );
+
real Get_Timing_Info( real );
+
void Update_Timing_Info( real*, real* );
/* from io_tools.h */
int Get_Atom_Type( reax_interaction*, char* );
+
char *Get_Element( reax_system*, int );
+
char *Get_Atom_Name( reax_system*, int );
+
int Allocate_Tokenizer_Space( char**, char**, char*** );
+
void Deallocate_Tokenizer_Space( char **, char **, char *** );
+
int Tokenize( char*, char*** );
/* from lammps */
-void *smalloc( long, char* );
-void *scalloc( int, int, char* );
-void sfree( void*, char* );
+void *smalloc( size_t, const char * );
+
+void* srealloc( void *, size_t, const char * );
+
+void* scalloc( size_t, size_t, const char * );
+
+void sfree( void *, const char * );
+
#endif
diff --git a/sPuReMD/src/two_body_interactions.c b/sPuReMD/src/two_body_interactions.c
index d5379777..49fa32bb 100644
--- a/sPuReMD/src/two_body_interactions.c
+++ b/sPuReMD/src/two_body_interactions.c
@@ -436,8 +436,8 @@ void LR_vdW_Coulomb( reax_system *system, control_params *control,
/* vdWaals calculations */
- powr_vdW1 = POW(r_ij, p_vdW1);
- powgi_vdW1 = POW( 1.0 / twbp->gamma_w, p_vdW1);
+ powr_vdW1 = POW( r_ij, p_vdW1 );
+ powgi_vdW1 = POW( 1.0 / twbp->gamma_w, p_vdW1 );
fn13 = POW( powr_vdW1 + powgi_vdW1, p_vdW1i );
exp1 = EXP( twbp->alpha * (1.0 - fn13 / twbp->r_vdW) );
@@ -449,7 +449,8 @@ void LR_vdW_Coulomb( reax_system *system, control_params *control,
Tap, r_ij, fn13, twbp->D, Tap * twbp->D * (exp1 - 2.0 * exp2),
powgi_vdW1, p_vdW1, twbp->alpha, twbp->r_vdW, exp1, exp2); */
- dfn13 = POW( powr_vdW1 + powgi_vdW1, p_vdW1i - 1.0) * POW(r_ij, p_vdW1 - 2.0);
+ dfn13 = POW( powr_vdW1 + powgi_vdW1, p_vdW1i - 1.0 )
+ * POW( r_ij, p_vdW1 - 2.0 );
lr->CEvd = dTap * twbp->D * (exp1 - 2 * exp2) -
Tap * twbp->D * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2) * dfn13;
--
GitLab