diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 22bfedc9cf8b1e157ffe35335c9ea210a097b405..89eefe4eb13db3484de4ae717bf13103e84d82de 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -27,38 +27,64 @@
/* allocate space for atoms */
void PreAllocate_Space( reax_system *system, control_params *control,
- static_storage *workspace )
+ static_storage *workspace, int n, int first_run )
{
int i;
- system->atoms = scalloc( system->N, sizeof(reax_atom),
- "PreAllocate_Space::system->atoms" );
- workspace->orig_id = scalloc( system->N, sizeof(int),
- "PreAllocate_Space::workspace->orid_id" );
-
- /* bond restriction info */
- if ( control->restrict_bonds )
+ if ( first_run == TRUE )
{
- workspace->restricted = scalloc( system->N, sizeof(int),
- "PreAllocate_Space::workspace->restricted_atoms" );
+ system->atoms = scalloc( n, sizeof(reax_atom),
+ "PreAllocate_Space::system->atoms" );
+ workspace->orig_id = scalloc( n, sizeof(int),
+ "PreAllocate_Space::workspace->orid_id" );
+
+ /* bond restriction info */
+ if ( control->restrict_bonds )
+ {
+ workspace->restricted = scalloc( n, sizeof(int),
+ "PreAllocate_Space::workspace->restricted_atoms" );
+
+ workspace->restricted_list = scalloc( n, sizeof(int*),
+ "PreAllocate_Space::workspace->restricted_list" );
- workspace->restricted_list = scalloc( system->N, sizeof(int*),
- "PreAllocate_Space::workspace->restricted_list" );
+ for ( i = 0; i < n; ++i )
+ {
+ workspace->restricted_list[i] = scalloc( MAX_RESTRICT, sizeof(int),
+ "PreAllocate_Space::workspace->restricted_list[i]" );
+ }
+ }
+ }
+ else
+ {
+ system->atoms = srealloc( system->atoms, n * sizeof(reax_atom),
+ "PreAllocate_Space::system->atoms" );
+ workspace->orig_id = srealloc( workspace->orig_id, n * sizeof(int),
+ "PreAllocate_Space::workspace->orid_id" );
- for ( i = 0; i < system->N; ++i )
+ /* bond restriction info */
+ if ( control->restrict_bonds )
{
- workspace->restricted_list[i] = scalloc( MAX_RESTRICT, sizeof(int),
- "PreAllocate_Space::workspace->restricted_list[i]" );
+ workspace->restricted = srealloc( workspace->restricted, n * sizeof(int),
+ "PreAllocate_Space::workspace->restricted_atoms" );
+
+ workspace->restricted_list = srealloc( workspace->restricted_list, n * sizeof(int*),
+ "PreAllocate_Space::workspace->restricted_list" );
+
+ for ( i = 0; i < n; ++i )
+ {
+ workspace->restricted_list[i] = srealloc( workspace->restricted_list[i], MAX_RESTRICT * sizeof(int),
+ "PreAllocate_Space::workspace->restricted_list[i]" );
+ }
}
}
}
-static void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs )
+static void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int n_max, int num_intrs )
{
Delete_List( TYP_FAR_NEIGHBOR, far_nbrs );
- Make_List( n, num_intrs, TYP_FAR_NEIGHBOR, far_nbrs );
+ Make_List( n, n_max, num_intrs, TYP_FAR_NEIGHBOR, far_nbrs );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "num_far = %d, far_nbrs = %d -> reallocating!\n",
@@ -72,10 +98,11 @@ static void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs
/* Dynamic allocation of memory for matrix in CSR format
*
* pH (output): pointer to sparse matrix for which to allocate
- * n: dimension of the matrix
+ * n: num. rows of the matrix
+ * n_max: max. num. rows of the matrix
* m: number of nonzeros to allocate space for in matrix
* */
-void Allocate_Matrix( sparse_matrix **pH, int n, int m )
+void Allocate_Matrix( sparse_matrix **pH, int n, int n_max, int m )
{
sparse_matrix *H;
@@ -83,9 +110,10 @@ void Allocate_Matrix( sparse_matrix **pH, int n, int m )
H = *pH;
H->n = n;
+ H->n_max = n_max;
H->m = m;
- H->start = smalloc( sizeof(unsigned int) * (n + 1), "Allocate_Matrix::H->start" );
+ H->start = smalloc( sizeof(unsigned int) * (n_max + 1), "Allocate_Matrix::H->start" );
H->j = smalloc( sizeof(unsigned int) * m, "Allocate_Matrix::H->j" );
H->val = smalloc( sizeof(real) * m, "Allocate_Matrix::H->val" );
}
@@ -104,15 +132,15 @@ void Deallocate_Matrix( sparse_matrix *H )
}
-static void Reallocate_Matrix( sparse_matrix **H, int n, int m )
+static void Reallocate_Matrix( sparse_matrix **H, int n, int n_max, int m )
{
Deallocate_Matrix( *H );
- Allocate_Matrix( H, n, m );
+ Allocate_Matrix( H, n, n_max, m );
}
-void Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
+void Allocate_HBond_List( int n, int num_h, int num_h_max, int *h_index, int *hb_top,
reax_list *hbonds )
{
int i, num_hbonds;
@@ -125,7 +153,7 @@ void Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
}
num_hbonds = hb_top[n - 1];
- Make_List( num_h, num_hbonds, TYP_HBOND, hbonds );
+ Make_List( num_h, num_h_max, num_hbonds, TYP_HBOND, hbonds );
for ( i = 0; i < n; ++i )
{
@@ -149,7 +177,8 @@ void Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
}
-static void Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *hbonds )
+static void Reallocate_HBonds_List( int n, int num_h, int num_h_max,
+ int *h_index, reax_list *hbonds )
{
int i;
int *hb_top;
@@ -165,13 +194,13 @@ static void Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *h
Delete_List( TYP_HBOND, hbonds );
- Allocate_HBond_List( n, num_h, h_index, hb_top, hbonds );
+ Allocate_HBond_List( n, num_h, num_h_max, h_index, hb_top, hbonds );
sfree( hb_top, "Reallocate_HBonds_List::hb_top" );
}
-void Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
+void Allocate_Bond_List( int n, int n_max, int *bond_top, reax_list *bonds )
{
int i, num_bonds;
@@ -183,7 +212,7 @@ void Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
}
num_bonds = bond_top[n - 1];
- Make_List( n, num_bonds, TYP_BOND, bonds );
+ Make_List( n, n_max, num_bonds, TYP_BOND, bonds );
Set_Start_Index( 0, 0, bonds );
Set_End_Index( 0, 0, bonds );
@@ -201,7 +230,8 @@ void Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
}
-static void Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int *est_3body )
+static void Reallocate_Bonds_List( int n, int n_max, reax_list *bonds,
+ int *num_bonds, int *est_3body )
{
int i;
int *bond_top;
@@ -221,7 +251,7 @@ static void Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int
Delete_List( TYP_BOND, bonds );
- Allocate_Bond_List( n, bond_top, bonds );
+ Allocate_Bond_List( n, n_max, bond_top, bonds );
*num_bonds = bond_top[n - 1];
sfree( bond_top, "Reallocate_Bonds_List::bond_top" );
@@ -242,13 +272,13 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
if ( realloc->num_far > 0 && nbr_flag == TRUE )
{
Reallocate_Neighbor_List( lists[FAR_NBRS],
- system->N, realloc->num_far * SAFE_ZONE );
+ system->N, system->N_max, realloc->num_far * SAFE_ZONE );
realloc->num_far = -1;
}
if ( realloc->Htop > 0 )
{
- Reallocate_Matrix( &workspace->H, system->N_cm,
+ Reallocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max,
realloc->Htop * SAFE_ZONE );
realloc->Htop = -1;
@@ -260,15 +290,15 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
if ( control->hbond_cut > 0.0 && realloc->hbonds > 0 )
{
- Reallocate_HBonds_List( system->N, workspace->num_H, workspace->hbond_index,
- lists[HBONDS] );
+ Reallocate_HBonds_List( system->N, workspace->num_H, workspace->num_H_max,
+ workspace->hbond_index, lists[HBONDS] );
realloc->hbonds = -1;
}
num_bonds = est_3body = -1;
if ( realloc->bonds > 0 )
{
- Reallocate_Bonds_List( system->N, lists[BONDS], &num_bonds, &est_3body );
+ Reallocate_Bonds_List( system->N, system->N_max, lists[BONDS], &num_bonds, &est_3body );
realloc->bonds = -1;
realloc->num_3body = MAX( realloc->num_3body, est_3body );
}
@@ -283,14 +313,15 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
}
realloc->num_3body *= SAFE_ZONE;
- Make_List( num_bonds, realloc->num_3body,
+ Make_List( num_bonds, num_bonds, realloc->num_3body,
TYP_THREE_BODY, lists[THREE_BODIES] );
realloc->num_3body = -1;
+
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "reallocating 3 bodies\n" );
- fprintf( stderr, "reallocated - num_bonds: %d\n", num_bonds );
- fprintf( stderr, "reallocated - num_3body: %d\n", realloc->num_3body );
- fprintf( stderr, "reallocated 3body memory: %ldMB\n",
+ fprintf( stderr, "[INFO] reallocating 3 bodies\n" );
+ fprintf( stderr, "[INFO] reallocated - num_bonds: %d\n", num_bonds );
+ fprintf( stderr, "[INFO] reallocated - num_3body: %d\n", realloc->num_3body );
+ fprintf( stderr, "[INFO] reallocated 3body memory: %ldMB\n",
realloc->num_3body * sizeof(three_body_interaction_data) /
(1024 * 1024) );
#endif
@@ -299,16 +330,15 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
if ( realloc->gcell_atoms > -1 )
{
#if defined(DEBUG_FOCUS)
- fprintf(stderr, "reallocating gcell: g->max_atoms: %d\n", g->max_atoms);
+ fprintf( stderr, "[INFO] reallocating gcell: g->max_atoms: %d\n", g->max_atoms );
#endif
- for ( i = 0; i < g->ncell[0]; i++ )
+ for ( i = 0; i < g->ncell_max[0]; i++ )
{
- for ( j = 0; j < g->ncell[1]; j++ )
+ for ( j = 0; j < g->ncell_max[1]; j++ )
{
- for ( k = 0; k < g->ncell[2]; k++ )
+ for ( k = 0; k < g->ncell_max[2]; k++ )
{
- // reallocate g->atoms
sfree( g->atoms[i][j][k], "Reallocate::g->atoms[i][j][k]" );
g->atoms[i][j][k] = scalloc( workspace->realloc.gcell_atoms, sizeof(int),
"Reallocate::g->atoms[i][j][k]" );
diff --git a/sPuReMD/src/allocate.h b/sPuReMD/src/allocate.h
index ae5270d849b78b516a634d705b741f9b93c42d97..e4cd681904a6328ac8bd553c4eb689ad63ad9a4e 100644
--- a/sPuReMD/src/allocate.h
+++ b/sPuReMD/src/allocate.h
@@ -25,17 +25,17 @@
#include "reax_types.h"
-void PreAllocate_Space( reax_system*, control_params*, static_storage* );
+void PreAllocate_Space( reax_system*, control_params*, static_storage*, int, int );
void Reallocate( reax_system*, control_params*, static_storage*, reax_list**, int );
-void Allocate_Matrix( sparse_matrix**, int, int );
+void Allocate_Matrix( sparse_matrix**, int, int, int );
void Deallocate_Matrix( sparse_matrix* );
-void Allocate_HBond_List( int, int, int*, int*, reax_list* );
+void Allocate_HBond_List( int, int, int, int*, int*, reax_list* );
-void Allocate_Bond_List( int, int*, reax_list* );
+void Allocate_Bond_List( int, int, int*, reax_list* );
#endif
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 193de2194aedcefdfb5c2723a60509d4dcb81236..57bc2e0cde0ab6e1d823fbbed4ce0fbf34226b07 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -1186,16 +1186,16 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
if ( workspace->L == NULL )
{
- Allocate_Matrix( &workspace->L, Hptr->n, fillin );
- Allocate_Matrix( &workspace->U, Hptr->n, fillin );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, fillin );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, fillin );
}
else if ( workspace->L->m < fillin )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
- Allocate_Matrix( &workspace->L, Hptr->n, fillin );
- Allocate_Matrix( &workspace->U, Hptr->n, fillin );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, fillin );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, fillin );
}
break;
@@ -1207,16 +1207,16 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
if ( workspace->L == NULL )
{
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
else if ( workspace->L->m < Hptr->m )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
break;
@@ -1228,8 +1228,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
{
/* safest storage estimate is ILU(0) (same as
* lower triangular portion of H), could improve later */
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
else if ( workspace->L->m < Hptr->m )
{
@@ -1238,8 +1238,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
/* safest storage estimate is ILU(0) (same as
* lower triangular portion of H), could improve later */
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
break;
@@ -1317,16 +1317,16 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
if ( workspace->L == NULL )
{
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
else if ( workspace->L->m < Hptr->m )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
break;
@@ -1344,8 +1344,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
{
/* safest storage estimate is ILU(0) (same as
* lower triangular portion of H), could improve later */
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
else if ( workspace->L->m < Hptr->m )
{
@@ -1354,8 +1354,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
/* safest storage estimate is ILU(0) (same as
* lower triangular portion of H), could improve later */
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
break;
@@ -1435,16 +1435,16 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
if ( workspace->L == NULL )
{
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
else if ( workspace->L->m < Hptr->m )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
break;
@@ -1464,8 +1464,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
{
/* safest storage estimate is ILU(0) (same as
* lower triangular portion of H), could improve later */
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
else if ( workspace->L->m < Hptr->m )
{
@@ -1473,8 +1473,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
Deallocate_Matrix( workspace->U );
/* factors have sparsity pattern as H */
- Allocate_Matrix( &workspace->L, Hptr->n, Hptr->m );
- Allocate_Matrix( &workspace->U, Hptr->n, Hptr->m );
+ Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
+ Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
}
break;
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index 1549998a74faebd8e808910096cc2ae1df80b60a..f536dc36c369e4acf8750a66b0b5866b97e015fd 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -26,7 +26,7 @@
#include "tool_box.h"
-void Read_Force_Field( FILE* fp, reax_interaction* reax )
+void Read_Force_Field( FILE *fp, reax_interaction *reax, int first_run )
{
char *s;
char **tmp;
@@ -54,15 +54,23 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
/* reading the number of global parameters */
n = atoi(tmp[0]);
- if (n < 1)
+ if ( n < 1 )
{
fprintf( stderr, "[WARNING] number of globals in ffield file is 0!\n" );
return;
}
+ if ( first_run == TRUE )
+ {
+ reax->gp.l = (real*) smalloc( sizeof(real) * n,
+ "Read_Force_Field::reax->gp-l" );
+ }
+ else if ( reax->gp.n_global < n )
+ {
+ reax->gp.l = (real*) srealloc( reax->gp.l, sizeof(real) * n,
+ "Read_Force_Field::reax->gp-l" );
+ }
reax->gp.n_global = n;
- reax->gp.l = (real*) smalloc( sizeof(real) * n,
- "Read_Force_Field::reax->gp-l" );
/* see reax_types.h for mapping between l[i] and the lambdas used in ff */
for ( i = 0; i < n; i++ )
@@ -78,61 +86,121 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
/* next line is number of atom types and some comments */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- reax->num_atom_types = atoi(tmp[0]);
+ n = atoi(tmp[0]);
/* 3 lines of comments */
- fgets(s, MAX_LINE, fp);
- fgets(s, MAX_LINE, fp);
- fgets(s, MAX_LINE, fp);
+ fgets( s, MAX_LINE, fp );
+ fgets( s, MAX_LINE, fp );
+ fgets( s, MAX_LINE, fp );
- /* Allocating structures in reax_interaction */
- reax->sbp = (single_body_parameters*) scalloc( reax->num_atom_types, sizeof(single_body_parameters),
- "Read_Force_Field::reax->sbp" );
- reax->tbp = (two_body_parameters**) scalloc( reax->num_atom_types, sizeof(two_body_parameters*),
- "Read_Force_Field::reax->tbp" );
- reax->thbp = (three_body_header***) scalloc( reax->num_atom_types, sizeof(three_body_header**),
- "Read_Force_Field::reax->thbp" );
- reax->hbp = (hbond_parameters***) scalloc( reax->num_atom_types, sizeof(hbond_parameters**),
- "Read_Force_Field::reax->hbp" );
- reax->fbp = (four_body_header****) scalloc( reax->num_atom_types, sizeof(four_body_header***),
- "Read_Force_Field::reax->fbp" );
- tor_flag = (char****) scalloc( reax->num_atom_types, sizeof(char***),
+ if ( first_run == TRUE )
+ {
+ /* Allocating structures in reax_interaction */
+ reax->sbp = (single_body_parameters*) scalloc( n, sizeof(single_body_parameters),
+ "Read_Force_Field::reax->sbp" );
+ reax->tbp = (two_body_parameters**) scalloc( n, sizeof(two_body_parameters*),
+ "Read_Force_Field::reax->tbp" );
+ reax->thbp = (three_body_header***) scalloc( n, sizeof(three_body_header**),
+ "Read_Force_Field::reax->thbp" );
+ reax->hbp = (hbond_parameters***) scalloc( n, sizeof(hbond_parameters**),
+ "Read_Force_Field::reax->hbp" );
+ reax->fbp = (four_body_header****) scalloc( n, sizeof(four_body_header***),
+ "Read_Force_Field::reax->fbp" );
+
+ for ( i = 0; i < n; i++ )
+ {
+ reax->tbp[i] = (two_body_parameters*) scalloc( n, sizeof(two_body_parameters),
+ "Read_Force_Field::reax->tbp[i]" );
+ reax->thbp[i] = (three_body_header**) scalloc( n, sizeof(three_body_header*),
+ "Read_Force_Field::reax->thbp[i]" );
+ reax->hbp[i] = (hbond_parameters**) scalloc( n, sizeof(hbond_parameters*),
+ "Read_Force_Field::reax->hbp[i]" );
+ reax->fbp[i] = (four_body_header***) scalloc( n, sizeof(four_body_header**),
+ "Read_Force_Field::reax->fbp[i]" );
+
+ for ( j = 0; j < n; j++ )
+ {
+ reax->thbp[i][j] = (three_body_header*) scalloc( n, sizeof(three_body_header),
+ "Read_Force_Field::reax->thbp[i][j]" );
+ reax->hbp[i][j] = (hbond_parameters*) scalloc( n, sizeof(hbond_parameters),
+ "Read_Force_Field::reax->hbp[i][j]" );
+ reax->fbp[i][j] = (four_body_header**) scalloc( n, sizeof(four_body_header*),
+ "Read_Force_Field::reax->fbp[i][j]" );
+
+ for ( k = 0; k < n; k++ )
+ {
+ reax->fbp[i][j][k] = (four_body_header*) scalloc( n, sizeof(four_body_header),
+ "Read_Force_Field::reax->fbp[i][j][k]" );
+ }
+ }
+ }
+ }
+ else if ( reax->num_atom_types < n )
+ {
+ /* Allocating structures in reax_interaction */
+ reax->sbp = (single_body_parameters*) srealloc( reax->sbp, n * sizeof(single_body_parameters),
+ "Read_Force_Field::reax->sbp" );
+ reax->tbp = (two_body_parameters**) srealloc( reax->tbp, n * sizeof(two_body_parameters*),
+ "Read_Force_Field::reax->tbp" );
+ reax->thbp = (three_body_header***) srealloc( reax->thbp, n * sizeof(three_body_header**),
+ "Read_Force_Field::reax->thbp" );
+ reax->hbp = (hbond_parameters***) srealloc( reax->hbp, n * sizeof(hbond_parameters**),
+ "Read_Force_Field::reax->hbp" );
+ reax->fbp = (four_body_header****) srealloc( reax->fbp, n * sizeof(four_body_header***),
+ "Read_Force_Field::reax->fbp" );
+
+ for ( i = 0; i < n; i++ )
+ {
+ reax->tbp[i] = (two_body_parameters*) srealloc( reax->tbp[i], n * sizeof(two_body_parameters),
+ "Read_Force_Field::reax->tbp[i]" );
+ reax->thbp[i] = (three_body_header**) srealloc( reax->thbp[i], n * sizeof(three_body_header*),
+ "Read_Force_Field::reax->thbp[i]" );
+ reax->hbp[i] = (hbond_parameters**) srealloc( reax->hbp[i], n * sizeof(hbond_parameters*),
+ "Read_Force_Field::reax->hbp[i]" );
+ reax->fbp[i] = (four_body_header***) srealloc( reax->fbp[i], n * sizeof(four_body_header**),
+ "Read_Force_Field::reax->fbp[i]" );
+
+ for ( j = 0; j < n; j++ )
+ {
+ reax->thbp[i][j] = (three_body_header*) srealloc( reax->thbp[i][j], n * sizeof(three_body_header),
+ "Read_Force_Field::reax->thbp[i][j]" );
+ reax->hbp[i][j] = (hbond_parameters*) srealloc( reax->hbp[i][j], n * sizeof(hbond_parameters),
+ "Read_Force_Field::reax->hbp[i][j]" );
+ reax->fbp[i][j] = (four_body_header**) srealloc( reax->fbp[i][j], n * sizeof(four_body_header*),
+ "Read_Force_Field::reax->fbp[i][j]" );
+
+ for ( k = 0; k < n; k++ )
+ {
+ reax->fbp[i][j][k] = (four_body_header*) srealloc( reax->fbp[i][j][k], n * sizeof(four_body_header),
+ "Read_Force_Field::reax->fbp[i][j][k]" );
+ }
+ }
+ }
+ }
+
+ tor_flag = (char****) smalloc( n * sizeof(char***),
"Read_Force_Field::tor_flag" );
- for ( i = 0; i < reax->num_atom_types; i++ )
+ for ( i = 0; i < n; i++ )
{
- reax->tbp[i] = (two_body_parameters*) scalloc( reax->num_atom_types, sizeof(two_body_parameters),
- "Read_Force_Field::reax->tbp[i]" );
- reax->thbp[i] = (three_body_header**) scalloc( reax->num_atom_types, sizeof(three_body_header*),
- "Read_Force_Field::reax->thbp[i]" );
- reax->hbp[i] = (hbond_parameters**) scalloc( reax->num_atom_types, sizeof(hbond_parameters*),
- "Read_Force_Field::reax->hbp[i]" );
- reax->fbp[i] = (four_body_header***) scalloc( reax->num_atom_types, sizeof(four_body_header**),
- "Read_Force_Field::reax->fbp[i]" );
- tor_flag[i] = (char***) scalloc( reax->num_atom_types, sizeof(char**),
+ tor_flag[i] = (char***) smalloc( n * sizeof(char**),
"Read_Force_Field::tor_flag[i]" );
- for ( j = 0; j < reax->num_atom_types; j++ )
+ for ( j = 0; j < n; j++ )
{
- reax->thbp[i][j] = (three_body_header*) scalloc( reax->num_atom_types, sizeof(three_body_header),
- "Read_Force_Field::reax->thbp[i][j]" );
- reax->hbp[i][j] = (hbond_parameters*) scalloc( reax->num_atom_types, sizeof(hbond_parameters),
- "Read_Force_Field::reax->hbp[i][j]" );
- reax->fbp[i][j] = (four_body_header**) scalloc( reax->num_atom_types, sizeof(four_body_header*),
- "Read_Force_Field::reax->fbp[i][j]" );
- tor_flag[i][j] = (char**) scalloc( reax->num_atom_types, sizeof(char*),
+ tor_flag[i][j] = (char**) smalloc( n * sizeof(char*),
"Read_Force_Field::tor_flag[i][j]" );
- for (k = 0; k < reax->num_atom_types; k++)
+ for ( k = 0; k < n; k++ )
{
- reax->fbp[i][j][k] = (four_body_header*) scalloc( reax->num_atom_types, sizeof(four_body_header),
- "Read_Force_Field::reax->fbp[i][j][k]" );
- tor_flag[i][j][k] = (char*) scalloc( reax->num_atom_types, sizeof(char),
+ tor_flag[i][j][k] = (char*) smalloc( n * sizeof(char),
"Read_Force_Field::tor_flag[i][j][k]" );
}
}
}
+ reax->num_atom_types = n;
+
// vdWaals type: 1: Shielded Morse, no inner-wall
// 2: inner wall, no shielding
// 3: inner wall+shielding
@@ -340,7 +408,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
/* a line of comments */
fgets(s, MAX_LINE, fp);
- for (i = 0; i < l; i++)
+ for ( i = 0; i < l; i++ )
{
/* line 1 */
fgets( s, MAX_LINE, fp );
@@ -349,7 +417,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
j = atoi(tmp[0]) - 1;
k = atoi(tmp[1]) - 1;
- if (j < reax->num_atom_types && k < reax->num_atom_types)
+ if ( j < reax->num_atom_types && k < reax->num_atom_types )
{
val = atof(tmp[2]);
@@ -481,35 +549,35 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
}
val = atof(tmp[3]);
- if (val > 0.0)
+ if ( val > 0.0 )
{
reax->tbp[j][k].r_vdW = 2.0 * val;
reax->tbp[k][j].r_vdW = 2.0 * val;
}
val = atof(tmp[4]);
- if (val > 0.0)
+ if ( val > 0.0 )
{
reax->tbp[j][k].alpha = val;
reax->tbp[k][j].alpha = val;
}
val = atof(tmp[5]);
- if (val > 0.0)
+ if ( val > 0.0 )
{
reax->tbp[j][k].r_s = val;
reax->tbp[k][j].r_s = val;
}
val = atof(tmp[6]);
- if (val > 0.0)
+ if ( val > 0.0 )
{
reax->tbp[j][k].r_p = val;
reax->tbp[k][j].r_p = val;
}
val = atof(tmp[7]);
- if (val > 0.0)
+ if ( val > 0.0 )
{
reax->tbp[j][k].r_pp = val;
reax->tbp[k][j].r_pp = val;
@@ -545,9 +613,9 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
k = atoi(tmp[1]) - 1;
m = atoi(tmp[2]) - 1;
- if (j < reax->num_atom_types &&
- k < reax->num_atom_types &&
- m < reax->num_atom_types)
+ if ( j < reax->num_atom_types
+ && k < reax->num_atom_types
+ && m < reax->num_atom_types )
{
cnt = reax->thbp[j][k][m].cnt;
reax->thbp[j][k][m].cnt++;
@@ -623,12 +691,12 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
n = atoi(tmp[3]) - 1;
/* this means the entry is not in compact form */
- if (j >= 0 && n >= 0)
+ if ( j >= 0 && n >= 0 )
{
- if (j < reax->num_atom_types &&
- k < reax->num_atom_types &&
- m < reax->num_atom_types &&
- n < reax->num_atom_types)
+ if ( j < reax->num_atom_types
+ && k < reax->num_atom_types
+ && m < reax->num_atom_types
+ && n < reax->num_atom_types )
{
/* these flags ensure that this entry take precedence
* over the compact form entries */
@@ -679,7 +747,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
// reax->fbp[p][k][m][o].cnt++;
// reax->fbp[o][m][k][p].cnt++;
- if (tor_flag[p][k][m][o] == 0)
+ if ( tor_flag[p][k][m][o] == 0 )
{
reax->fbp[p][k][m][o].prm[0].V1 = atof(tmp[4]);
reax->fbp[p][k][m][o].prm[0].V2 = atof(tmp[5]);
@@ -688,7 +756,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
reax->fbp[p][k][m][o].prm[0].p_cot1 = atof(tmp[8]);
}
- if (tor_flag[o][m][k][p] == 0)
+ if ( tor_flag[o][m][k][p] == 0 )
{
reax->fbp[o][m][k][p].prm[0].V1 = atof(tmp[4]);
reax->fbp[o][m][k][p].prm[0].V2 = atof(tmp[5]);
diff --git a/sPuReMD/src/ffield.h b/sPuReMD/src/ffield.h
index 611a705d7afc3bf823995b75efcf9f0e2778e23b..47b0eb6b4d9fcec9e7fc792f1bd15b97cd022c5d 100644
--- a/sPuReMD/src/ffield.h
+++ b/sPuReMD/src/ffield.h
@@ -25,7 +25,7 @@
#include "reax_types.h"
-void Read_Force_Field( FILE*, reax_interaction* );
+void Read_Force_Field( FILE *, reax_interaction *, int );
#endif
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 4a1650c29c97ea0c9e6086960d5438cf9acead1c..e08760b1195677562152fa0846d9721dd8ca8083 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -1408,7 +1408,7 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
far_nbrs = lists[FAR_NBRS];
- for ( i = 0; i < system->N; ++i )
+ for ( i = 0; i < far_nbrs->n; ++i )
{
atom_i = &system->atoms[i];
type_i = atom_i->type;
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index a60acddca39f52518eb2da8e2af523ec08a019b6..0fef6e0336a0e219f4c7ae17d9d3f56de47142f5 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -229,12 +229,14 @@ static int Read_Box_Info( reax_system *system, FILE *fp, int geo_format )
* fp: file pointer to open geometry file
* geo_format: type of geometry file
* */
-static void Count_Atoms( reax_system *system, FILE *fp, int geo_format )
+static int Count_Atoms( reax_system *system, FILE *fp, int geo_format )
{
char line[MAX_LINE];
+ int n;
+
+ n = 0;
fseek( fp, 0, SEEK_SET );
- system->N = 0;
/* increment number of atoms for each line denoting an atom desc */
switch ( geo_format )
@@ -247,7 +249,7 @@ static void Count_Atoms( reax_system *system, FILE *fp, int geo_format )
if ( strncmp( line, "ATOM", 4 ) == 0
|| strncmp( line, "HETATM", 6 ) == 0 )
{
- system->N++;
+ ++n;
}
}
@@ -259,7 +261,7 @@ static void Count_Atoms( reax_system *system, FILE *fp, int geo_format )
/* second line contains integer count
* of the number of atoms */
- fscanf( fp, " %d", &system->N );
+ fscanf( fp, " %d", &n );
break;
@@ -269,20 +271,22 @@ static void Count_Atoms( reax_system *system, FILE *fp, int geo_format )
break;
}
- assert( system->N > 0 );
+ assert( n > 0 );
#if defined(DEBUG)
fprintf( stderr, "[INFO] Count_Atoms: " );
fprintf( stderr, "N = %d\n", system->N );
#endif
+
+ return n;
}
void Read_Geo( const char * const geo_file, reax_system* system, control_params *control,
- simulation_data *data, static_storage *workspace )
+ simulation_data *data, static_storage *workspace, int first_run )
{
FILE *geo;
- int i, serial, top;
+ int i, n, serial, top;
rvec x;
char element[3], name[9];
reax_atom *atom;
@@ -297,8 +301,12 @@ void Read_Geo( const char * const geo_file, reax_system* system, control_params
}
/* count atoms and allocate storage */
- Count_Atoms( system, geo, CUSTOM );
- PreAllocate_Space( system, control, workspace );
+ n = Count_Atoms( system, geo, CUSTOM );
+ if ( first_run == TRUE || n > system->N )
+ {
+ PreAllocate_Space( system, control, workspace, n, first_run );
+ }
+ system->N = n;
/* parse atom info lines (3rd line and beyond) */
top = 0;
@@ -331,9 +339,9 @@ void Read_Geo( const char * const geo_file, reax_system* system, control_params
void Read_PDB( const char * const pdb_file, reax_system* system, control_params *control,
- simulation_data *data, static_storage *workspace )
+ simulation_data *data, static_storage *workspace, int first_run )
{
- int i, c, c1, pdb_serial, top;
+ int i, n, c, c1, pdb_serial, top;
FILE *pdb;
char **tmp;
char *s, *s1;
@@ -361,8 +369,12 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
exit( INVALID_GEO );
}
- Count_Atoms( system, pdb, PDB );
- PreAllocate_Space( system, control, workspace );
+ n = Count_Atoms( system, pdb, PDB );
+ if ( first_run == TRUE || n > system->N )
+ {
+ PreAllocate_Space( system, control, workspace, n, first_run );
+ }
+ system->N = n;
/* reset file pointer to beginning of file
* and parse the PDB file */
@@ -601,7 +613,7 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
void Read_BGF( const char * const bgf_file, reax_system* system, control_params *control,
- simulation_data *data, static_storage *workspace )
+ simulation_data *data, static_storage *workspace, int first_run )
{
FILE *bgf;
char **tokens;
@@ -614,7 +626,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
char chain_id;
char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
char *endptr;
- int i, atom_cnt, token_cnt, bgf_serial, ratom, crystx_found;
+ int i, n, atom_cnt, token_cnt, bgf_serial, ratom, crystx_found;
endptr = NULL;
ratom = 0;
@@ -626,7 +638,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
&tokens, MAX_TOKENS, MAX_TOKEN_LEN );
/* count number of atoms in the BGF file */
- system->N = 0;
+ n = 0;
line[0] = 0;
while ( fgets( line, MAX_LINE, bgf ) )
@@ -637,7 +649,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
if ( strncmp( tokens[0], "ATOM", 4 ) == 0
|| strncmp( tokens[0], "HETATM", 6 ) == 0 )
{
- ++system->N;
+ ++n;
}
line[0] = 0;
@@ -650,7 +662,11 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
sfclose( bgf, "Read_BGF::bgf" );
- PreAllocate_Space( system, control, workspace );
+ if ( first_run == TRUE || n > system->N )
+ {
+ PreAllocate_Space( system, control, workspace, n, first_run );
+ }
+ system->N = n;
workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
"Read_BGF::workspace->map_serials" );
diff --git a/sPuReMD/src/geo_tools.h b/sPuReMD/src/geo_tools.h
index b72c8bd1357d2bf4663070b1b9a7971dd0785ead..af91cb3c349e50c6c30b2191d8db2a943a8cf856 100644
--- a/sPuReMD/src/geo_tools.h
+++ b/sPuReMD/src/geo_tools.h
@@ -25,17 +25,17 @@
#include "reax_types.h"
-void Read_Geo( const char * const, reax_system*, control_params*,
- simulation_data*, static_storage* );
+void Read_Geo( const char * const, reax_system *, control_params *,
+ simulation_data *, static_storage *, int );
-void Read_PDB( const char * const, reax_system*, control_params*,
- simulation_data*, static_storage* );
+void Read_PDB( const char * const, reax_system *, control_params *,
+ simulation_data *, static_storage *, int );
-void Read_BGF( const char * const, reax_system*, control_params*,
- simulation_data*, static_storage* );
+void Read_BGF( const char * const, reax_system *, control_params *,
+ simulation_data *, static_storage *, int );
-void Write_PDB( reax_system*, reax_list*, simulation_data*,
- control_params*, static_storage*, output_controls* );
+void Write_PDB( reax_system *, reax_list *, simulation_data *,
+ control_params *, static_storage *, output_controls * );
#endif
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index ede82a292d8fba1e0c19386a35da68bd7a508006..7e6d9b5fc4fb4c2b3f5bff11da92c50bbd19a007 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -38,9 +38,11 @@ static int Estimate_GCell_Population( reax_system* system )
for ( l = 0; l < system->N; l++ )
{
- i = (int)(system->atoms[l].x[0] * g->inv_len[0]);
- j = (int)(system->atoms[l].x[1] * g->inv_len[1]);
- k = (int)(system->atoms[l].x[2] * g->inv_len[2]);
+ /* NOTE: this assumes the atom positions
+ * are within the simulation box boundaries */
+ i = (int) (system->atoms[l].x[0] * g->inv_len[0]);
+ 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",
@@ -67,109 +69,197 @@ static int Estimate_GCell_Population( reax_system* system )
}
-static void Allocate_Space_for_Grid( reax_system *system )
+static void Allocate_Space_for_Grid( reax_system *system, int alloc )
{
- int i, j, k, l;
+ int i, j, k, l, max_atoms;
grid *g;
g = &system->g;
g->max_nbrs = (2 * g->spread[0] + 1)
* (2 * g->spread[1] + 1) * (2 * g->spread[2] + 1) + 3;
- /* allocate space for the new grid */
- g->atoms = (int****) scalloc( g->ncell[0], sizeof( int*** ),
- "Allocate_Space_for_Grid::g->atoms" );
- g->top = (int***) scalloc( g->ncell[0], sizeof( int** ),
- "Allocate_Space_for_Grid::g->top" );
- g->mark = (int***) scalloc( g->ncell[0], sizeof( int** ),
- "Allocate_Space_for_Grid::g->mark" );
- g->start = (int***) scalloc( g->ncell[0], sizeof( int** ),
- "Allocate_Space_for_Grid::g->start" );
- g->end = (int***) scalloc( g->ncell[0], sizeof( int** ),
- "Allocate_Space_for_Grid::g->end" );
- g->nbrs = (ivec****) scalloc( g->ncell[0], sizeof( ivec*** ),
- "Allocate_Space_for_Grid::g->nbrs" );
- g->nbrs_cp = (rvec****) scalloc( g->ncell[0], sizeof( rvec*** ),
- "Allocate_Space_for_Grid::g->nbrs_cp" );
-
- for ( i = 0; i < g->ncell[0]; i++ )
+ if ( alloc == TRUE )
{
- g->atoms[i] = (int***) scalloc( g->ncell[1], sizeof( int** ),
- "Allocate_Space_for_Grid::g->atoms[i]" );
- g->top [i] = (int**) scalloc( g->ncell[1], sizeof( int* ),
- "Allocate_Space_for_Grid::g->top[i]" );
- g->mark[i] = (int**) scalloc( g->ncell[1], sizeof( int* ),
- "Allocate_Space_for_Grid::g->mark[i]" );
- g->start[i] = (int**) scalloc( g->ncell[1], sizeof( int* ),
- "Allocate_Space_for_Grid::g->start[i]" );
- g->end[i] = (int**) scalloc( g->ncell[1], sizeof( int* ),
- "Allocate_Space_for_Grid::g->end[i]" );
- g->nbrs[i] = (ivec***) scalloc( g->ncell[1], sizeof( ivec** ),
- "Allocate_Space_for_Grid::g->nbrs[i]" );
- g->nbrs_cp[i] = (rvec***) scalloc( g->ncell[1], sizeof( rvec** ),
- "Allocate_Space_for_Grid::g->nbrs_cp[i]" );
+ /* allocate space for the new grid */
+ g->atoms = (int****) scalloc( g->ncell_max[0], sizeof( int*** ),
+ "Allocate_Space_for_Grid::g->atoms" );
- for ( j = 0; j < g->ncell[1]; j++ )
+ for ( i = 0; i < g->ncell_max[0]; i++ )
{
- g->atoms[i][j] = (int**) scalloc( g->ncell[2], sizeof( int* ),
- "Allocate_Space_for_Grid::g->atoms[i][j]" );
- g->top[i][j] = (int*) scalloc( g->ncell[2], sizeof( int ),
- "Allocate_Space_for_Grid::g->top[i][j]" );
- g->mark[i][j] = (int*) scalloc( g->ncell[2], sizeof( int ),
- "Allocate_Space_for_Grid::g->mark[i][j]" );
- g->start[i][j] = (int*) scalloc( g->ncell[2], sizeof( int ),
- "Allocate_Space_for_Grid::g->start[i][j]" );
- g->end[i][j] = (int*) scalloc( g->ncell[2], sizeof( int ),
- "Allocate_Space_for_Grid::g->end[i][j]" );
- g->nbrs[i][j] = (ivec**) scalloc( g->ncell[2], sizeof( ivec* ),
- "Allocate_Space_for_Grid::g->nbrs[i][j]" );
- g->nbrs_cp[i][j] = (rvec**) scalloc( g->ncell[2], sizeof( rvec* ),
- "Allocate_Space_for_Grid::g->nbrs_cp[i][j]" );
+ g->atoms[i] = (int***) scalloc( g->ncell_max[1], sizeof( int** ),
+ "Allocate_Space_for_Grid::g->atoms[i]" );
- for ( k = 0; k < g->ncell[2]; k++ )
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ g->atoms[i][j] = (int**) scalloc( g->ncell_max[2], sizeof( int* ),
+ "Allocate_Space_for_Grid::g->atoms[i][j]" );
+
+ }
+
+ g->top = (int***) scalloc( g->ncell_max[0], sizeof( int** ),
+ "Allocate_Space_for_Grid::g->top" );
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ {
+ g->top[i] = (int**) scalloc( g->ncell_max[1], sizeof( int* ),
+ "Allocate_Space_for_Grid::g->top[i]" );
+
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ g->top[i][j] = (int*) scalloc( g->ncell_max[2], sizeof( int ),
+ "Allocate_Space_for_Grid::g->top[i][j]" );
+ }
+
+ g->mark = (int***) scalloc( g->ncell_max[0], sizeof( int** ),
+ "Allocate_Space_for_Grid::g->mark" );
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ {
+ g->mark[i] = (int**) scalloc( g->ncell_max[1], sizeof( int* ),
+ "Allocate_Space_for_Grid::g->mark[i]" );
+
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ g->mark[i][j] = (int*) scalloc( g->ncell_max[2], sizeof( int ),
+ "Allocate_Space_for_Grid::g->mark[i][j]" );
+ }
+
+ g->start = (int***) scalloc( g->ncell_max[0], sizeof( int** ),
+ "Allocate_Space_for_Grid::g->start" );
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ {
+ g->start[i] = (int**) scalloc( g->ncell_max[1], sizeof( int* ),
+ "Allocate_Space_for_Grid::g->start[i]" );
+
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ g->start[i][j] = (int*) scalloc( g->ncell_max[2], sizeof( int ),
+ "Allocate_Space_for_Grid::g->start[i][j]" );
+ }
+
+ g->end = (int***) scalloc( g->ncell_max[0], sizeof( int** ),
+ "Allocate_Space_for_Grid::g->end" );
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ {
+ g->end[i] = (int**) scalloc( g->ncell_max[1], sizeof( int* ),
+ "Allocate_Space_for_Grid::g->end[i]" );
+
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ g->end[i][j] = (int*) scalloc( g->ncell_max[2], sizeof( int ),
+ "Allocate_Space_for_Grid::g->end[i][j]" );
+ }
+
+ g->nbrs = (ivec****) scalloc( g->ncell_max[0], sizeof( ivec*** ),
+ "Allocate_Space_for_Grid::g->nbrs" );
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ {
+ g->nbrs[i] = (ivec***) scalloc( g->ncell_max[1], sizeof( ivec** ),
+ "Allocate_Space_for_Grid::g->nbrs[i]" );
+
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ {
+ g->nbrs[i][j] = (ivec**) scalloc( g->ncell_max[2], sizeof( ivec* ),
+ "Allocate_Space_for_Grid::g->nbrs[i][j]" );
+
+ for ( k = 0; k < g->ncell_max[2]; k++ )
+ g->nbrs[i][j][k] = (ivec*) smalloc( g->max_nbrs * sizeof( ivec ),
+ "Allocate_Space_for_Grid::g->nbrs[i][j]][k]" );
+ }
+ }
+
+ g->nbrs_cp = (rvec****) scalloc( g->ncell_max[0], sizeof( rvec*** ),
+ "Allocate_Space_for_Grid::g->nbrs_cp" );
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ {
+ g->nbrs_cp[i] = (rvec***) scalloc( g->ncell_max[1], sizeof( rvec** ),
+ "Allocate_Space_for_Grid::g->nbrs_cp[i]" );
+
+ for ( j = 0; j < g->ncell_max[1]; j++ )
{
+ g->nbrs_cp[i][j] = (rvec**) scalloc( g->ncell_max[2], sizeof( rvec* ),
+ "Allocate_Space_for_Grid::g->nbrs_cp[i][j]" );
+
+ for ( k = 0; k < g->ncell_max[2]; k++ )
+ g->nbrs_cp[i][j][k] = (rvec*) smalloc( g->max_nbrs * sizeof( rvec ),
+ "Allocate_Space_for_Grid::g->nbrs_cp[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;
+
+ 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->mark[i][j][k] = 0;
+
+ 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->start[i][j][k] = 0;
+
+ 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->end[i][j][k] = 0;
- g->nbrs[i][j][k] = (ivec*) smalloc( g->max_nbrs * sizeof( ivec ),
- "Allocate_Space_for_Grid::g->nbrs[i][j]][k]" );
- g->nbrs_cp[i][j][k] = (rvec*) smalloc( g->max_nbrs * sizeof( rvec ),
- "Allocate_Space_for_Grid::g->nbrs_cp[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++ )
for ( l = 0; l < g->max_nbrs; ++l )
{
g->nbrs[i][j][k][l][0] = -1;
g->nbrs[i][j][k][l][1] = -1;
g->nbrs[i][j][k][l][2] = -1;
+ }
+ for ( i = 0; i < g->ncell[0]; i++ )
+ for ( j = 0; j < g->ncell[1]; j++ )
+ for ( k = 0; k < g->ncell[2]; k++ )
+ for ( l = 0; l < g->max_nbrs; ++l )
+ {
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 );
+ max_atoms = Estimate_GCell_Population( system );
- for ( i = 0; i < g->ncell[0]; i++ )
+ if ( alloc == TRUE )
+ {
+ g->max_atoms = MAX( max_atoms, g->max_atoms );
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ for ( k = 0; k < g->ncell_max[2]; k++ )
+ g->atoms[i][j][k] = (int*) smalloc( g->max_atoms * sizeof( int ),
+ "Allocate_Space_for_Grid::g->atoms[i][j]][k]" );
+ }
+ /* case: grid large enough but max. atoms per grid cells insufficient */
+ else if ( g->max_atoms > 0 && g->max_atoms < max_atoms )
{
+ g->max_atoms = max_atoms;
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ for ( k = 0; k < g->ncell_max[2]; k++ )
+ sfree( g->atoms[i][j][k], "Allocate_Space_for_Grid::g->atoms[i][j]][k]" );
+
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ for ( k = 0; k < g->ncell_max[2]; k++ )
+ g->atoms[i][j][k] = (int*) smalloc( g->max_atoms * sizeof( int ),
+ "Allocate_Space_for_Grid::g->atoms[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->atoms[i][j][k] = (int*) smalloc( g->max_atoms * sizeof( int ),
- "Allocate_Space_for_Grid::g->atoms[i][j]][k]" );
-
for ( l = 0; l < g->max_atoms; ++l )
- {
g->atoms[i][j][k][l] = -1;
- }
- }
- }
- }
-
}
@@ -178,11 +268,11 @@ static void Deallocate_Grid_Space( grid *g )
int i, j, k;
/* deallocate the old grid */
- for ( i = 0; i < g->ncell[0]; i++ )
+ for ( i = 0; i < g->ncell_max[0]; i++ )
{
- for ( j = 0; j < g->ncell[1]; j++ )
+ for ( j = 0; j < g->ncell_max[1]; j++ )
{
- for ( k = 0; k < g->ncell[2]; k++ )
+ for ( k = 0; k < g->ncell_max[2]; 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]" );
@@ -372,36 +462,53 @@ static void Find_Neighbor_Grid_Cells( grid *g )
}
-void Setup_Grid( reax_system* system )
+void Setup_Grid( reax_system* system, int first_run )
{
- int d;
- ivec ncell;
+ int d, alloc;
grid *g;
simulation_box *my_box;
+ alloc = FALSE;
g = &system->g;
my_box = &system->box;
/* determine number of grid cells in each direction */
- ivec_rScale( ncell, 1.0 / g->cell_size, my_box->box_norms );
+ ivec_rScale( g->ncell, 1.0 / g->cell_size, my_box->box_norms );
for ( d = 0; d < 3; ++d )
{
- if ( ncell[d] <= 0 )
+ if ( g->ncell[d] <= 0 )
{
- ncell[d] = 1;
+ g->ncell[d] = 1;
}
}
- /* find the number of grid cells */
- g->total = ncell[0] * ncell[1] * ncell[2];
- ivec_Copy( g->ncell, ncell );
+ if ( first_run == TRUE )
+ {
+ g->ncell_max[0] = (int) CEIL( SAFE_ZONE * g->ncell[0] );
+ g->ncell_max[1] = (int) CEIL( SAFE_ZONE * g->ncell[1] );
+ g->ncell_max[2] = (int) CEIL( SAFE_ZONE * g->ncell[2] );
+ g->max_atoms = 0;
+
+ alloc = TRUE;
+ }
+ else if ( first_run == FALSE && (g->ncell[0] > g->ncell_max[0]
+ || g->ncell[1] > g->ncell_max[1] || g->ncell[2] > g->ncell_max[2]) )
+ {
+ Deallocate_Grid_Space( g );
+
+ g->ncell_max[0] = (int) CEIL( SAFE_ZONE * MAX( g->ncell[0], g->ncell_max[0] ) );
+ g->ncell_max[1] = (int) CEIL( SAFE_ZONE * MAX( g->ncell[1], g->ncell_max[1] ) );
+ g->ncell_max[2] = (int) CEIL( SAFE_ZONE * MAX( g->ncell[2], g->ncell_max[2] ) );
+
+ alloc = TRUE;
+ }
/* compute cell lengths */
rvec_iDivide( g->len, my_box->box_norms, g->ncell );
rvec_Invert( g->inv_len, g->len );
- Allocate_Space_for_Grid( system );
+ Allocate_Space_for_Grid( system, alloc );
Find_Neighbor_Grid_Cells( g );
#if defined(DEBUG_FOCUS)
@@ -473,14 +580,13 @@ void Update_Grid( reax_system* system )
Deallocate_Grid_Space( g );
/* update number of grid cells */
- g->total = ncell[0] * ncell[1] * ncell[2];
ivec_Copy( g->ncell, ncell );
/* update cell lengths */
rvec_iDivide( g->len, my_box->box_norms, g->ncell );
rvec_Invert( g->inv_len, g->len );
- Allocate_Space_for_Grid( system );
+ Allocate_Space_for_Grid( system, TRUE );
Find_Neighbor_Grid_Cells( g );
#if defined(DEBUG_FOCUS)
@@ -539,9 +645,9 @@ void Bin_Atoms( reax_system* system, static_storage *workspace )
}
/* reallocation check */
- if ( max_atoms >= (int)(g->max_atoms * SAFE_ZONE) )
+ if ( max_atoms >= (int) CEIL( g->max_atoms * SAFE_ZONE ) )
{
- workspace->realloc.gcell_atoms = MAX( (int)(max_atoms * SAFE_ZONE),
+ workspace->realloc.gcell_atoms = MAX( (int) CEIL( max_atoms * SAFE_ZONE ),
MIN_GCELL_POPL );
}
}
diff --git a/sPuReMD/src/grid.h b/sPuReMD/src/grid.h
index abd0f401bf44462682bd5003be40f63a4abfc7c7..cdc030a318af63e25d333698c9fe4f781d9b0084 100644
--- a/sPuReMD/src/grid.h
+++ b/sPuReMD/src/grid.h
@@ -25,13 +25,13 @@
#include "reax_types.h"
-void Setup_Grid( reax_system* );
+void Setup_Grid( reax_system *, int );
-void Update_Grid( reax_system* );
+void Update_Grid( reax_system * );
-void Bin_Atoms( reax_system*, static_storage* );
+void Bin_Atoms( reax_system *, static_storage * );
-void Finalize_Grid( reax_system* );
+void Finalize_Grid( reax_system * );
void Cluster_Atoms( reax_system *, static_storage *, control_params * );
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 48b8bdef6a8b7dbb35df1bc3e8e4ef942b7b3e74..9809a3bfd57f19b6e48217b2d181bcaa812c290b 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -99,7 +99,7 @@ static void Generate_Initial_Velocities( reax_system *system,
static void Init_System( reax_system *system, control_params *control,
- simulation_data *data )
+ simulation_data *data, int first_run )
{
int i;
rvec dx;
@@ -110,7 +110,7 @@ static void Init_System( reax_system *system, control_params *control,
}
Compute_Total_Mass( system, data );
- Compute_Center_of_Mass( system, data, stderr );
+ Compute_Center_of_Mass( system, data );
/* just fit the atoms to the periodic box */
if ( control->reposition_atoms == 0 )
@@ -154,17 +154,17 @@ static void Init_System( reax_system *system, control_params *control,
Generate_Initial_Velocities( system, control, control->T_init );
}
- Setup_Grid( system );
+ Setup_Grid( system, first_run );
}
static void Init_Simulation_Data( reax_system *system, control_params *control,
simulation_data *data, output_controls *out_control,
- evolve_function *Evolve )
+ evolve_function *Evolve, int alloc )
{
#if defined(_OPENMP)
- if ( control->ensemble == sNPT || control->ensemble == iNPT
- || control->ensemble == aNPT || control->compute_pressure == TRUE )
+ if ( alloc == TRUE && (control->ensemble == sNPT || control->ensemble == iNPT
+ || control->ensemble == aNPT || control->compute_pressure == TRUE) )
{
data->press_local = smalloc( sizeof( rtensor ) * control->num_threads,
"Init_Simulation_Data::data->press_local" );
@@ -323,60 +323,66 @@ static void Init_Taper( control_params *control, static_storage *workspace )
static void Init_Workspace( reax_system *system, control_params *control,
- static_storage *workspace )
+ static_storage *workspace, int alloc )
{
int i;
- /* hydrogen bond list */
- workspace->hbond_index = smalloc( system->N * sizeof( int ),
- "Init_Workspace::workspace->hbond_index" );
-
- /* bond order related storage */
- workspace->total_bond_order = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->bond_order" );
- workspace->Deltap = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->Deltap" );
- workspace->Deltap_boc = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->Deltap_boc" );
- workspace->dDeltap_self = smalloc( system->N * sizeof( rvec ),
- "Init_Workspace::workspace->dDeltap_self" );
-
- workspace->Delta = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->Delta" );
- workspace->Delta_lp = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->Delta_lp" );
- workspace->Delta_lp_temp = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->Delta_lp_temp" );
- workspace->dDelta_lp = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->dDelta_lp" );
- workspace->dDelta_lp_temp = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->dDelta_lp_temp" );
- workspace->Delta_e = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->Delta_e" );
- workspace->Delta_boc = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->Delta_boc" );
- workspace->nlp = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->nlp" );
- workspace->nlp_temp = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->nlp_temp" );
- workspace->Clp = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->Clp" );
- workspace->CdDelta = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->CdDelta" );
- workspace->vlpex = smalloc( system->N * sizeof( real ),
- "Init_Workspace::workspace->vlpex" );
+ if ( alloc == TRUE )
+ {
+ /* hydrogen bond list */
+ workspace->hbond_index = smalloc( system->N_max * sizeof( int ),
+ "Init_Workspace::workspace->hbond_index" );
+
+ /* bond order related storage */
+ workspace->total_bond_order = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->bond_order" );
+ workspace->Deltap = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->Deltap" );
+ workspace->Deltap_boc = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->Deltap_boc" );
+ workspace->dDeltap_self = smalloc( system->N_max * sizeof( rvec ),
+ "Init_Workspace::workspace->dDeltap_self" );
+
+ workspace->Delta = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->Delta" );
+ workspace->Delta_lp = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->Delta_lp" );
+ workspace->Delta_lp_temp = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->Delta_lp_temp" );
+ workspace->dDelta_lp = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->dDelta_lp" );
+ workspace->dDelta_lp_temp = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->dDelta_lp_temp" );
+ workspace->Delta_e = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->Delta_e" );
+ workspace->Delta_boc = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->Delta_boc" );
+ workspace->nlp = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->nlp" );
+ workspace->nlp_temp = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->nlp_temp" );
+ workspace->Clp = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->Clp" );
+ workspace->CdDelta = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->CdDelta" );
+ workspace->vlpex = smalloc( system->N_max * sizeof( real ),
+ "Init_Workspace::workspace->vlpex" );
+ }
/* charge method storage */
switch ( control->charge_method )
{
case QEQ_CM:
system->N_cm = system->N;
+ system->N_cm_max = system->N_max;
break;
case EE_CM:
system->N_cm = system->N + 1;
+ system->N_cm_max = system->N_max + 1;
break;
case ACKS2_CM:
system->N_cm = 2 * system->N + 2;
+ system->N_cm_max = 2 * system->N_max + 2;
break;
default:
fprintf( stderr, "[ERROR] Unknown charge method type. Terminating...\n" );
@@ -384,44 +390,47 @@ static void Init_Workspace( reax_system *system, control_params *control,
break;
}
- /* sparse matrices */
- workspace->H = NULL;
- workspace->H_full = NULL;
- workspace->H_sp = NULL;
- workspace->H_p = NULL;
- workspace->H_spar_patt = NULL;
- workspace->H_spar_patt_full = NULL;
- workspace->H_app_inv = NULL;
- workspace->L = NULL;
- workspace->U = NULL;
- workspace->Hdia_inv = NULL;
- if ( control->cm_solver_pre_comp_type == ICHOLT_PC
- || (control->cm_solver_pre_comp_type == ILUT_PC && control->cm_solver_pre_comp_droptol > 0.0 )
- || control->cm_solver_pre_comp_type == ILUTP_PC
- || control->cm_solver_pre_comp_type == FG_ILUT_PC )
- {
- workspace->droptol = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->droptol" );
- }
-
- workspace->b_s = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->b_s" );
- workspace->b_t = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->b_t" );
- workspace->b_prc = scalloc( system->N_cm * 2, sizeof( real ),
- "Init_Workspace::workspace->b_prc" );
- workspace->b_prm = scalloc( system->N_cm * 2, sizeof( real ),
- "Init_Workspace::workspace->b_prm" );
- workspace->s = scalloc( 5, sizeof( real* ),
- "Init_Workspace::workspace->s" );
- workspace->t = scalloc( 5, sizeof( real* ),
- "Init_Workspace::workspace->t" );
- for ( i = 0; i < 5; ++i )
- {
- workspace->s[i] = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->s[i]" );
- workspace->t[i] = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->t[i]" );
+ if ( alloc == TRUE )
+ {
+ /* sparse matrices */
+ workspace->H = NULL;
+ workspace->H_full = NULL;
+ workspace->H_sp = NULL;
+ workspace->H_p = NULL;
+ workspace->H_spar_patt = NULL;
+ workspace->H_spar_patt_full = NULL;
+ workspace->H_app_inv = NULL;
+ workspace->L = NULL;
+ workspace->U = NULL;
+ workspace->Hdia_inv = NULL;
+ if ( control->cm_solver_pre_comp_type == ICHOLT_PC
+ || (control->cm_solver_pre_comp_type == ILUT_PC && control->cm_solver_pre_comp_droptol > 0.0 )
+ || control->cm_solver_pre_comp_type == ILUTP_PC
+ || control->cm_solver_pre_comp_type == FG_ILUT_PC )
+ {
+ workspace->droptol = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->droptol" );
+ }
+
+ workspace->b_s = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->b_s" );
+ workspace->b_t = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->b_t" );
+ workspace->b_prc = scalloc( system->N_cm_max * 2, sizeof( real ),
+ "Init_Workspace::workspace->b_prc" );
+ workspace->b_prm = scalloc( system->N_cm_max * 2, sizeof( real ),
+ "Init_Workspace::workspace->b_prm" );
+ workspace->s = scalloc( 5, sizeof( real* ),
+ "Init_Workspace::workspace->s" );
+ workspace->t = scalloc( 5, sizeof( real* ),
+ "Init_Workspace::workspace->t" );
+ for ( i = 0; i < 5; ++i )
+ {
+ workspace->s[i] = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->s[i]" );
+ workspace->t[i] = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->t[i]" );
+ }
}
switch ( control->charge_method )
@@ -469,274 +478,283 @@ static void Init_Workspace( reax_system *system, control_params *control,
break;
}
- switch ( control->cm_solver_type )
+ if ( alloc == TRUE )
{
- case GMRES_S:
- case GMRES_H_S:
- workspace->y = scalloc( control->cm_solver_restart + 1, sizeof( real ),
- "Init_Workspace::workspace->y" );
- workspace->z = scalloc( control->cm_solver_restart + 1, sizeof( real ),
- "Init_Workspace::workspace->z" );
- workspace->g = scalloc( control->cm_solver_restart + 1, sizeof( real ),
- "Init_Workspace::workspace->g" );
- workspace->h = scalloc( control->cm_solver_restart + 1, sizeof( real*),
- "Init_Workspace::workspace->h" );
- workspace->hs = scalloc( control->cm_solver_restart + 1, sizeof( real ),
- "Init_Workspace::workspace->hs" );
- workspace->hc = scalloc( control->cm_solver_restart + 1, sizeof( real ),
- "Init_Workspace::workspace->hc" );
- workspace->rn = scalloc( control->cm_solver_restart + 1, sizeof( real*),
- "Init_Workspace::workspace->rn" );
- workspace->v = scalloc( control->cm_solver_restart + 1, sizeof( real*),
- "Init_Workspace::workspace->v" );
-
- for ( i = 0; i < control->cm_solver_restart + 1; ++i )
- {
- workspace->h[i] = scalloc( control->cm_solver_restart + 1, sizeof( real ),
- "Init_Workspace::workspace->h[i]" );
- workspace->rn[i] = scalloc( system->N_cm * 2, sizeof( real ),
- "Init_Workspace::workspace->rn[i]" );
- workspace->v[i] = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->v[i]" );
- }
-
- workspace->r = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->r" );
- workspace->d = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->d" );
- workspace->q = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->q" );
- workspace->p = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->p" );
- break;
-
- case CG_S:
- workspace->r = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->r" );
- workspace->d = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->d" );
- workspace->q = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->q" );
- workspace->p = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->p" );
- break;
-
- case SDM_S:
- workspace->r = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->r" );
- workspace->d = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->d" );
- workspace->q = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->q" );
- break;
-
- case BiCGStab_S:
- workspace->r = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->r" );
- workspace->r_hat = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->r_hat" );
- workspace->d = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->d" );
- workspace->q = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->q" );
- workspace->q_hat = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->q_hat" );
- workspace->p = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->p" );
- workspace->y = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->y" );
- workspace->z = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->z" );
- workspace->g = scalloc( system->N_cm, sizeof( real ),
- "Init_Workspace::workspace->g" );
- break;
-
- default:
- fprintf( stderr, "[ERROR] Unknown charge method linear solver type. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
- }
+ switch ( control->cm_solver_type )
+ {
+ case GMRES_S:
+ case GMRES_H_S:
+ workspace->y = scalloc( control->cm_solver_restart + 1, sizeof( real ),
+ "Init_Workspace::workspace->y" );
+ workspace->z = scalloc( control->cm_solver_restart + 1, sizeof( real ),
+ "Init_Workspace::workspace->z" );
+ workspace->g = scalloc( control->cm_solver_restart + 1, sizeof( real ),
+ "Init_Workspace::workspace->g" );
+ workspace->h = scalloc( control->cm_solver_restart + 1, sizeof( real*),
+ "Init_Workspace::workspace->h" );
+ workspace->hs = scalloc( control->cm_solver_restart + 1, sizeof( real ),
+ "Init_Workspace::workspace->hs" );
+ workspace->hc = scalloc( control->cm_solver_restart + 1, sizeof( real ),
+ "Init_Workspace::workspace->hc" );
+ workspace->rn = scalloc( control->cm_solver_restart + 1, sizeof( real*),
+ "Init_Workspace::workspace->rn" );
+ workspace->v = scalloc( control->cm_solver_restart + 1, sizeof( real*),
+ "Init_Workspace::workspace->v" );
+
+ for ( i = 0; i < control->cm_solver_restart + 1; ++i )
+ {
+ workspace->h[i] = scalloc( control->cm_solver_restart + 1, sizeof( real ),
+ "Init_Workspace::workspace->h[i]" );
+ workspace->rn[i] = scalloc( system->N_cm_max * 2, sizeof( real ),
+ "Init_Workspace::workspace->rn[i]" );
+ workspace->v[i] = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->v[i]" );
+ }
+
+ workspace->r = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->r" );
+ workspace->d = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->d" );
+ workspace->q = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->q" );
+ workspace->p = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->p" );
+ break;
+
+ case CG_S:
+ workspace->r = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->r" );
+ workspace->d = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->d" );
+ workspace->q = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->q" );
+ workspace->p = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->p" );
+ break;
+
+ case SDM_S:
+ workspace->r = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->r" );
+ workspace->d = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->d" );
+ workspace->q = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->q" );
+ break;
+
+ case BiCGStab_S:
+ workspace->r = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->r" );
+ workspace->r_hat = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->r_hat" );
+ workspace->d = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->d" );
+ workspace->q = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->q" );
+ workspace->q_hat = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->q_hat" );
+ workspace->p = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->p" );
+ workspace->y = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->y" );
+ workspace->z = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->z" );
+ workspace->g = scalloc( system->N_cm_max, sizeof( real ),
+ "Init_Workspace::workspace->g" );
+ break;
+
+ default:
+ fprintf( stderr, "[ERROR] Unknown charge method linear solver type. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
- /* SpMV related */
#if defined(_OPENMP)
- workspace->b_local = smalloc( control->num_threads * system->N_cm * sizeof(real),
- "Init_Workspace::b_local" );
+ /* SpMV related */
+ workspace->b_local = smalloc( control->num_threads * system->N_cm_max * sizeof(real),
+ "Init_Workspace::b_local" );
#endif
+ }
/* level scheduling related */
workspace->levels_L = 1;
workspace->levels_U = 1;
- if ( control->cm_solver_pre_app_type == TRI_SOLVE_LEVEL_SCHED_PA ||
- control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
+ if ( alloc == TRUE )
{
- workspace->row_levels_L = smalloc( system->N_cm * sizeof(unsigned int),
- "Init_Workspace::row_levels_L" );
- workspace->level_rows_L = smalloc( system->N_cm * sizeof(unsigned int),
- "Init_Workspace::level_rows_L" );
- workspace->level_rows_cnt_L = smalloc( (system->N_cm + 1) * sizeof(unsigned int),
- "Init_Workspace::level_rows_cnt_L" );
- workspace->row_levels_U = smalloc( system->N_cm * sizeof(unsigned int),
- "Init_Workspace::row_levels_U" );
- workspace->level_rows_U = smalloc( system->N_cm * sizeof(unsigned int),
- "Init_Workspace::level_rows_U" );
- workspace->level_rows_cnt_U = smalloc( (system->N_cm + 1) * sizeof(unsigned int),
- "Init_Workspace::level_rows_cnt_U" );
- workspace->top = smalloc( (system->N_cm + 1) * sizeof(unsigned int),
- "Init_Workspace::top" );
- }
- else
- {
- workspace->row_levels_L = NULL;
- workspace->level_rows_L = NULL;
- workspace->level_rows_cnt_L = NULL;
- workspace->row_levels_U = NULL;
- workspace->level_rows_U = NULL;
- workspace->level_rows_cnt_U = NULL;
- workspace->top = NULL;
+ if ( control->cm_solver_pre_app_type == TRI_SOLVE_LEVEL_SCHED_PA ||
+ control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
+ {
+ workspace->row_levels_L = smalloc( system->N_cm_max * sizeof(unsigned int),
+ "Init_Workspace::row_levels_L" );
+ workspace->level_rows_L = smalloc( system->N_cm_max * sizeof(unsigned int),
+ "Init_Workspace::level_rows_L" );
+ workspace->level_rows_cnt_L = smalloc( (system->N_cm_max + 1) * sizeof(unsigned int),
+ "Init_Workspace::level_rows_cnt_L" );
+ workspace->row_levels_U = smalloc( system->N_cm_max * sizeof(unsigned int),
+ "Init_Workspace::row_levels_U" );
+ workspace->level_rows_U = smalloc( system->N_cm_max * sizeof(unsigned int),
+ "Init_Workspace::level_rows_U" );
+ workspace->level_rows_cnt_U = smalloc( (system->N_cm_max + 1) * sizeof(unsigned int),
+ "Init_Workspace::level_rows_cnt_U" );
+ workspace->top = smalloc( (system->N_cm_max + 1) * sizeof(unsigned int),
+ "Init_Workspace::top" );
+ }
+ else
+ {
+ workspace->row_levels_L = NULL;
+ workspace->level_rows_L = NULL;
+ workspace->level_rows_cnt_L = NULL;
+ workspace->row_levels_U = NULL;
+ workspace->level_rows_U = NULL;
+ workspace->level_rows_cnt_U = NULL;
+ workspace->top = NULL;
+ }
}
/* graph coloring related */
workspace->recolor_cnt = 0;
- if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
- {
- workspace->color = smalloc( sizeof(unsigned int) * system->N_cm,
- "Init_Workspace::color" );
- workspace->to_color = smalloc( sizeof(unsigned int) * system->N_cm,
- "Init_Workspace::to_color" );
- workspace->conflict = smalloc( sizeof(unsigned int) * system->N_cm,
- "setup_graph_coloring::conflict" );
- workspace->conflict_cnt = smalloc( sizeof(unsigned int) * (control->num_threads + 1),
- "Init_Workspace::conflict_cnt" );
- workspace->recolor = smalloc( sizeof(unsigned int) * system->N_cm,
- "Init_Workspace::recolor" );
- workspace->color_top = smalloc( sizeof(unsigned int) * (system->N_cm + 1),
- "Init_Workspace::color_top" );
- workspace->permuted_row_col = smalloc( sizeof(unsigned int) * system->N_cm,
- "Init_Workspace::premuted_row_col" );
- workspace->permuted_row_col_inv = smalloc( sizeof(unsigned int) * system->N_cm,
- "Init_Workspace::premuted_row_col_inv" );
- }
- else
+ if ( alloc == TRUE )
{
- workspace->color = NULL;
- workspace->to_color = NULL;
- workspace->conflict = NULL;
- workspace->conflict_cnt = NULL;
- workspace->recolor = NULL;
- workspace->color_top = NULL;
- workspace->permuted_row_col = NULL;
- workspace->permuted_row_col_inv = NULL;
- }
+ if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
+ {
+ workspace->color = smalloc( sizeof(unsigned int) * system->N_cm_max,
+ "Init_Workspace::color" );
+ workspace->to_color = smalloc( sizeof(unsigned int) * system->N_cm_max,
+ "Init_Workspace::to_color" );
+ workspace->conflict = smalloc( sizeof(unsigned int) * system->N_cm_max,
+ "setup_graph_coloring::conflict" );
+ workspace->conflict_cnt = smalloc( sizeof(unsigned int) * (control->num_threads + 1),
+ "Init_Workspace::conflict_cnt" );
+ workspace->recolor = smalloc( sizeof(unsigned int) * system->N_cm_max,
+ "Init_Workspace::recolor" );
+ workspace->color_top = smalloc( sizeof(unsigned int) * (system->N_cm_max + 1),
+ "Init_Workspace::color_top" );
+ workspace->permuted_row_col = smalloc( sizeof(unsigned int) * system->N_cm_max,
+ "Init_Workspace::premuted_row_col" );
+ workspace->permuted_row_col_inv = smalloc( sizeof(unsigned int) * system->N_cm_max,
+ "Init_Workspace::premuted_row_col_inv" );
+ }
+ else
+ {
+ workspace->color = NULL;
+ workspace->to_color = NULL;
+ workspace->conflict = NULL;
+ workspace->conflict_cnt = NULL;
+ workspace->recolor = NULL;
+ workspace->color_top = NULL;
+ workspace->permuted_row_col = NULL;
+ workspace->permuted_row_col_inv = NULL;
+ }
- /* graph coloring related OR ILUTP preconditioner */
- if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA
- || control->cm_solver_pre_comp_type == ILUTP_PC )
- {
- workspace->y_p = smalloc( sizeof(real) * system->N_cm, "Init_Workspace::y_p" );
- workspace->x_p = smalloc( sizeof(real) * system->N_cm, "Init_Workspace::x_p" );
- }
- else
- {
- workspace->y_p = NULL;
- workspace->x_p = NULL;
- }
+ /* graph coloring related OR ILUTP preconditioner */
+ if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA
+ || control->cm_solver_pre_comp_type == ILUTP_PC )
+ {
+ workspace->y_p = smalloc( sizeof(real) * system->N_cm_max, "Init_Workspace::y_p" );
+ workspace->x_p = smalloc( sizeof(real) * system->N_cm_max, "Init_Workspace::x_p" );
+ }
+ else
+ {
+ workspace->y_p = NULL;
+ workspace->x_p = NULL;
+ }
- /* Jacobi iteration related */
- if ( control->cm_solver_pre_app_type == JACOBI_ITER_PA )
- {
- workspace->Dinv_L = smalloc( sizeof(real) * system->N_cm,
- "Init_Workspace::Dinv_L" );
- workspace->Dinv_U = smalloc( sizeof(real) * system->N_cm,
- "Init_Workspace::Dinv_U" );
- workspace->Dinv_b = smalloc( sizeof(real) * system->N_cm,
- "Init_Workspace::Dinv_b" );
- workspace->rp = smalloc( sizeof(real) * system->N_cm,
- "Init_Workspace::rp" );
- workspace->rp2 = smalloc( sizeof(real) * system->N_cm,
- "Init_Workspace::rp2" );
- }
- else
- {
- workspace->Dinv_L = NULL;
- workspace->Dinv_U = NULL;
- workspace->Dinv_b = NULL;
- workspace->rp = NULL;
- workspace->rp2 = NULL;
- }
+ /* Jacobi iteration related */
+ if ( control->cm_solver_pre_app_type == JACOBI_ITER_PA )
+ {
+ workspace->Dinv_L = smalloc( sizeof(real) * system->N_cm_max,
+ "Init_Workspace::Dinv_L" );
+ workspace->Dinv_U = smalloc( sizeof(real) * system->N_cm_max,
+ "Init_Workspace::Dinv_U" );
+ workspace->Dinv_b = smalloc( sizeof(real) * system->N_cm_max,
+ "Init_Workspace::Dinv_b" );
+ workspace->rp = smalloc( sizeof(real) * system->N_cm_max,
+ "Init_Workspace::rp" );
+ workspace->rp2 = smalloc( sizeof(real) * system->N_cm_max,
+ "Init_Workspace::rp2" );
+ }
+ else
+ {
+ workspace->Dinv_L = NULL;
+ workspace->Dinv_U = NULL;
+ workspace->Dinv_b = NULL;
+ workspace->rp = NULL;
+ workspace->rp2 = NULL;
+ }
- /* ILUTP preconditioner related */
- if ( control->cm_solver_pre_comp_type == ILUTP_PC )
- {
- workspace->perm_ilutp = smalloc( sizeof( int ) * system->N_cm,
- "Init_Workspace::workspace->perm_ilutp" );
- }
- else
- {
- workspace->perm_ilutp = NULL;
- }
+ /* ILUTP preconditioner related */
+ if ( control->cm_solver_pre_comp_type == ILUTP_PC )
+ {
+ workspace->perm_ilutp = smalloc( sizeof( int ) * system->N_cm_max,
+ "Init_Workspace::workspace->perm_ilutp" );
+ }
+ else
+ {
+ workspace->perm_ilutp = NULL;
+ }
- /* integrator storage */
- workspace->a = smalloc( system->N * sizeof( rvec ),
- "Init_Workspace::workspace->a" );
- workspace->f_old = smalloc( system->N * sizeof( rvec ),
- "Init_Workspace::workspace->f_old" );
- workspace->v_const = smalloc( system->N * sizeof( rvec ),
- "Init_Workspace::workspace->v_const" );
+ /* integrator storage */
+ workspace->a = smalloc( system->N_max * sizeof( rvec ),
+ "Init_Workspace::workspace->a" );
+ workspace->f_old = smalloc( system->N_max * sizeof( rvec ),
+ "Init_Workspace::workspace->f_old" );
+ workspace->v_const = smalloc( system->N_max * sizeof( rvec ),
+ "Init_Workspace::workspace->v_const" );
#if defined(_OPENMP)
- workspace->f_local = smalloc( control->num_threads * system->N * sizeof( rvec ),
- "Init_Workspace::workspace->f_local" );
+ workspace->f_local = smalloc( control->num_threads * system->N_max * sizeof( rvec ),
+ "Init_Workspace::workspace->f_local" );
#endif
- /* storage for analysis */
- if ( control->molec_anal || control->diffusion_coef )
- {
- workspace->mark = scalloc( system->N, sizeof(int),
- "Init_Workspace::workspace->mark" );
- workspace->old_mark = scalloc( system->N, sizeof(int),
- "Init_Workspace::workspace->old_mark" );
- }
- else
- {
- workspace->mark = workspace->old_mark = NULL;
- }
+ /* storage for analysis */
+ if ( control->molec_anal || control->diffusion_coef )
+ {
+ workspace->mark = scalloc( system->N_max, sizeof(int),
+ "Init_Workspace::workspace->mark" );
+ workspace->old_mark = scalloc( system->N_max, sizeof(int),
+ "Init_Workspace::workspace->old_mark" );
+ }
+ else
+ {
+ workspace->mark = workspace->old_mark = NULL;
+ }
- if ( control->diffusion_coef )
- {
- workspace->x_old = scalloc( system->N, sizeof( rvec ),
- "Init_Workspace::workspace->x_old" );
- }
- else
- {
- workspace->x_old = NULL;
+ if ( control->diffusion_coef )
+ {
+ workspace->x_old = scalloc( system->N_max, sizeof( rvec ),
+ "Init_Workspace::workspace->x_old" );
+ }
+ else
+ {
+ workspace->x_old = NULL;
+ }
}
#if defined(TEST_FORCES)
- workspace->dDelta = smalloc( system->N * sizeof( rvec ),
+ workspace->dDelta = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->dDelta" );
- workspace->f_ele = smalloc( system->N * sizeof( rvec ),
+ workspace->f_ele = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_ele" );
- workspace->f_vdw = smalloc( system->N * sizeof( rvec ),
+ workspace->f_vdw = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_vdw" );
- workspace->f_be = smalloc( system->N * sizeof( rvec ),
+ workspace->f_be = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_be" );
- workspace->f_lp = smalloc( system->N * sizeof( rvec ),
+ workspace->f_lp = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_lp" );
- workspace->f_ov = smalloc( system->N * sizeof( rvec ),
+ workspace->f_ov = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_ov" );
- workspace->f_un = smalloc( system->N * sizeof( rvec ),
+ workspace->f_un = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_un" );
- workspace->f_ang = smalloc( system->N * sizeof( rvec ),
+ workspace->f_ang = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_ang" );
- workspace->f_coa = smalloc( system->N * sizeof( rvec ),
+ workspace->f_coa = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_coa" );
- workspace->f_pen = smalloc( system->N * sizeof( rvec ),
+ workspace->f_pen = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_pen" );
- workspace->f_hb = smalloc( system->N * sizeof( rvec ),
+ workspace->f_hb = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_hb" );
- workspace->f_tor = smalloc( system->N * sizeof( rvec ),
+ workspace->f_tor = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_tor" );
- workspace->f_con = smalloc( system->N * sizeof( rvec ),
+ workspace->f_con = smalloc( system->N_max * sizeof( rvec ),
"Init_Workspace::workspace->f_con" );
#endif
@@ -756,22 +774,27 @@ static void Init_Workspace( reax_system *system, control_params *control,
static void Init_Lists( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+ reax_list **lists, output_controls *out_control,
+ int first_run, int alloc )
{
- int i, num_nbrs, num_bonds, num_3body, Htop, max_nnz;
+ int i, num_nbrs, num_bonds, num_hbonds, num_3body, Htop, max_nnz;
int *hb_top, *bond_top;
-#if defined(DEBUG_FOCUS)
- int num_hbonds;
-#endif
num_nbrs = Estimate_Num_Neighbors( system, control, workspace, lists );
- Make_List( system->N, num_nbrs, TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "memory allocated: far_nbrs = %ldMB\n",
- num_nbrs * sizeof(far_neighbor_data) / (1024 * 1024) );
-#endif
+ if ( first_run == TRUE )
+ {
+ Make_List( system->N, system->N_max, num_nbrs, TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
+ }
+ else if ( alloc == TRUE || lists[FAR_NBRS]->total_intrs < num_nbrs )
+ {
+ Delete_List( TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
+ Make_List( system->N, system->N_max, num_nbrs, TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
+ }
+ else
+ {
+ lists[FAR_NBRS]->n = system->N;
+ }
Generate_Neighbor_Lists( system, control, data, workspace, lists, out_control );
@@ -800,12 +823,41 @@ static void Init_Lists( reax_system *system, control_params *control,
break;
}
- Allocate_Matrix( &workspace->H, system->N_cm, max_nnz );
- /* TODO: better estimate for H_sp?
- * If so, need to refactor Estimate_Storage_Sizes
- * to use various cut-off distances as parameters
- * (non-bonded, hydrogen, 3body, etc.) */
- Allocate_Matrix( &workspace->H_sp, system->N_cm, max_nnz );
+ if ( first_run == TRUE )
+ {
+ Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max, max_nnz );
+ }
+ else if ( alloc == TRUE || workspace->H->m < max_nnz )
+ {
+ Deallocate_Matrix( workspace->H );
+ Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max, max_nnz );
+ }
+ else
+ {
+ workspace->H->n = system->N_cm;
+ }
+
+ if ( first_run == TRUE )
+ {
+ /* TODO: better estimate for H_sp?
+ * If so, need to refactor Estimate_Storage_Sizes
+ * to use various cut-off distances as parameters
+ * (non-bonded, hydrogen, 3body, etc.) */
+ Allocate_Matrix( &workspace->H_sp, system->N_cm, system->N_cm_max, max_nnz );
+ }
+ else if ( alloc == TRUE || workspace->H_sp->m < max_nnz )
+ {
+ Deallocate_Matrix( workspace->H_sp );
+ /* TODO: better estimate for H_sp?
+ * If so, need to refactor Estimate_Storage_Sizes
+ * to use various cut-off distances as parameters
+ * (non-bonded, hydrogen, 3body, etc.) */
+ Allocate_Matrix( &workspace->H_sp, system->N_cm, system->N_cm_max, max_nnz );
+ }
+ else
+ {
+ workspace->H_sp->n = system->N_cm;
+ }
workspace->num_H = 0;
if ( control->hbond_cut > 0.0 )
@@ -829,21 +881,64 @@ static void Init_Lists( reax_system *system, control_params *control,
}
else
{
- Allocate_HBond_List( system->N, workspace->num_H, workspace->hbond_index,
- hb_top, lists[HBONDS] );
+ num_hbonds = 0;
+ for ( i = 0; i < system->N; ++i )
+ {
+ num_hbonds += hb_top[i];
+ }
+
+ if ( first_run == TRUE )
+ {
+ workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
+
+ Allocate_HBond_List( system->N, workspace->num_H, workspace->num_H_max,
+ workspace->hbond_index, hb_top, lists[HBONDS] );
+ }
+ else if ( workspace->num_H_max < workspace->num_H
+ || lists[HBONDS]->total_intrs < num_hbonds )
+ {
+ if ( workspace->num_H_max < workspace->num_H )
+ {
+ workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
+ }
+
+ Delete_List( TYP_HBOND, lists[HBONDS] );
+ Allocate_HBond_List( system->N, workspace->num_H, workspace->num_H_max,
+ workspace->hbond_index, hb_top, lists[HBONDS] );
+ }
+ else
+ {
+ lists[HBONDS]->n = workspace->num_H;
+ }
}
#if defined(DEBUG_FOCUS)
- num_hbonds = hb_top[system->N - 1];
fprintf( stderr, "estimated storage - num_hbonds: %d\n", num_hbonds );
fprintf( stderr, "memory allocated: hbonds = %ldMB\n",
num_hbonds * sizeof(hbond_data) / (1024 * 1024) );
#endif
}
+ num_bonds = 0;
+ for ( i = 0; i < system->N; ++i )
+ {
+ num_bonds += bond_top[i];
+ }
+
/* bonds list */
- Allocate_Bond_List( system->N, bond_top, lists[BONDS] );
- num_bonds = bond_top[system->N - 1];
+ if ( first_run == TRUE )
+ {
+ Allocate_Bond_List( system->N, system->N_max, bond_top, lists[BONDS] );
+ }
+ else if ( alloc == TRUE || lists[BONDS]->total_intrs < num_bonds )
+ {
+ Delete_List( TYP_BOND, lists[BONDS] );
+ Allocate_Bond_List( system->N, system->N_max, bond_top, lists[BONDS] );
+ }
+ else
+ {
+ lists[BONDS]->n = system->N;
+ }
#if defined(DEBUG_FOCUS)
fprintf( stderr, "estimated storage - num_bonds: %d\n", num_bonds );
@@ -852,7 +947,29 @@ static void Init_Lists( reax_system *system, control_params *control,
#endif
/* 3bodies list */
- Make_List( num_bonds, num_3body, TYP_THREE_BODY, lists[THREE_BODIES] );
+ if ( first_run == TRUE )
+ {
+ Make_List( num_bonds, num_bonds, num_3body, TYP_THREE_BODY, lists[THREE_BODIES] );
+ }
+ else if ( lists[THREE_BODIES]->n_max < num_bonds
+ || lists[THREE_BODIES]->total_intrs < num_3body )
+ {
+ Delete_List( TYP_THREE_BODY, lists[THREE_BODIES] );
+ if ( lists[THREE_BODIES]->n_max < num_bonds )
+ {
+ Make_List( num_bonds, num_bonds, num_3body,
+ TYP_THREE_BODY, lists[THREE_BODIES] );
+ }
+ else
+ {
+ Make_List( num_bonds, lists[THREE_BODIES]->n_max, num_3body,
+ TYP_THREE_BODY, lists[THREE_BODIES] );
+ }
+ }
+ else
+ {
+ lists[THREE_BODIES]->n = num_bonds;
+ }
#if defined(DEBUG_FOCUS)
fprintf( stderr, "estimated storage - num_3body: %d\n", num_3body );
@@ -885,12 +1002,12 @@ static void Init_Lists( reax_system *system, control_params *control,
static void Init_Out_Controls( reax_system *system, control_params *control,
- static_storage *workspace, output_controls *out_control )
+ static_storage *workspace, output_controls *out_control, int output_enabled )
{
#define TEMP_SIZE (1000)
char temp[TEMP_SIZE];
- if ( out_control->write_steps > 0 )
+ if ( output_enabled == TRUE && out_control->write_steps > 0 )
{
strncpy( temp, control->sim_name, TEMP_SIZE - 5 );
temp[TEMP_SIZE - 5] = '\0';
@@ -903,7 +1020,7 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
out_control->trj = NULL;
}
- if ( out_control->log_update_freq > 0 )
+ if ( output_enabled == TRUE && out_control->log_update_freq > 0 )
{
strncpy( temp, control->sim_name, TEMP_SIZE - 5 );
temp[TEMP_SIZE - 5] = '\0';
@@ -940,8 +1057,8 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
out_control->log = NULL;
}
- if ( control->ensemble == sNPT || control->ensemble == iNPT
- || control->ensemble == aNPT || control->compute_pressure == TRUE )
+ if ( output_enabled == TRUE && (control->ensemble == sNPT || control->ensemble == iNPT
+ || control->ensemble == aNPT || control->compute_pressure == TRUE) )
{
strncpy( temp, control->sim_name, TEMP_SIZE - 5 );
temp[TEMP_SIZE - 5] = '\0';
@@ -964,7 +1081,7 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
}
/* Init molecular analysis file */
- if ( control->molec_anal )
+ if ( output_enabled == TRUE && control->molec_anal )
{
snprintf( temp, TEMP_SIZE, "%.*s.mol", TEMP_SIZE - 5, control->sim_name );
out_control->mol = sfopen( temp, "w" );
@@ -980,7 +1097,7 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
out_control->ign = NULL;
}
- if ( control->dipole_anal )
+ if ( output_enabled == TRUE && control->dipole_anal )
{
strncpy( temp, control->sim_name, TEMP_SIZE - 5 );
temp[TEMP_SIZE - 5] = '\0';
@@ -995,7 +1112,7 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
out_control->dpl = NULL;
}
- if ( control->diffusion_coef )
+ if ( output_enabled == TRUE && control->diffusion_coef )
{
strncpy( temp, control->sim_name, TEMP_SIZE - 6 );
temp[TEMP_SIZE - 6] = '\0';
@@ -1148,12 +1265,8 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
void Initialize( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace, reax_list **lists,
output_controls *out_control, evolve_function *Evolve,
- const int output_enabled )
+ int output_enabled, int first_run, int alloc )
{
-#if defined(DEBUG)
- real start, end;
-#endif
-
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(control)
{
@@ -1166,18 +1279,15 @@ void Initialize( reax_system *system, control_params *control,
Randomize( );
- Init_System( system, control, data );
+ Init_System( system, control, data, first_run );
- Init_Simulation_Data( system, control, data, out_control, Evolve );
+ Init_Simulation_Data( system, control, data, out_control, Evolve, alloc );
- Init_Workspace( system, control, workspace );
+ Init_Workspace( system, control, workspace, alloc );
- Init_Lists( system, control, data, workspace, lists, out_control );
+ Init_Lists( system, control, data, workspace, lists, out_control, first_run, alloc );
- if ( output_enabled == TRUE )
- {
- Init_Out_Controls( system, control, workspace, out_control );
- }
+ Init_Out_Controls( system, control, workspace, out_control, output_enabled );
/* These are done in forces.c, only forces.c can see all those functions */
Init_Bonded_Force_Functions( control );
@@ -1188,27 +1298,13 @@ void Initialize( reax_system *system, control_params *control,
if ( control->tabulate )
{
-#if defined(DEBUG)
- start = Get_Time( );
-#endif
-
Make_LR_Lookup_Table( system, control, workspace );
-
-#if defined(DEBUG)
- end = Get_Timing_Info( start );
-
- fprintf( stderr, "[INFO] LR lookup table calculation time: %f\n", end );
-#endif
}
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "data structures have been initialized...\n" );
-#endif
}
static void Finalize_System( reax_system *system, control_params *control,
- simulation_data *data )
+ simulation_data *data, int reset )
{
int i, j, k;
reax_interaction *reax;
@@ -1217,35 +1313,38 @@ static void Finalize_System( reax_system *system, control_params *control,
Finalize_Grid( system );
- sfree( reax->gp.l, "Finalize_System::reax->gp.l" );
-
- for ( i = 0; i < reax->num_atom_types; i++ )
+ if ( reset == FALSE )
{
- for ( j = 0; j < reax->num_atom_types; j++ )
+ sfree( reax->gp.l, "Finalize_System::reax->gp.l" );
+
+ for ( i = 0; i < reax->num_atom_types; i++ )
{
- for ( k = 0; k < reax->num_atom_types; k++ )
+ for ( j = 0; j < reax->num_atom_types; j++ )
{
- sfree( reax->fbp[i][j][k], "Finalize_System::reax->fbp[i][j][k]" );
+ for ( k = 0; k < reax->num_atom_types; k++ )
+ {
+ sfree( reax->fbp[i][j][k], "Finalize_System::reax->fbp[i][j][k]" );
+ }
+
+ 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]" );
}
- 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]" );
+ 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]" );
}
- 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]" );
- }
-
- 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" );
+ 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" );
- sfree( system->atoms, "Finalize_System::system->atoms" );
+ sfree( system->atoms, "Finalize_System::system->atoms" );
+ }
}
@@ -1256,14 +1355,14 @@ static void Finalize_Simulation_Data( reax_system *system, control_params *contr
if ( control->ensemble == sNPT || control->ensemble == iNPT
|| control->ensemble == aNPT || control->compute_pressure == TRUE )
{
- sfree( data->press_local, "Finalize_Workspace::workspace->press_local" );
+ sfree( data->press_local, "Finalize_Simulation_Data::data->press_local" );
}
#endif
}
static void Finalize_Workspace( reax_system *system, control_params *control,
- static_storage *workspace )
+ static_storage *workspace, int reset )
{
int i;
@@ -1285,9 +1384,9 @@ static void Finalize_Workspace( reax_system *system, control_params *control,
sfree( workspace->CdDelta, "Finalize_Workspace::workspace->CdDelta" );
sfree( workspace->vlpex, "Finalize_Workspace::workspace->vlpex" );
- if ( control->geo_format == BGF
+ if ( reset == FALSE && (control->geo_format == BGF
|| control->geo_format == ASCII_RESTART
- || control->geo_format == BINARY_RESTART )
+ || control->geo_format == BINARY_RESTART) )
{
sfree( workspace->map_serials, "Finalize_Workspace::workspace->map_serials" );
}
@@ -1472,19 +1571,22 @@ static void Finalize_Workspace( reax_system *system, control_params *control,
sfree( workspace->x_old, "Finalize_Workspace::workspace->x_old" );
}
- sfree( workspace->orig_id, "Finalize_Workspace::workspace->orig_id" );
-
- /* space for keeping restriction info, if any */
- if ( control->restrict_bonds )
+ if ( reset == FALSE )
{
- for ( i = 0; i < system->N; ++i )
+ sfree( workspace->orig_id, "Finalize_Workspace::workspace->orig_id" );
+
+ /* space for keeping restriction info, if any */
+ if ( control->restrict_bonds )
{
- sfree( workspace->restricted_list[i],
- "Finalize_Workspace::workspace->restricted_list[i]" );
- }
+ for ( i = 0; i < system->N; ++i )
+ {
+ sfree( workspace->restricted_list[i],
+ "Finalize_Workspace::workspace->restricted_list[i]" );
+ }
- sfree( workspace->restricted, "Finalize_Workspace::workspace->restricted" );
- sfree( workspace->restricted_list, "Finalize_Workspace::workspace->restricted_list" );
+ sfree( workspace->restricted, "Finalize_Workspace::workspace->restricted" );
+ sfree( workspace->restricted_list, "Finalize_Workspace::workspace->restricted_list" );
+ }
}
#if defined(TEST_FORCES)
@@ -1600,7 +1702,7 @@ static void Finalize_Out_Controls( reax_system *system, control_params *control,
*/
void Finalize( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace, reax_list **lists,
- output_controls *out_control, const int output_enabled )
+ output_controls *out_control, int output_enabled, int reset )
{
if ( control->tabulate )
{
@@ -1614,9 +1716,9 @@ void Finalize( reax_system *system, control_params *control,
Finalize_Lists( control, lists );
- Finalize_Workspace( system, control, workspace );
+ Finalize_Workspace( system, control, workspace, reset );
Finalize_Simulation_Data( system, control, data, out_control );
- Finalize_System( system, control, data );
+ Finalize_System( system, control, data, reset );
}
diff --git a/sPuReMD/src/init_md.h b/sPuReMD/src/init_md.h
index 98d5cc90902422075a705c71de6508df5c8a36e2..26b3b0d5229f4111964c50af27b966661d1afcd8 100644
--- a/sPuReMD/src/init_md.h
+++ b/sPuReMD/src/init_md.h
@@ -27,10 +27,10 @@
void Initialize( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls*, evolve_function*,
- const int );
+ int, int, int );
void Finalize( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls*, const int );
+ static_storage*, reax_list**, output_controls*, int, int );
#endif
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 9fc07f2d3efbf8e0b37c8ac996261df31a329bd1..44ddce77c8107bf7363809d72cf713eb15937901 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -57,7 +57,7 @@ static sparse_matrix * create_test_mat( void )
unsigned int i, n;
sparse_matrix *H_test;
- Allocate_Matrix( &H_test, 3, 6 );
+ Allocate_Matrix( &H_test, 3, 3, 6 );
//3x3, SPD, store lower half
i = 0;
@@ -169,15 +169,15 @@ static void compute_full_sparse_matrix( const sparse_matrix * const A,
if ( *A_full == NULL )
{
- Allocate_Matrix( A_full, A->n, 2 * A->m - A->n );
+ Allocate_Matrix( A_full, A->n, A->n_max, 2 * A->m - A->n );
}
else if ( (*A_full)->m < 2 * A->m - A->n )
{
Deallocate_Matrix( *A_full );
- Allocate_Matrix( A_full, A->n, 2 * A->m - A->n );
+ Allocate_Matrix( A_full, A->n, A->n_max, 2 * A->m - A->n );
}
- Allocate_Matrix( &A_t, A->n, A->m );
+ Allocate_Matrix( &A_t, A->n, A->n_max, A->m );
/* Set up the sparse matrix data structure for A. */
Transpose( A, A_t );
@@ -233,12 +233,12 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
if ( *A_spar_patt == NULL )
{
- Allocate_Matrix( A_spar_patt, A->n, A->m );
+ Allocate_Matrix( A_spar_patt, A->n, A->n_max, A->m );
}
else if ( ((*A_spar_patt)->m) < (A->m) )
{
Deallocate_Matrix( *A_spar_patt );
- Allocate_Matrix( A_spar_patt, A->n, A->m );
+ Allocate_Matrix( A_spar_patt, A->n, A->n_max, A->m );
}
list = smalloc( sizeof(real) * A->start[A->n],
@@ -352,7 +352,8 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
{
/* A_app_inv has the same sparsity pattern
* as A_spar_patt_full (omit non-zero values) */
- Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m );
+ Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->n_max,
+ (*A_spar_patt_full)->m );
}
else if ( (*A_app_inv)->m < A->m )
{
@@ -360,7 +361,8 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
/* A_app_inv has the same sparsity pattern
* as A_spar_patt_full (omit non-zero values) */
- Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m );
+ Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->n_max,
+ (*A_spar_patt_full)->m );
}
sfree( list, "setup_sparse_approx_inverse::list" );
@@ -1130,8 +1132,8 @@ real FG_ICHOLT( const sparse_matrix * const A, const real * droptol,
start = Get_Time( );
- Allocate_Matrix( &DAD, A->n, A->m );
- Allocate_Matrix( &U_T_temp, A->n, A->m );
+ Allocate_Matrix( &DAD, A->n, A->n_max, A->m );
+ Allocate_Matrix( &U_T_temp, A->n, A->n_max, A->m );
D = smalloc( sizeof(real) * A->n, "FG_ICHOLT::D" );
D_inv = smalloc( sizeof(real) * A->n, "FG_ICHOLT::D_inv" );
@@ -1338,9 +1340,9 @@ real FG_ILUT( const sparse_matrix * const A, const real * droptol,
start = Get_Time( );
- Allocate_Matrix( &DAD, A->n, A->m );
- Allocate_Matrix( &L_temp, A->n, A->m );
- Allocate_Matrix( &U_T_temp, A->n, A->m );
+ Allocate_Matrix( &DAD, A->n, A->n_max, A->m );
+ Allocate_Matrix( &L_temp, A->n, A->n_max, A->m );
+ Allocate_Matrix( &U_T_temp, A->n, A->n_max, A->m );
D = smalloc( sizeof(real) * A->n, "FG_ILUT::D" );
D_inv = smalloc( sizeof(real) * A->n, "FG_ILUT::D_inv" );
@@ -1980,7 +1982,7 @@ void Transpose_I( sparse_matrix * const A )
{
sparse_matrix * A_t;
- Allocate_Matrix( &A_t, A->n, A->m );
+ Allocate_Matrix( &A_t, A->n, A->n_max, A->m );
Transpose( A, A_t );
@@ -2487,7 +2489,7 @@ static void permute_matrix( const static_storage * const workspace,
int i, pj, nr, nc;
sparse_matrix *LUtemp;
- Allocate_Matrix( &LUtemp, LU->n, LU->m );
+ Allocate_Matrix( &LUtemp, LU->n, LU->n_max, LU->m );
for ( i = 0; i < LU->n + 1; ++i )
workspace->color_top[i] = 0;
@@ -2623,12 +2625,12 @@ void setup_graph_coloring( const control_params * const control,
{
if ( *H_p == NULL )
{
- Allocate_Matrix( H_p, H->n, H->m );
+ Allocate_Matrix( H_p, H->n, H->n_max, H->m );
}
else if ( (*H_p)->m < H->m )
{
Deallocate_Matrix( *H_p );
- Allocate_Matrix( H_p, H->n, H->m );
+ Allocate_Matrix( H_p, H->n, H->n_max, H->m );
}
compute_full_sparse_matrix( H, H_full );
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index ad7bbc626fa343eadda8f481274f7d9fde0b0aed..f5414e93a3424e97364dda0382e201eb82ad1a72 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -24,17 +24,20 @@
#include "tool_box.h"
-void Make_List( int n, int total_intrs, int type, reax_list* l )
+void Make_List( int n, int n_max, int total_intrs, int type, reax_list* l )
{
assert( n > 0 );
+ assert( n_max > 0 );
+ assert( n_max >= n );
assert( total_intrs > 0 );
assert( l != NULL );
l->n = n;
+ l->n_max = n_max;
l->total_intrs = total_intrs;
- l->index = smalloc( n * sizeof(int), "Make_List::l->index" );
- l->end_index = smalloc( n * sizeof(int), "Make_List::l->end_index" );
+ l->index = smalloc( n_max * sizeof(int), "Make_List::l->index" );
+ l->end_index = smalloc( n_max * sizeof(int), "Make_List::l->end_index" );
switch ( type )
{
diff --git a/sPuReMD/src/list.h b/sPuReMD/src/list.h
index c83d8546b1386100bdaa69599757253b862cb2f1..3a0bcdcb39d1f57bd76eaaab3de8149a1f8935e1 100644
--- a/sPuReMD/src/list.h
+++ b/sPuReMD/src/list.h
@@ -25,7 +25,7 @@
#include "reax_types.h"
-void Make_List( int, int, int, reax_list* );
+void Make_List( int, int, int, int, reax_list* );
void Delete_List( int, reax_list* );
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 734ad9b5eb8cbcad32f5c03b331868d156a6b199..f1fac619f57aea003770c24365ea9b3749f42e04 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -786,14 +786,14 @@ struct grid
/* max. num. of neighbors for a grid cell;
* used for memory allocation purposes */
int max_nbrs;
- /**/
- int total;
/* lengths of each dimesion of a grid cell */
real cell_size;
/**/
ivec spread;
- /* num. of cells along each dimension for entire grid */
+ /* num. of cells along each dimension for entire grid for this simulation */
ivec ncell;
+ /* max. num. of cells along each dimension for entire grid across all simulations */
+ ivec ncell_max;
/* lengths of each dimesion of a grid cell */
rvec len;
/* multiplicative inverse of length of each dimesion of a grid cell */
@@ -820,10 +820,14 @@ struct grid
struct reax_system
{
- /* number of local (non-periodic image) atoms */
+ /* number of local (non-periodic image) atoms for the current simulation */
int N;
+ /* max. number of local (non-periodic image) atoms across all simulations */
+ int N_max;
/* dimension of the N x N sparse charge method matrix H */
int N_cm;
+ /* max. dimension of the N x N sparse charge method matrix H across all simulations */
+ int N_cm_max;
/* atom info */
reax_atom *atoms;
/* atomic interaction parameters */
@@ -1122,6 +1126,8 @@ struct reax_timing
struct simulation_data
{
+ /* integer ID uniquely identifying the simulation */
+ int sim_id;
/* current simulation step number (0-based) */
int step;
/* last simulation step number for restarted runs (0-based) */
@@ -1349,8 +1355,10 @@ struct bond_data
*/
struct sparse_matrix
{
- /* number of rows */
+ /* active number of rows for this simulation */
unsigned int n;
+ /* max. number of rows across all simulations */
+ unsigned int n_max;
/* number of nonzeros (NNZ) ALLOCATED */
unsigned int m;
/* row pointer (last element contains ACTUAL NNZ) */
@@ -1530,8 +1538,12 @@ struct static_storage
/* permutation for ILUTP */
unsigned int *perm_ilutp;
+ /* num. hydrogen atoms for this simulation */
int num_H;
- int *hbond_index; // for hydrogen bonds
+ /* max. num. hydrogen atoms across all simulations */
+ int num_H_max;
+ /* for hydrogen bonds */
+ int *hbond_index;
rvec *v_const;
rvec *f_old;
@@ -1586,15 +1598,17 @@ struct static_storage
/* interaction lists */
struct reax_list
{
- /* num. entries in list */
+ /* num. active entries in list for this simulation */
int n;
- /* sum of max. interactions per atom */
+ /* max. num. entries in list across all simulations */
+ int n_max;
+ /* total interactions across all entries which can be stored in the list */
int total_intrs;
/* starting position of atom's interactions */
int *index;
/* ending position of atom's interactions */
int *end_index;
- /* max. num. of interactions per atom */
+ /* max. num. of interactions per list entry */
int *max_intrs;
/* interaction list (polymorphic via union dereference) */
// union
@@ -1701,20 +1715,24 @@ struct output_controls
/* Handle for working with an instance of the sPuReMD library */
struct spuremd_handle
{
- /* System info. struct pointer */
+ /* System info. */
reax_system *system;
- /* Control parameters struct pointer */
+ /* Control parameters */
control_params *control;
- /* Atomic simulation data struct pointer */
+ /* Atomic simulation data */
simulation_data *data;
- /* Internal workspace struct pointer */
+ /* Internal workspace */
static_storage *workspace;
- /* Reax interaction list struct pointer */
+ /* Reax interaction list */
reax_list **lists;
- /* Output controls struct pointer */
+ /* Output controls */
output_controls *out_control;
+ /* TRUE for first simulation, FALSE otherwise */
+ int first_run;
/* TRUE if file I/O for simulation output enabled, FALSE otherwise */
int output_enabled;
+ /* TRUE if reallocation is required due to num. atoms increasing, FALSE otherwise */
+ int realloc;
/* Callback for getting simulation state at the end of each time step */
callback_function callback;
};
diff --git a/sPuReMD/src/restart.c b/sPuReMD/src/restart.c
index 0cdbcd4dbb184efcef6c4e80883c2c15986c2c50..32b3a62244b6cbdd1e62c30c0134edc21c393cc3 100644
--- a/sPuReMD/src/restart.c
+++ b/sPuReMD/src/restart.c
@@ -73,7 +73,7 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
void Read_Binary_Restart( const char * const fname, reax_system *system,
control_params *control, simulation_data *data,
- static_storage *workspace )
+ static_storage *workspace, int first_run )
{
int i;
FILE *fres;
@@ -87,7 +87,6 @@ void Read_Binary_Restart( const char * const fname, reax_system *system,
fread( &res_header, sizeof(restart_header), 1, fres );
data->prev_steps = res_header.step;
- system->N = res_header.N;
data->therm.T = res_header.T;
data->therm.xi = res_header.xi;
data->therm.v_xi = res_header.v_xi;
@@ -109,10 +108,17 @@ void Read_Binary_Restart( const char * const fname, reax_system *system,
system->box.box[2][0], system->box.box[2][1], system->box.box[2][2] );
#endif
- PreAllocate_Space( system, control, workspace );
+ if ( first_run == TRUE || res_header.N > system->N )
+ {
+ PreAllocate_Space( system, control, workspace, res_header.N, first_run );
+ }
+ system->N = res_header.N;
- workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
- "Read_Binary_Restart::workspace->map_serials" );
+ if ( first_run == TRUE )
+ {
+ workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
+ "Read_Binary_Restart::workspace->map_serials" );
+ }
for ( i = 0; i < MAX_ATOM_ID; ++i )
{
@@ -185,9 +191,9 @@ void Write_ASCII_Restart( reax_system *system, control_params *control,
void Read_ASCII_Restart( const char * const fname, reax_system *system,
control_params *control, simulation_data *data,
- static_storage *workspace )
+ static_storage *workspace, int first_run )
{
- int i;
+ int i, n;
FILE *fres;
reax_atom *p_atom;
@@ -195,7 +201,7 @@ void Read_ASCII_Restart( const char * const fname, reax_system *system,
/* parse header of restart file */
fscanf( fres, READ_RESTART_HEADER,
- &data->prev_steps, &system->N, &data->therm.T, &data->therm.xi,
+ &data->prev_steps, &n, &data->therm.T, &data->therm.xi,
&data->therm.v_xi, &data->therm.v_xi_old, &data->therm.G_xi,
&system->box.box[0][0], &system->box.box[0][1], &system->box.box[0][2],
&system->box.box[1][0], &system->box.box[1][1], &system->box.box[1][2],
@@ -216,10 +222,17 @@ void Read_ASCII_Restart( const char * const fname, reax_system *system,
system->box.box[2][0], system->box.box[2][1], system->box.box[2][2] );
#endif
- PreAllocate_Space( system, control, workspace );
+ if ( first_run == TRUE || n > system->N )
+ {
+ PreAllocate_Space( system, control, workspace, n, first_run );
+ }
+ system->N = n;
- workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
- "Read_ASCII_Restart::workspace->map_serials" );
+ if ( first_run == TRUE )
+ {
+ workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
+ "Read_ASCII_Restart::workspace->map_serials" );
+ }
for ( i = 0; i < MAX_ATOM_ID; ++i )
{
diff --git a/sPuReMD/src/restart.h b/sPuReMD/src/restart.h
index cfa36c5afaa05bebdddf4ba72cd8f9f354cf3361..41e221eaafc0c99b558b0116344cb83b45cecef8 100644
--- a/sPuReMD/src/restart.h
+++ b/sPuReMD/src/restart.h
@@ -53,12 +53,12 @@ typedef struct
#define READ_RESTART_HEADER " %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf"
#define READ_RESTART_LINE " %d %d %s %lf %lf %lf %lf %lf %lf"
-void Write_Restart( reax_system*, control_params*,
- simulation_data*, static_storage*, output_controls* );
+void Write_Restart( reax_system *, control_params *,
+ simulation_data *, static_storage *, output_controls * );
-void Read_Binary_Restart( const char * const, reax_system*, control_params*,
- simulation_data*, static_storage* );
-void Read_ASCII_Restart( const char * const, reax_system*, control_params*,
- simulation_data*, static_storage* );
+void Read_Binary_Restart( const char * const, reax_system *, control_params *,
+ simulation_data *, static_storage *, int );
+void Read_ASCII_Restart( const char * const, reax_system *, control_params *,
+ simulation_data *, static_storage *, int );
#endif
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index d009d160e52e11de7f6d654ed0faf78d8eef53a6..453b0263837bd685e61815365bcda06d445f2e7a 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -39,6 +39,9 @@
#include "vector.h"
+/* Handles additional entire geometry calculations after
+ * perturbing atom positions during a simulation step
+ */
static void Post_Evolve( reax_system * const system, control_params * const control,
simulation_data * const data, static_storage * const workspace,
reax_list ** const lists, output_controls * const out_control )
@@ -51,7 +54,7 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
&& data->step % control->remove_CoM_vel == 0 )
{
/* compute velocity of the center of mass */
- Compute_Center_of_Mass( system, data, out_control->prs );
+ Compute_Center_of_Mass( system, data );
for ( i = 0; i < system->N; i++ )
{
@@ -83,44 +86,47 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
}
-static void Read_System( const char * const geo_file,
- const char * const ffield_file,
- const char * const control_file,
- reax_system * const system,
- control_params * const control,
- simulation_data * const data,
- static_storage * const workspace,
- output_controls * const out_control )
+/* Parse input files
+ *
+ * geo_file: file containing geometry info of the structure to simulate
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
+static void Read_Input_Files( const char * const geo_file,
+ const char * const ffield_file, const char * const control_file,
+ reax_system * const system, control_params * const control,
+ simulation_data * const data, static_storage * const workspace,
+ output_controls * const out_control, int first_run )
{
FILE *ffield, *ctrl;
ffield = sfopen( ffield_file, "r" );
ctrl = sfopen( control_file, "r" );
- Read_Force_Field( ffield, &system->reax_param );
+ Read_Force_Field( ffield, &system->reax_param, first_run );
Read_Control_File( ctrl, system, control, out_control );
if ( control->geo_format == CUSTOM )
{
- Read_Geo( geo_file, system, control, data, workspace );
+ Read_Geo( geo_file, system, control, data, workspace, first_run );
}
else if ( control->geo_format == PDB )
{
- Read_PDB( geo_file, system, control, data, workspace );
+ Read_PDB( geo_file, system, control, data, workspace, first_run );
}
else if ( control->geo_format == BGF )
{
- Read_BGF( geo_file, system, control, data, workspace );
+ Read_BGF( geo_file, system, control, data, workspace, first_run );
}
else if ( control->geo_format == ASCII_RESTART )
{
- Read_ASCII_Restart( geo_file, system, control, data, workspace );
+ Read_ASCII_Restart( geo_file, system, control, data, workspace, first_run );
control->restart = TRUE;
}
else if ( control->geo_format == BINARY_RESTART )
{
- Read_Binary_Restart( geo_file, system, control, data, workspace );
+ Read_Binary_Restart( geo_file, system, control, data, workspace, first_run );
control->restart = TRUE;
}
else
@@ -129,16 +135,22 @@ static void Read_System( const char * const geo_file,
exit( INVALID_GEO );
}
- sfclose( ffield, "Read_System::ffield" );
- sfclose( ctrl, "Read_System::ctrl" );
+ sfclose( ffield, "Read_Input_Files::ffield" );
+ sfclose( ctrl, "Read_Input_Files::ctrl" );
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "input files have been read...\n" );
Print_Box( &system->box, stderr );
#endif
}
+/* Allocate top-level data structures and parse input files
+ * for the first simulation
+ *
+ * geo_file: file containing geometry info of the structure to simulate
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
void* setup( const char * const geo_file, const char * const ffield_file,
const char * const control_file )
{
@@ -169,19 +181,28 @@ void* setup( const char * const geo_file, const char * const ffield_file,
spmd_handle->out_control = smalloc( sizeof(output_controls),
"Setup::spmd_handle->out_control" );
+ spmd_handle->first_run = TRUE;
spmd_handle->output_enabled = TRUE;
+ spmd_handle->realloc = TRUE;
spmd_handle->callback = NULL;
+ spmd_handle->data->sim_id = 0;
- /* parse geometry file */
- Read_System( geo_file, ffield_file, control_file,
+ Read_Input_Files( geo_file, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->workspace,
- spmd_handle->out_control );
+ spmd_handle->out_control, spmd_handle->first_run );
+
+ spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
return (void*) spmd_handle;
}
+/* Setup callback function to be run after each simulation step
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * callback: function pointer to attach for callback
+ */
int setup_callback( const void * const handle, const callback_function callback )
{
int ret;
@@ -201,6 +222,10 @@ int setup_callback( const void * const handle, const callback_function callback
}
+/* Run the simulation according to the prescribed parameters
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ */
int simulate( const void * const handle )
{
int steps, ret;
@@ -216,7 +241,10 @@ int simulate( const void * const handle )
Initialize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists,
spmd_handle->out_control, &Evolve,
- spmd_handle->output_enabled );
+ spmd_handle->output_enabled,
+ spmd_handle->first_run, spmd_handle->realloc );
+
+ spmd_handle->realloc = FALSE;
/* compute f_0 */
//if( control.restart == FALSE ) {
@@ -248,7 +276,10 @@ int simulate( const void * const handle )
{
Compute_Total_Energy( spmd_handle->data );
}
+ }
+ if ( spmd_handle->output_enabled == TRUE )
+ {
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
}
@@ -338,6 +369,10 @@ int simulate( const void * const handle )
}
+/* Deallocate all data structures post-simulation
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ */
int cleanup( const void * const handle )
{
int i, ret;
@@ -351,7 +386,7 @@ int cleanup( const void * const handle )
Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
- spmd_handle->output_enabled );
+ spmd_handle->output_enabled, FALSE );
sfree( spmd_handle->out_control, "cleanup::spmd_handle->out_control" );
for ( i = 0; i < LIST_N; ++i )
@@ -373,6 +408,59 @@ int cleanup( const void * const handle )
}
+/* Reset for the next simulation by parsing input files and triggering
+ * reallocation if more space is needed
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * geo_file: file containing geometry info of the structure to simulate
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
+int reset( const void * const handle, const char * const geo_file,
+ const char * const ffield_file, const char * const control_file )
+{
+ int ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ spmd_handle->first_run = FALSE;
+ spmd_handle->realloc = FALSE;
+ spmd_handle->data->sim_id++;
+
+ Read_Input_Files( geo_file, ffield_file, control_file,
+ spmd_handle->system, spmd_handle->control,
+ spmd_handle->data, spmd_handle->workspace,
+ spmd_handle->out_control, spmd_handle->first_run );
+
+ if ( spmd_handle->system->N > spmd_handle->system->N_max )
+ {
+ /* deallocate everything which needs more space
+ * (i.e., structures whose space is a function of the number of atoms),
+ * except for data structures allocated while parsing input files */
+ Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
+ spmd_handle->output_enabled, TRUE );
+
+ spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
+ spmd_handle->realloc = TRUE;
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+/* Getter for atom info
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ */
reax_atom* get_atoms( const void * const handle )
{
spuremd_handle *spmd_handle;
@@ -390,6 +478,11 @@ reax_atom* get_atoms( const void * const handle )
}
+/* Setter for writing output to files
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * enabled: TRUE enables writing output to files, FALSE otherwise
+ */
int set_output_enabled( const void * const handle, const int enabled )
{
int ret;
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index 64e263d1475e2d01ab571432f58db630cff974aa..980d3448981a8a9ed6c05e23391a46c84f046fc8 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -42,6 +42,9 @@ int simulate( const void * const );
int cleanup( const void * const );
+int reset( const void * const, const char * const,
+ const char * const, const char * const );
+
reax_atom* get_atoms( const void * const );
int set_output_enabled( const void * const, const int );
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index 5a9acd17641abc02825d6c9c0874cb3a3eeb0960..05f820003cea38039d2375b91a4e220e04c1ae4a 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -74,8 +74,7 @@ void Compute_Total_Mass( reax_system *system, simulation_data *data )
}
-void Compute_Center_of_Mass( reax_system *system, simulation_data *data,
- FILE *fout )
+void Compute_Center_of_Mass( reax_system *system, simulation_data *data )
{
int i;
real m, xx, xy, xz, yy, yz, zz, det;
diff --git a/sPuReMD/src/system_props.h b/sPuReMD/src/system_props.h
index 3b8d6a98f25188dd5a1dd604a2a90e15a5622880..26e1594b8cce64b7ee85db5cc3dee5c40f3442e6 100644
--- a/sPuReMD/src/system_props.h
+++ b/sPuReMD/src/system_props.h
@@ -29,7 +29,7 @@ void Temperature_Control( control_params*, simulation_data*, output_controls* );
void Compute_Total_Mass( reax_system*, simulation_data* );
-void Compute_Center_of_Mass( reax_system*, simulation_data*, FILE* );
+void Compute_Center_of_Mass( reax_system*, simulation_data* );
void Compute_Kinetic_Energy( reax_system*, simulation_data* );