From 78362ca3a07ba7cfd5233340f04b65e349ea00a5 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Tue, 15 Dec 2020 00:46:07 -0500
Subject: [PATCH] sPuReMD: finalize sPuReMD multiple consecutive simulation
code (for integration with other codebases).
---
sPuReMD/src/charges.c | 366 ++++++++-------------------------------
sPuReMD/src/charges.h | 7 +-
sPuReMD/src/ffield.c | 66 +++++--
sPuReMD/src/forces.c | 8 +-
sPuReMD/src/forces.h | 8 +-
sPuReMD/src/init_md.c | 50 +++---
sPuReMD/src/init_md.h | 3 +
sPuReMD/src/integrate.c | 25 ++-
sPuReMD/src/integrate.h | 26 +--
sPuReMD/src/lin_alg.c | 43 ++---
sPuReMD/src/lin_alg.h | 4 +-
sPuReMD/src/list.c | 22 +++
sPuReMD/src/reax_types.h | 21 ++-
sPuReMD/src/spuremd.c | 72 ++++++--
sPuReMD/src/spuremd.h | 5 +-
tools/driver.py | 29 +++-
16 files changed, 338 insertions(+), 417 deletions(-)
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 57bc2e0c..c9fde2b3 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -491,8 +491,8 @@ static void Spline_Extrapolate_Charges_EE( const reax_system * const system,
/* Compute preconditioner for QEq
*/
static void Compute_Preconditioner_QEq( const reax_system * const system,
- const control_params * const control,
- simulation_data * const data, static_storage * const workspace )
+ const control_params * const control, simulation_data * const data,
+ static_storage * const workspace, int realloc )
{
real time;
sparse_matrix *Hptr;
@@ -513,7 +513,8 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
time = Get_Time( );
if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
- setup_graph_coloring( control, workspace, Hptr, &workspace->H_full, &workspace->H_p );
+ setup_graph_coloring( control, workspace, Hptr, &workspace->H_full,
+ &workspace->H_p, realloc );
Sort_Matrix_Rows( workspace->H_p );
Hptr = workspace->H_p;
}
@@ -612,259 +613,11 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
}
-/* Compute preconditioner for EE
- */
-//static void Compute_Preconditioner_EE( const reax_system * const system,
-// const control_params * const control,
-// simulation_data * const data, static_storage * const workspace )
-//{
-// int i, top;
-// static real * ones = NULL, * x = NULL, * y = NULL;
-// sparse_matrix *Hptr;
-//
-// Hptr = workspace->H_EE;
-//
-//#if defined(TEST_MAT)
-// Hptr = create_test_mat( );
-//#endif
-//
-// if ( ones == NULL )
-// {
-// ones = smalloc( system->N * sizeof(real), "Compute_Preconditioner_EE::ones" );
-// x = smalloc( system->N * sizeof(real), "Compute_Preconditioner_EE::x" );
-// y = smalloc( system->N * sizeof(real), "Compute_Preconditioner_EE::y" );
-//
-// for ( i = 0; i < system->N; ++i )
-// {
-// ones[i] = 1.0;
-// }
-// }
-//
-// switch ( control->cm_solver_pre_comp_type )
-// {
-// case JACOBI_PC:
-// data->timing.cm_solver_pre_comp += jacobi( Hptr, workspace->Hdia_inv );
-// break;
-//
-// case ICHOLT_PC:
-// fprintf( stderr, "[ERROR] ICHOLT is not supported for indefinite, symmetric matrices of EE. Terminating...\n" );
-// exit( INVALID_INPUT );
-// break;
-//
-// case ILUT_PC:
-// if ( control->cm_solver_pre_comp_droptol > 0.0 )
-// {
-// data->timing.cm_solver_pre_comp +=
-// ILUT( Hptr, workspace->droptol, workspace->L, workspace->U );
-// }
-// else
-// {
-// data->timing.cm_solver_pre_comp +=
-// ILU( Hptr, workspace->L, workspace->U );
-// }
-// break;
-//
-// case ILUTP_PC:
-// data->timing.cm_solver_pre_comp +=
-// ILUTP( Hptr, workspace->droptol, workspace->L, workspace->U );
-// break;
-//
-// case FG_ILUT_PC:
-// data->timing.cm_solver_pre_comp +=
-// FG_ILUT( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
-// workspace->L_EE, workspace->U_EE );
-// break;
-//
-// default:
-// fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method (%d). Terminating...\n",
-// control->cm_solver_pre_comp_type );
-// exit( INVALID_INPUT );
-// break;
-// }
-//
-// if ( control->cm_solver_pre_comp_type != JACOBI_PC )
-// {
-// switch ( control->cm_solver_pre_app_type )
-// {
-// case TRI_SOLVE_PA:
-// tri_solve( workspace->L_EE, ones, x, LOWER );
-// Transpose_I( workspace->U_EE );
-// tri_solve( workspace->U_EE, ones, y, LOWER );
-// Transpose_I( workspace->U_EE );
-//
-// memcpy( workspace->L->start, workspace->L_EE->start, sizeof(unsigned int) * (system->N + 1) );
-// memcpy( workspace->L->j, workspace->L_EE->j, sizeof(unsigned int) * workspace->L_EE->start[workspace->L_EE->n] );
-// memcpy( workspace->L->val, workspace->L_EE->val, sizeof(real) * workspace->L_EE->start[workspace->L_EE->n] );
-//
-// top = workspace->L->start[system->N];
-// for ( i = 0; i < system->N; ++i )
-// {
-// workspace->L->j[top] = i;
-// workspace->L->val[top] = x[i];
-// ++top;
-// }
-//
-// workspace->L->j[top] = system->N_cm - 1;
-// workspace->L->val[top] = 1.0;
-// ++top;
-//
-// workspace->L->start[system->N_cm] = top;
-//
-// top = 0;
-// for ( i = 0; i < system->N; ++i )
-// {
-// workspace->U->start[i] = top;
-// memcpy( workspace->U->j + top, workspace->U_EE->j + workspace->U_EE->start[i],
-// sizeof(unsigned int) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
-// memcpy( workspace->U->val + top, workspace->U_EE->val + workspace->U_EE->start[i],
-// sizeof(real) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
-// top += (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]);
-//
-// workspace->U->j[top] = system->N_cm - 1;
-// workspace->U->val[top] = y[i];
-// ++top;
-// }
-//
-// workspace->U->start[system->N_cm - 1] = top;
-//
-// workspace->U->j[top] = system->N_cm - 1;
-// workspace->U->val[top] = -Dot( x, y, system->N );
-// ++top;
-//
-// workspace->U->start[system->N_cm] = top;
-// break;
-//
-// case TRI_SOLVE_LEVEL_SCHED_PA:
-// tri_solve_level_sched( workspace->L_EE, ones, x, LOWER, TRUE );
-// Transpose_I( workspace->U_EE );
-// tri_solve_level_sched( workspace->U_EE, ones, y, LOWER, TRUE );
-// Transpose_I( workspace->U_EE );
-//
-// memcpy( workspace->L->start, workspace->L_EE->start, sizeof(unsigned int) * (system->N + 1) );
-// memcpy( workspace->L->j, workspace->L_EE->j, sizeof(unsigned int) * workspace->L_EE->start[workspace->L_EE->n] );
-// memcpy( workspace->L->val, workspace->L_EE->val, sizeof(real) * workspace->L_EE->start[workspace->L_EE->n] );
-//
-// top = workspace->L->start[system->N];
-// for ( i = 0; i < system->N; ++i )
-// {
-// workspace->L->j[top] = i;
-// workspace->L->val[top] = x[i];
-// ++top;
-// }
-//
-// workspace->L->j[top] = system->N_cm - 1;
-// workspace->L->val[top] = 1.0;
-// ++top;
-//
-// workspace->L->start[system->N_cm] = top;
-//
-// top = 0;
-// for ( i = 0; i < system->N; ++i )
-// {
-// workspace->U->start[i] = top;
-// memcpy( workspace->U->j + top, workspace->U_EE->j + workspace->U_EE->start[i],
-// sizeof(unsigned int) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
-// memcpy( workspace->U->val + top, workspace->U_EE->val + workspace->U_EE->start[i],
-// sizeof(real) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
-// top += (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]);
-//
-// workspace->U->j[top] = system->N_cm - 1;
-// workspace->U->val[top] = y[i];
-// ++top;
-// }
-//
-// workspace->U->start[system->N_cm - 1] = top;
-//
-// workspace->U->j[top] = system->N_cm - 1;
-// workspace->U->val[top] = -Dot( x, y, system->N );
-// ++top;
-//
-// workspace->U->start[system->N_cm] = top;
-// break;
-//
-// //TODO: add Jacobi iter, etc.?
-// default:
-// tri_solve( workspace->L_EE, ones, x, LOWER );
-// Transpose_I( workspace->U_EE );
-// tri_solve( workspace->U_EE, ones, y, LOWER );
-// Transpose_I( workspace->U_EE );
-//
-// memcpy( workspace->L->start, workspace->L_EE->start, sizeof(unsigned int) * (system->N + 1) );
-// memcpy( workspace->L->j, workspace->L_EE->j, sizeof(unsigned int) * workspace->L_EE->start[workspace->L_EE->n] );
-// memcpy( workspace->L->val, workspace->L_EE->val, sizeof(real) * workspace->L_EE->start[workspace->L_EE->n] );
-//
-// top = workspace->L->start[system->N];
-// for ( i = 0; i < system->N; ++i )
-// {
-// workspace->L->j[top] = i;
-// workspace->L->val[top] = x[i];
-// ++top;
-// }
-//
-// workspace->L->j[top] = system->N_cm - 1;
-// workspace->L->val[top] = 1.0;
-// ++top;
-//
-// workspace->L->start[system->N_cm] = top;
-//
-// top = 0;
-// for ( i = 0; i < system->N; ++i )
-// {
-// workspace->U->start[i] = top;
-// memcpy( workspace->U->j + top, workspace->U_EE->j + workspace->U_EE->start[i],
-// sizeof(unsigned int) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
-// memcpy( workspace->U->val + top, workspace->U_EE->val + workspace->U_EE->start[i],
-// sizeof(real) * (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]) );
-// top += (workspace->U_EE->start[i + 1] - workspace->U_EE->start[i]);
-//
-// workspace->U->j[top] = system->N_cm - 1;
-// workspace->U->val[top] = y[i];
-// ++top;
-// }
-//
-// workspace->U->start[system->N_cm - 1] = top;
-//
-// workspace->U->j[top] = system->N_cm - 1;
-// workspace->U->val[top] = -Dot( x, y, system->N );
-// ++top;
-//
-// workspace->U->start[system->N_cm] = top;
-// break;
-// }
-// }
-//
-//#if defined(DEBUG)
-// if ( control->cm_solver_pre_comp_type == ICHOLT_PC
-// || control->cm_solver_pre_comp_type == ILUT_PC
-// || control->cm_solver_pre_comp_type == ILUTP_PC
-// || control->cm_solver_pre_comp_type == FG_ILUT_PC )
-// {
-// fprintf( stderr, "[INFO] condest = %f\n", condest(workspace->L) );
-//
-//#if defined(DEBUG_FOCUS)
-//#define SIZE (1000)
-// char fname[SIZE];
-// snprintf( fname, SIZE, "%s.L%d.out", control->sim_name, data->step );
-// Print_Sparse_Matrix2( workspace->L, fname, NULL );
-// snprintf( fname, SIZE, "%s.U%d.out", control->sim_name, data->step );
-// Print_Sparse_Matrix2( workspace->U, fname, NULL );
-//
-// fprintf( stderr, "icholt-" );
-// snprintf( fname, SIZE, "%s.L%d.out", control->sim_name, data->step );
-// Print_Sparse_Matrix2( workspace->L, fname, NULL );
-// Print_Sparse_Matrix( U );
-//#undef SIZE
-//#endif
-// }
-//#endif
-//}
-
-
/* Compute preconditioner for EE
*/
static void Compute_Preconditioner_EE( const reax_system * const system,
- const control_params * const control,
- simulation_data * const data, static_storage * const workspace )
+ const control_params * const control, simulation_data * const data,
+ static_storage * const workspace, int realloc )
{
real time;
sparse_matrix *Hptr;
@@ -891,7 +644,8 @@ static void Compute_Preconditioner_EE( const reax_system * const system,
time = Get_Time( );
if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
- setup_graph_coloring( control, workspace, Hptr, &workspace->H_full, &workspace->H_p );
+ setup_graph_coloring( control, workspace, Hptr, &workspace->H_full,
+ &workspace->H_p, realloc );
Sort_Matrix_Rows( workspace->H_p );
Hptr = workspace->H_p;
}
@@ -1002,8 +756,8 @@ static void Compute_Preconditioner_EE( const reax_system * const system,
/* Compute preconditioner for ACKS2
*/
static void Compute_Preconditioner_ACKS2( const reax_system * const system,
- const control_params * const control,
- simulation_data * const data, static_storage * const workspace )
+ const control_params * const control, simulation_data * const data,
+ static_storage * const workspace, int realloc )
{
real time;
sparse_matrix *Hptr;
@@ -1031,7 +785,8 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
time = Get_Time( );
if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
- setup_graph_coloring( control, workspace, Hptr, &workspace->H_full, &workspace->H_p );
+ setup_graph_coloring( control, workspace, Hptr, &workspace->H_full,
+ &workspace->H_p, realloc );
Sort_Matrix_Rows( workspace->H_p );
Hptr = workspace->H_p;
}
@@ -1141,8 +896,8 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
static void Setup_Preconditioner_QEq( const reax_system * const system,
- const control_params * const control,
- simulation_data * const data, static_storage * const workspace )
+ const control_params * const control, simulation_data * const data,
+ static_storage * const workspace, int realloc )
{
int fillin;
real time;
@@ -1174,7 +929,13 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
case JACOBI_PC:
if ( workspace->Hdia_inv == NULL )
{
- workspace->Hdia_inv = scalloc( Hptr->n, sizeof( real ),
+ workspace->Hdia_inv = scalloc( Hptr->n_max, sizeof( real ),
+ "Setup_Preconditioner_QEq::workspace->Hdia_inv" );
+ }
+ else if ( realloc == TRUE )
+ {
+ workspace->Hdia_inv = srealloc( workspace->Hdia_inv,
+ sizeof( real ) * Hptr->n_max,
"Setup_Preconditioner_QEq::workspace->Hdia_inv" );
}
break;
@@ -1189,7 +950,7 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
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 )
+ else if ( workspace->L->m < fillin || realloc == TRUE )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
@@ -1210,7 +971,7 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
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 )
+ else if ( workspace->L->m < Hptr->m || realloc == TRUE )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
@@ -1231,7 +992,7 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
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 )
+ else if ( workspace->L->m < Hptr->m || realloc == TRUE )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
@@ -1244,9 +1005,10 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
break;
case SAI_PC:
- setup_sparse_approx_inverse( Hptr, &workspace->H_full, &workspace->H_spar_patt,
- &workspace->H_spar_patt_full, &workspace->H_app_inv,
- control->cm_solver_pre_comp_sai_thres );
+ setup_sparse_approx_inverse( Hptr, &workspace->H_full,
+ &workspace->H_spar_patt, &workspace->H_spar_patt_full,
+ &workspace->H_app_inv, control->cm_solver_pre_comp_sai_thres,
+ realloc );
break;
default:
@@ -1261,8 +1023,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
/* Setup routines before computing the preconditioner for EE
*/
static void Setup_Preconditioner_EE( const reax_system * const system,
- const control_params * const control,
- simulation_data * const data, static_storage * const workspace )
+ const control_params * const control, simulation_data * const data,
+ static_storage * const workspace, int realloc )
{
real time;
sparse_matrix *Hptr;
@@ -1293,7 +1055,13 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
case JACOBI_PC:
if ( workspace->Hdia_inv == NULL )
{
- workspace->Hdia_inv = scalloc( Hptr->n, sizeof( real ),
+ workspace->Hdia_inv = scalloc( Hptr->n_max, sizeof( real ),
+ "Setup_Preconditioner_EE::workspace->Hdiv_inv" );
+ }
+ else if ( realloc == TRUE )
+ {
+ workspace->Hdia_inv = srealloc( workspace->Hdia_inv,
+ sizeof( real ) * Hptr->n_max,
"Setup_Preconditioner_EE::workspace->Hdiv_inv" );
}
break;
@@ -1320,7 +1088,7 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
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 )
+ else if ( workspace->L->m < Hptr->m || realloc == TRUE )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
@@ -1347,7 +1115,7 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
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 )
+ else if ( workspace->L->m < Hptr->m || realloc == TRUE )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
@@ -1360,9 +1128,10 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
break;
case SAI_PC:
- setup_sparse_approx_inverse( Hptr, &workspace->H_full, &workspace->H_spar_patt,
- &workspace->H_spar_patt_full, &workspace->H_app_inv,
- control->cm_solver_pre_comp_sai_thres );
+ setup_sparse_approx_inverse( Hptr, &workspace->H_full,
+ &workspace->H_spar_patt, &workspace->H_spar_patt_full,
+ &workspace->H_app_inv, control->cm_solver_pre_comp_sai_thres,
+ realloc );
break;
default:
@@ -1377,8 +1146,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
/* Setup routines before computing the preconditioner for ACKS2
*/
static void Setup_Preconditioner_ACKS2( const reax_system * const system,
- const control_params * const control,
- simulation_data * const data, static_storage * const workspace )
+ const control_params * const control, simulation_data * const data,
+ static_storage * const workspace, int realloc )
{
real time;
sparse_matrix *Hptr;
@@ -1409,7 +1178,13 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
case JACOBI_PC:
if ( workspace->Hdia_inv == NULL )
{
- workspace->Hdia_inv = scalloc( Hptr->n, sizeof( real ),
+ workspace->Hdia_inv = scalloc( Hptr->n_max, sizeof( real ),
+ "Setup_Preconditioner_ACKS2::workspace->Hdiv_inv" );
+ }
+ else if ( realloc == TRUE )
+ {
+ workspace->Hdia_inv = srealloc( workspace->Hdia_inv,
+ sizeof( real ) * Hptr->n_max,
"Setup_Preconditioner_ACKS2::workspace->Hdiv_inv" );
}
break;
@@ -1438,7 +1213,7 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
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 )
+ else if ( workspace->L->m < Hptr->m || realloc == TRUE )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
@@ -1467,7 +1242,7 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
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 )
+ else if ( workspace->L->m < Hptr->m || realloc == TRUE )
{
Deallocate_Matrix( workspace->L );
Deallocate_Matrix( workspace->U );
@@ -1479,9 +1254,10 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
break;
case SAI_PC:
- setup_sparse_approx_inverse( Hptr, &workspace->H_full, &workspace->H_spar_patt,
- &workspace->H_spar_patt_full, &workspace->H_app_inv,
- control->cm_solver_pre_comp_sai_thres );
+ setup_sparse_approx_inverse( Hptr, &workspace->H_full,
+ &workspace->H_spar_patt, &workspace->H_spar_patt_full,
+ &workspace->H_app_inv, control->cm_solver_pre_comp_sai_thres,
+ realloc );
break;
default:
@@ -1572,7 +1348,7 @@ static void Calculate_Charges_ACKS2( const reax_system * const system,
*/
static void QEq( reax_system * const system, control_params * const control,
simulation_data * const data, static_storage * const workspace,
- const output_controls * const out_control )
+ const output_controls * const out_control, int realloc )
{
int iters, refactor;
@@ -1580,9 +1356,9 @@ static void QEq( reax_system * const system, control_params * const control,
if ( refactor == TRUE )
{
- Setup_Preconditioner_QEq( system, control, data, workspace );
+ Setup_Preconditioner_QEq( system, control, data, workspace, realloc );
- Compute_Preconditioner_QEq( system, control, data, workspace );
+ Compute_Preconditioner_QEq( system, control, data, workspace, realloc );
}
switch ( control->cm_init_guess_type )
@@ -1684,7 +1460,7 @@ static void QEq( reax_system * const system, control_params * const control,
*/
static void EE( reax_system * const system, control_params * const control,
simulation_data * const data, static_storage * const workspace,
- const output_controls * const out_control )
+ const output_controls * const out_control, int realloc )
{
int iters, refactor;
@@ -1692,9 +1468,9 @@ static void EE( reax_system * const system, control_params * const control,
if ( refactor == TRUE )
{
- Setup_Preconditioner_EE( system, control, data, workspace );
+ Setup_Preconditioner_EE( system, control, data, workspace, realloc );
- Compute_Preconditioner_EE( system, control, data, workspace );
+ Compute_Preconditioner_EE( system, control, data, workspace, realloc );
}
switch ( control->cm_init_guess_type )
@@ -1778,7 +1554,7 @@ static void EE( reax_system * const system, control_params * const control,
*/
static void ACKS2( reax_system * const system, control_params * const control,
simulation_data * const data, static_storage * const workspace,
- const output_controls * const out_control )
+ const output_controls * const out_control, int realloc )
{
int iters, refactor;
@@ -1786,9 +1562,9 @@ static void ACKS2( reax_system * const system, control_params * const control,
if ( refactor == TRUE )
{
- Setup_Preconditioner_ACKS2( system, control, data, workspace );
+ Setup_Preconditioner_ACKS2( system, control, data, workspace, realloc );
- Compute_Preconditioner_ACKS2( system, control, data, workspace );
+ Compute_Preconditioner_ACKS2( system, control, data, workspace, realloc );
}
// Print_Linear_System( system, control, workspace, data->step );
@@ -1879,7 +1655,7 @@ static void ACKS2( reax_system * const system, control_params * const control,
void Compute_Charges( reax_system * const system, control_params * const control,
simulation_data * const data, static_storage * const workspace,
- const output_controls * const out_control )
+ const output_controls * const out_control, int realloc )
{
#if defined(DEBUG_FOCUS)
#define SIZE (200)
@@ -1908,15 +1684,15 @@ void Compute_Charges( reax_system * const system, control_params * const control
switch ( control->charge_method )
{
case QEQ_CM:
- QEq( system, control, data, workspace, out_control );
+ QEq( system, control, data, workspace, out_control, realloc );
break;
case EE_CM:
- EE( system, control, data, workspace, out_control );
+ EE( system, control, data, workspace, out_control, realloc );
break;
case ACKS2_CM:
- ACKS2( system, control, data, workspace, out_control );
+ ACKS2( system, control, data, workspace, out_control, realloc );
break;
default:
diff --git a/sPuReMD/src/charges.h b/sPuReMD/src/charges.h
index 3c438331..b242fffe 100644
--- a/sPuReMD/src/charges.h
+++ b/sPuReMD/src/charges.h
@@ -25,9 +25,10 @@
#include "reax_types.h"
-void Compute_Charges( reax_system* const, control_params* const, simulation_data* const,
- static_storage* const, const output_controls* const );
+void Compute_Charges( reax_system * const, control_params * const,
+ simulation_data * const, static_storage * const,
+ const output_controls * const, int );
-int is_refactoring_step ( control_params* const, simulation_data* const );
+int is_refactoring_step( control_params * const, simulation_data * const );
#endif
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index f536dc36..13ab1dba 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -40,7 +40,7 @@ void Read_Force_Field( FILE *fp, reax_interaction *reax, int first_run )
{
s = smalloc( sizeof(char) * MAX_LINE, "Read_Force_Field::s" );
tmp = smalloc( sizeof(char*) * MAX_TOKENS, "Read_Force_Field::tmp" );
- for (i = 0; i < MAX_TOKENS; i++)
+ for ( i = 0; i < MAX_TOKENS; i++ )
{
tmp[i] = smalloc( sizeof(char) * MAX_TOKEN_LEN, "Read_Force_Field::tmp[i]" );
}
@@ -64,11 +64,15 @@ void Read_Force_Field( FILE *fp, reax_interaction *reax, int first_run )
{
reax->gp.l = (real*) smalloc( sizeof(real) * n,
"Read_Force_Field::reax->gp-l" );
+
+ reax->gp.max_n_global = n;
}
- else if ( reax->gp.n_global < n )
+ else if ( reax->gp.max_n_global < n )
{
reax->gp.l = (real*) srealloc( reax->gp.l, sizeof(real) * n,
"Read_Force_Field::reax->gp-l" );
+
+ reax->gp.max_n_global = n;
}
reax->gp.n_global = n;
@@ -134,48 +138,78 @@ void Read_Force_Field( FILE *fp, reax_interaction *reax, int first_run )
}
}
}
+
+ reax->max_num_atom_types = n;
}
- else if ( reax->num_atom_types < n )
+ else if ( reax->max_num_atom_types < n )
{
+ for ( i = 0; i < reax->max_num_atom_types; i++ )
+ {
+ for ( j = 0; j < reax->max_num_atom_types; j++ )
+ {
+ for ( k = 0; k < reax->max_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->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" );
+
/* Allocating structures in reax_interaction */
- reax->sbp = (single_body_parameters*) srealloc( reax->sbp, n * sizeof(single_body_parameters),
+ reax->sbp = (single_body_parameters*) scalloc( n, sizeof(single_body_parameters),
"Read_Force_Field::reax->sbp" );
- reax->tbp = (two_body_parameters**) srealloc( reax->tbp, n * sizeof(two_body_parameters*),
+ reax->tbp = (two_body_parameters**) scalloc( n, sizeof(two_body_parameters*),
"Read_Force_Field::reax->tbp" );
- reax->thbp = (three_body_header***) srealloc( reax->thbp, n * sizeof(three_body_header**),
+ reax->thbp = (three_body_header***) scalloc( n, sizeof(three_body_header**),
"Read_Force_Field::reax->thbp" );
- reax->hbp = (hbond_parameters***) srealloc( reax->hbp, n * sizeof(hbond_parameters**),
+ reax->hbp = (hbond_parameters***) scalloc( n, sizeof(hbond_parameters**),
"Read_Force_Field::reax->hbp" );
- reax->fbp = (four_body_header****) srealloc( reax->fbp, n * sizeof(four_body_header***),
+ 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*) srealloc( reax->tbp[i], n * sizeof(two_body_parameters),
+ reax->tbp[i] = (two_body_parameters*) scalloc( 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*),
+ reax->thbp[i] = (three_body_header**) scalloc( n, sizeof(three_body_header*),
"Read_Force_Field::reax->thbp[i]" );
- reax->hbp[i] = (hbond_parameters**) srealloc( reax->hbp[i], n * sizeof(hbond_parameters*),
+ reax->hbp[i] = (hbond_parameters**) scalloc( 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**),
+ 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*) srealloc( reax->thbp[i][j], n * sizeof(three_body_header),
+ 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*) srealloc( reax->hbp[i][j], n * sizeof(hbond_parameters),
+ 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**) srealloc( reax->fbp[i][j], n * sizeof(four_body_header*),
+ 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*) srealloc( reax->fbp[i][j][k], n * sizeof(four_body_header),
+ reax->fbp[i][j][k] = (four_body_header*) scalloc( n, sizeof(four_body_header),
"Read_Force_Field::reax->fbp[i][j][k]" );
}
}
}
+
+ reax->max_num_atom_types = n;
}
tor_flag = (char****) smalloc( n * sizeof(char***),
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index e08760b1..2608c65b 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -138,7 +138,7 @@ static void Compute_Bonded_Forces( reax_system *system, control_params *control,
static void Compute_NonBonded_Forces( 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 realloc )
{
real t_start, t_elapsed;
@@ -150,7 +150,7 @@ static void Compute_NonBonded_Forces( reax_system *system, control_params *contr
#endif
t_start = Get_Time( );
- Compute_Charges( system, control, data, workspace, out_control );
+ Compute_Charges( system, control, data, workspace, out_control, realloc );
t_elapsed = Get_Timing_Info( t_start );
data->timing.cm += t_elapsed;
@@ -1513,7 +1513,7 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
void Compute_Forces( 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 realloc )
{
real t_start, t_elapsed;
@@ -1536,7 +1536,7 @@ void Compute_Forces( reax_system *system, control_params *control,
t_start = Get_Time( );
Compute_NonBonded_Forces( system, control, data, workspace,
- lists, out_control );
+ lists, out_control, realloc );
t_elapsed = Get_Timing_Info( t_start );
data->timing.nonb += t_elapsed;
diff --git a/sPuReMD/src/forces.h b/sPuReMD/src/forces.h
index 78587a50..982ba286 100644
--- a/sPuReMD/src/forces.h
+++ b/sPuReMD/src/forces.h
@@ -27,11 +27,11 @@
void Init_Bonded_Force_Functions( control_params* );
-void Compute_Forces( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+void Compute_Forces( reax_system *, control_params *, simulation_data *,
+ static_storage *, reax_list **, output_controls *, int );
-void Estimate_Storage_Sizes( reax_system*, control_params*, reax_list**,
- int*, int*, int*, int* );
+void Estimate_Storage_Sizes( reax_system *, control_params *, reax_list **,
+ int *, int *, int *, int * );
#endif
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 9809a3bf..d86ceb48 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -160,10 +160,10 @@ static void Init_System( reax_system *system, control_params *control,
static void Init_Simulation_Data( reax_system *system, control_params *control,
simulation_data *data, output_controls *out_control,
- evolve_function *Evolve, int alloc )
+ evolve_function *Evolve, int realloc )
{
#if defined(_OPENMP)
- if ( alloc == TRUE && (control->ensemble == sNPT || control->ensemble == iNPT
+ if ( realloc == TRUE && (control->ensemble == sNPT || control->ensemble == iNPT
|| control->ensemble == aNPT || control->compute_pressure == TRUE) )
{
data->press_local = smalloc( sizeof( rtensor ) * control->num_threads,
@@ -323,11 +323,11 @@ static void Init_Taper( control_params *control, static_storage *workspace )
static void Init_Workspace( reax_system *system, control_params *control,
- static_storage *workspace, int alloc )
+ static_storage *workspace, int realloc )
{
int i;
- if ( alloc == TRUE )
+ if ( realloc == TRUE )
{
/* hydrogen bond list */
workspace->hbond_index = smalloc( system->N_max * sizeof( int ),
@@ -390,7 +390,7 @@ static void Init_Workspace( reax_system *system, control_params *control,
break;
}
- if ( alloc == TRUE )
+ if ( realloc == TRUE )
{
/* sparse matrices */
workspace->H = NULL;
@@ -478,7 +478,7 @@ static void Init_Workspace( reax_system *system, control_params *control,
break;
}
- if ( alloc == TRUE )
+ if ( realloc == TRUE )
{
switch ( control->cm_solver_type )
{
@@ -578,7 +578,7 @@ static void Init_Workspace( reax_system *system, control_params *control,
/* level scheduling related */
workspace->levels_L = 1;
workspace->levels_U = 1;
- if ( alloc == TRUE )
+ if ( realloc == TRUE )
{
if ( control->cm_solver_pre_app_type == TRI_SOLVE_LEVEL_SCHED_PA ||
control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
@@ -612,7 +612,7 @@ static void Init_Workspace( reax_system *system, control_params *control,
/* graph coloring related */
workspace->recolor_cnt = 0;
- if ( alloc == TRUE )
+ if ( realloc == TRUE )
{
if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
@@ -775,7 +775,7 @@ 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,
- int first_run, int alloc )
+ int first_run, int realloc )
{
int i, num_nbrs, num_bonds, num_hbonds, num_3body, Htop, max_nnz;
int *hb_top, *bond_top;
@@ -786,7 +786,7 @@ static void Init_Lists( reax_system *system, control_params *control,
{
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 )
+ else if ( realloc == 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] );
@@ -827,7 +827,7 @@ static void Init_Lists( reax_system *system, control_params *control,
{
Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max, max_nnz );
}
- else if ( alloc == TRUE || workspace->H->m < max_nnz )
+ else if ( realloc == TRUE || workspace->H->m < max_nnz )
{
Deallocate_Matrix( workspace->H );
Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max, max_nnz );
@@ -845,7 +845,7 @@ static void Init_Lists( reax_system *system, control_params *control,
* (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 )
+ else if ( realloc == TRUE || workspace->H_sp->m < max_nnz )
{
Deallocate_Matrix( workspace->H_sp );
/* TODO: better estimate for H_sp?
@@ -887,7 +887,7 @@ static void Init_Lists( reax_system *system, control_params *control,
num_hbonds += hb_top[i];
}
- if ( first_run == TRUE )
+ if ( first_run == TRUE || lists[HBONDS]->allocated == FALSE )
{
workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
@@ -926,11 +926,11 @@ static void Init_Lists( reax_system *system, control_params *control,
}
/* bonds list */
- if ( first_run == TRUE )
+ if ( lists[BONDS]->allocated == FALSE )
{
Allocate_Bond_List( system->N, system->N_max, bond_top, lists[BONDS] );
}
- else if ( alloc == TRUE || lists[BONDS]->total_intrs < num_bonds )
+ else if ( realloc == 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] );
@@ -947,7 +947,7 @@ static void Init_Lists( reax_system *system, control_params *control,
#endif
/* 3bodies list */
- if ( first_run == TRUE )
+ if ( lists[THREE_BODIES]->allocated == FALSE )
{
Make_List( num_bonds, num_bonds, num_3body, TYP_THREE_BODY, lists[THREE_BODIES] );
}
@@ -1265,7 +1265,7 @@ 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,
- int output_enabled, int first_run, int alloc )
+ int output_enabled, int first_run, int realloc )
{
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(control)
@@ -1281,11 +1281,11 @@ void Initialize( reax_system *system, control_params *control,
Init_System( system, control, data, first_run );
- Init_Simulation_Data( system, control, data, out_control, Evolve, alloc );
+ Init_Simulation_Data( system, control, data, out_control, Evolve, realloc );
- Init_Workspace( system, control, workspace, alloc );
+ Init_Workspace( system, control, workspace, realloc );
- Init_Lists( system, control, data, workspace, lists, out_control, first_run, alloc );
+ Init_Lists( system, control, data, workspace, lists, out_control, first_run, realloc );
Init_Out_Controls( system, control, workspace, out_control, output_enabled );
@@ -1317,11 +1317,11 @@ static void Finalize_System( reax_system *system, control_params *control,
{
sfree( reax->gp.l, "Finalize_System::reax->gp.l" );
- for ( i = 0; i < reax->num_atom_types; i++ )
+ for ( i = 0; i < reax->max_num_atom_types; i++ )
{
- for ( j = 0; j < reax->num_atom_types; j++ )
+ for ( j = 0; j < reax->max_num_atom_types; j++ )
{
- for ( k = 0; k < reax->num_atom_types; k++ )
+ for ( k = 0; k < reax->max_num_atom_types; k++ )
{
sfree( reax->fbp[i][j][k], "Finalize_System::reax->fbp[i][j][k]" );
}
@@ -1610,7 +1610,7 @@ static void Finalize_Workspace( reax_system *system, control_params *control,
static void Finalize_Lists( control_params *control, reax_list **lists )
{
Delete_List( TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
- if ( control->hbond_cut > 0.0 )
+ if ( control->hbond_cut > 0.0 && lists[HBONDS]->index != NULL )
{
Delete_List( TYP_HBOND, lists[HBONDS] );
}
@@ -1624,7 +1624,7 @@ static void Finalize_Lists( control_params *control, reax_list **lists )
}
-static void Finalize_Out_Controls( reax_system *system, control_params *control,
+void Finalize_Out_Controls( reax_system *system, control_params *control,
static_storage *workspace, output_controls *out_control )
{
if ( out_control->write_steps > 0 )
diff --git a/sPuReMD/src/init_md.h b/sPuReMD/src/init_md.h
index 26b3b0d5..60c7f489 100644
--- a/sPuReMD/src/init_md.h
+++ b/sPuReMD/src/init_md.h
@@ -29,6 +29,9 @@ void Initialize( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls*, evolve_function*,
int, int, int );
+void Finalize_Out_Controls( reax_system *, control_params *,
+ static_storage *, output_controls * );
+
void Finalize( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls*, int, int );
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index c4fc219b..f23591a3 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -71,7 +71,7 @@ void Velocity_Verlet_NVE( reax_system *system, control_params *control,
lists, out_control );
}
- Compute_Forces( system, control, data, workspace, lists, out_control );
+ Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
for ( i = 0; i < system->N; i++ )
{
@@ -123,7 +123,7 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
}
Compute_Forces( system, control, data, workspace,
- lists, out_control );
+ lists, out_control, FALSE );
/* velocity verlet, 2nd part */
for ( i = 0; i < system->N; i++ )
@@ -213,7 +213,7 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system* system, control_params*
}
/* Calculate Forces at time (t + dt) */
- Compute_Forces( system, control, data, workspace, lists, out_control );
+ Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
/* Compute iteration constants for each atom's velocity */
for ( i = 0; i < system->N; ++i )
@@ -302,7 +302,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
lists, out_control );
}
- Compute_Forces( system, control, data, workspace, lists, out_control );
+ Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
/* velocity verlet, 2nd part */
for ( i = 0; i < system->N; i++ )
@@ -410,7 +410,7 @@ void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
lists, out_control );
}
- Compute_Forces( system, control, data, workspace, lists, out_control );
+ Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
/* velocity verlet, 2nd part */
for ( i = 0; i < system->N; i++ )
@@ -481,9 +481,8 @@ void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
/* BELOW FUNCTIONS ARE NOT BEING USED ANYMORE! */
/************************************************/
#if defined(ANISOTROPIC)
-void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
- control_params* control, simulation_data *data,
- static_storage *workspace, reax_list **lists,
+void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* control,
+ simulation_data *data, static_storage *workspace, reax_list **lists,
output_controls *out_control )
{
int i;
@@ -528,7 +527,7 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
// fprintf( out_control->log, "qeq-" );
// fflush( out_control->log );
- Compute_Forces( system, control, data, workspace, lists, out_control );
+ Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
fprintf(out_control->log, "forces\n");
fflush( out_control->log );
@@ -555,9 +554,8 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
}
-void Velocity_Verlet_Isotropic_NPT( reax_system* system,
- control_params* control, simulation_data *data,
- static_storage *workspace, reax_list **lists,
+void Velocity_Verlet_Isotropic_NPT( reax_system* system, control_params* control,
+ simulation_data *data, static_storage *workspace, reax_list **lists,
output_controls *out_control )
{
int i, itr;
@@ -633,8 +631,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
// fprintf( out_control->log, "qeq-" );
// fflush( out_control->log );
- Compute_Forces( system, control, data, workspace, lists,
- out_control );
+ Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
fprintf(out_control->log, "forces\n");
fflush( out_control->log );
diff --git a/sPuReMD/src/integrate.h b/sPuReMD/src/integrate.h
index 57935aa2..7207fa8f 100644
--- a/sPuReMD/src/integrate.h
+++ b/sPuReMD/src/integrate.h
@@ -26,27 +26,27 @@
#include "reax_types.h"
-void Velocity_Verlet_NVE( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+void Velocity_Verlet_NVE( reax_system *, control_params *, simulation_data *,
+ static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Berendsen_NVT( reax_system* , control_params* ,
+void Velocity_Verlet_Berendsen_NVT( reax_system *, control_params *,
simulation_data *, static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls* );
+void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system *, control_params *,
+ simulation_data *, static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls* );
+void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system *, control_params *,
+ simulation_data *, static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls* );
+void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system *, control_params *,
+ simulation_data *, static_storage *, reax_list **, output_controls * );
#if defined(ANISOTROPIC)
-void Velocity_Verlet_Nose_Hoover_NVT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls* );
+void Velocity_Verlet_Nose_Hoover_NVT( reax_system *, control_params *,
+ simulation_data *, static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Isotropic_NPT( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls* );
+void Velocity_Verlet_Isotropic_NPT( reax_system *, control_params *,
+ simulation_data *, static_storage *, reax_list **, output_controls * );
#endif
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 29f4eb8e..786a0413 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -162,7 +162,7 @@ void Sort_Matrix_Rows( sparse_matrix * const A )
* A has non-zero diagonals
* Each row of A has at least one non-zero (i.e., no rows with all zeros) */
static void compute_full_sparse_matrix( const sparse_matrix * const A,
- sparse_matrix ** A_full )
+ sparse_matrix ** A_full, int realloc )
{
int count, i, pj;
sparse_matrix *A_t;
@@ -171,7 +171,7 @@ static void compute_full_sparse_matrix( const sparse_matrix * const A,
{
Allocate_Matrix( A_full, A->n, A->n_max, 2 * A->m - A->n );
}
- else if ( (*A_full)->m < 2 * A->m - A->n )
+ else if ( (*A_full)->m < 2 * A->m - A->n || realloc == TRUE )
{
Deallocate_Matrix( *A_full );
Allocate_Matrix( A_full, A->n, A->n_max, 2 * A->m - A->n );
@@ -221,9 +221,10 @@ static void compute_full_sparse_matrix( const sparse_matrix * const A,
* Assumptions:
* A has non-zero diagonals
* Each row of A has at least one non-zero (i.e., no rows with all zeros) */
-void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix ** A_full,
- sparse_matrix ** A_spar_patt, sparse_matrix **A_spar_patt_full,
- sparse_matrix ** A_app_inv, const real filter )
+void setup_sparse_approx_inverse( const sparse_matrix * const A,
+ sparse_matrix ** A_full, sparse_matrix ** A_spar_patt,
+ sparse_matrix **A_spar_patt_full, sparse_matrix ** A_app_inv,
+ const real filter, int realloc )
{
int i, pj, size;
int left, right, k, p, turn;
@@ -235,7 +236,7 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
{
Allocate_Matrix( A_spar_patt, A->n, A->n_max, A->m );
}
- else if ( ((*A_spar_patt)->m) < (A->m) )
+ else if ( ((*A_spar_patt)->m) < (A->m) || realloc == TRUE )
{
Deallocate_Matrix( *A_spar_patt );
Allocate_Matrix( A_spar_patt, A->n, A->n_max, A->m );
@@ -345,24 +346,24 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
}
(*A_spar_patt)->start[A->n] = size;
- compute_full_sparse_matrix( A, A_full );
- compute_full_sparse_matrix( *A_spar_patt, A_spar_patt_full );
+ compute_full_sparse_matrix( A, A_full, realloc );
+ compute_full_sparse_matrix( *A_spar_patt, A_spar_patt_full, realloc );
if ( *A_app_inv == NULL )
{
/* 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)->n_max,
- (*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 )
+ else if ( (*A_app_inv)->m < A->m || realloc == TRUE )
{
Deallocate_Matrix( *A_app_inv );
/* 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)->n_max,
- (*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" );
@@ -706,7 +707,7 @@ real ILU( const sparse_matrix * const A, sparse_matrix * const L,
start = Get_Time( );
- compute_full_sparse_matrix( A, &A_full );
+ compute_full_sparse_matrix( A, &A_full, FALSE );
Ltop = 0;
Utop = 0;
@@ -812,7 +813,7 @@ real ILU( const sparse_matrix * const A, sparse_matrix * const L,
* droptol: row-wise tolerances used for dropping
* L / U: (output) triangular matrices, A \approx LU, CSR format */
real ILUT( const sparse_matrix * const A, const real * const droptol,
- sparse_matrix * const L, sparse_matrix * const U )
+ sparse_matrix * const L, sparse_matrix * const U )
{
unsigned int i, k, pj, Ltop, Utop, *nz_mask;
real *w, start;
@@ -825,7 +826,7 @@ real ILUT( const sparse_matrix * const A, const real * const droptol,
w = smalloc( sizeof(real) * A->n, "ILUT::w" );
nz_mask = smalloc( sizeof(unsigned int) * A->n, "ILUT::nz_mask" );
- compute_full_sparse_matrix( A, &A_full );
+ compute_full_sparse_matrix( A, &A_full, FALSE );
Ltop = 0;
Utop = 0;
@@ -969,7 +970,7 @@ real ILUT( const sparse_matrix * const A, const real * const droptol,
* droptol: row-wise tolerances used for dropping
* L / U: (output) triangular matrices, A \approx LU, CSR format */
real ILUTP( const sparse_matrix * const A, const real * const droptol,
- sparse_matrix * const L, sparse_matrix * const U )
+ sparse_matrix * const L, sparse_matrix * const U )
{
unsigned int i, k, pj, Ltop, Utop, *nz_mask, *perm, *perm_inv, pivot_j;
real *w, start, pivot_val;
@@ -983,7 +984,7 @@ real ILUTP( const sparse_matrix * const A, const real * const droptol,
perm = smalloc( sizeof(int) * A->n, "ILUTP::perm" );
perm_inv = smalloc( sizeof(int) * A->n, "ILUTP::perm_inv" );
- compute_full_sparse_matrix( A, &A_full );
+ compute_full_sparse_matrix( A, &A_full, FALSE );
Ltop = 0;
Utop = 0;
@@ -2628,19 +2629,19 @@ static void permute_matrix( const static_storage * const workspace,
*/
void setup_graph_coloring( const control_params * const control,
const static_storage * const workspace, const sparse_matrix * const H,
- sparse_matrix ** H_full, sparse_matrix ** H_p )
+ sparse_matrix ** H_full, sparse_matrix ** H_p, int realloc )
{
if ( *H_p == NULL )
{
Allocate_Matrix( H_p, H->n, H->n_max, H->m );
}
- else if ( (*H_p)->m < H->m )
+ else if ( (*H_p)->m < H->m || realloc == TRUE )
{
Deallocate_Matrix( *H_p );
Allocate_Matrix( H_p, H->n, H->n_max, H->m );
}
- compute_full_sparse_matrix( H, H_full );
+ compute_full_sparse_matrix( H, H_full, realloc );
graph_coloring( control, (static_storage *) workspace, *H_full, LOWER );
sort_rows_by_colors( workspace, (*H_full)->n, LOWER );
diff --git a/sPuReMD/src/lin_alg.h b/sPuReMD/src/lin_alg.h
index aa9defaf..1edf57b7 100644
--- a/sPuReMD/src/lin_alg.h
+++ b/sPuReMD/src/lin_alg.h
@@ -33,7 +33,7 @@ typedef enum
void Sort_Matrix_Rows( sparse_matrix * const );
void setup_sparse_approx_inverse( const sparse_matrix * const, sparse_matrix **,
- sparse_matrix **, sparse_matrix **, sparse_matrix **, const real );
+ sparse_matrix **, sparse_matrix **, sparse_matrix **, const real, int );
int Estimate_LU_Fill( const sparse_matrix * const, const real * const );
@@ -84,7 +84,7 @@ void jacobi_iter( const static_storage * const,
void setup_graph_coloring( const control_params * const,
const static_storage * const, const sparse_matrix * const,
- sparse_matrix **, sparse_matrix ** );
+ sparse_matrix **, sparse_matrix **, int );
int GMRES( const static_storage * const, const control_params * const,
simulation_data * const, const sparse_matrix * const,
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index f5414e93..8f69ca35 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -32,6 +32,14 @@ void Make_List( int n, int n_max, int total_intrs, int type, reax_list* l )
assert( total_intrs > 0 );
assert( l != NULL );
+ if ( l->allocated == TRUE )
+ {
+ fprintf( stderr, "[WARNING] attempted to allocate list which was already allocated."
+ " Returning without allocation...\n" );
+ return;
+ }
+
+ l->allocated = TRUE;
l->n = n;
l->n_max = n_max;
l->total_intrs = total_intrs;
@@ -91,6 +99,20 @@ void Make_List( int n, int n_max, int total_intrs, int type, reax_list* l )
void Delete_List( int type, reax_list* l )
{
+ assert( l != NULL );
+
+ if ( l->allocated == FALSE )
+ {
+ fprintf( stderr, "[WARNING] attempted to free list which was not allocated."
+ " Returning without deallocation...\n" );
+ return;
+ }
+
+ l->allocated = FALSE;
+ l->n = 0;
+ l->n_max = 0;
+ l->total_intrs = 0;
+
sfree( l->index, "Delete_List::l->index" );
sfree( l->end_index, "Delete_List::l->end_index" );
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 58decd75..a551ce07 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -449,11 +449,10 @@ typedef void (*print_interaction)(reax_system*, control_params*, simulation_data
/* function pointer for evolving the atomic system (i.e., updating the positions)
* given the pre-computed forces from the prescribed interactions */
typedef void (*evolve_function)( reax_system*, control_params*,
- simulation_data*, static_storage*,
- reax_list**, output_controls* );
+ simulation_data*, static_storage*, reax_list**, output_controls* );
/* function pointer for a callback function to be triggered after
* completion of a simulation step -- useful for, e.g., the Python wrapper */
-typedef void (*callback_function)( reax_atom*, simulation_data*, reax_list** );
+typedef void (*callback_function)( int, reax_atom*, simulation_data* );
/* function pointer for writing trajectory file header */
typedef int (*write_header_function)( reax_system*, control_params*,
static_storage*, output_controls* );
@@ -509,8 +508,13 @@ typedef int (*write_function)( FILE *, const char *, ... );
* l[38] = p_coa3 */
struct global_parameters
{
+ /* num. global parameters in the force field paramater file for this simulation */
int n_global;
+ /* max. num. global parameters in the force field paramater file across all simulations */
+ int max_n_global;
+ /* parameters, see list in comment above for mapping */
real* l;
+ /* van der Waals interaction type */
int vdw_type;
};
@@ -732,12 +736,21 @@ struct four_body_header
struct reax_interaction
{
+ /* num. atom types in force field parameter file for this simulation */
int num_atom_types;
+ /* max. num. atom types in force field parameter file across all simulations */
+ int max_num_atom_types;
+ /* */
global_parameters gp;
+ /* */
single_body_parameters *sbp;
+ /* */
two_body_parameters **tbp;
+ /* */
three_body_header ***thbp;
+ /* */
hbond_parameters ***hbp;
+ /* */
four_body_header ****fbp;
};
@@ -1617,6 +1630,8 @@ struct static_storage
/* interaction lists */
struct reax_list
{
+ /* 0 if struct members are NOT allocated, 1 otherwise */
+ int allocated;
/* num. active entries in list for this simulation */
int n;
/* max. num. entries in list across all simulations */
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 7fc793ec..dfd1f869 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -179,7 +179,7 @@ void* setup( const char * const geo_file, const char * const ffield_file,
{
spmd_handle->lists[i] = smalloc( sizeof(reax_list),
"Setup::spmd_handle->lists[i]" );
-// spmd_handle->lists[i]->allocated = FALSE;
+ spmd_handle->lists[i]->allocated = FALSE;
}
spmd_handle->out_control = smalloc( sizeof(output_controls),
"Setup::spmd_handle->out_control" );
@@ -258,7 +258,8 @@ int simulate( const void * const handle )
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
+ spmd_handle->realloc );
Compute_Kinetic_Energy( spmd_handle->system, spmd_handle->data );
@@ -301,8 +302,8 @@ int simulate( const void * const handle )
if ( spmd_handle->callback != NULL )
{
- spmd_handle->callback( spmd_handle->system->atoms, spmd_handle->data,
- spmd_handle->lists );
+ spmd_handle->callback( spmd_handle->system->N, spmd_handle->system->atoms,
+ spmd_handle->data );
}
//}
@@ -352,8 +353,8 @@ int simulate( const void * const handle )
if ( spmd_handle->callback != NULL )
{
- spmd_handle->callback( spmd_handle->system->atoms, spmd_handle->data,
- spmd_handle->lists );
+ spmd_handle->callback( spmd_handle->system->N, spmd_handle->system->atoms,
+ spmd_handle->data );
}
}
@@ -431,6 +432,13 @@ int reset( const void * const handle, const char * const geo_file,
{
spmd_handle = (spuremd_handle*) handle;
+ /* close files used in previous simulation */
+ if ( spmd_handle->output_enabled == TRUE )
+ {
+ Finalize_Out_Controls( spmd_handle->system, spmd_handle->control,
+ spmd_handle->workspace, spmd_handle->out_control );
+ }
+
spmd_handle->first_run = FALSE;
spmd_handle->realloc = FALSE;
spmd_handle->data->sim_id++;
@@ -460,24 +468,64 @@ int reset( const void * const handle, const char * const geo_file,
}
-/* Getter for atom info
+/* Getter for atom positions
*
* handle: pointer to wrapper struct with top-level data structures
+ * pos_x: x-coordinate of atom positions (allocated by caller)
+ * pos_y: y-coordinate of atom positions (allocated by caller)
+ * pos_z: z-coordinate of atom positions (allocated by caller)
*/
-reax_atom* get_atoms( const void * const handle )
+int get_atom_positions( const void * const handle, double * const pos_x,
+ double * const pos_y, double * const pos_z )
{
+ int i, ret;
spuremd_handle *spmd_handle;
- reax_atom *atoms;
- atoms = NULL;
+ ret = SPUREMD_FAILURE;
if ( handle != NULL )
{
spmd_handle = (spuremd_handle*) handle;
- atoms = spmd_handle->system->atoms;
+
+ for ( i = 0; i < spmd_handle->system->N; ++i )
+ {
+ pos_x[i] = spmd_handle->system->atoms[i].x[0];
+ pos_y[i] = spmd_handle->system->atoms[i].x[1];
+ pos_z[i] = spmd_handle->system->atoms[i].x[2];
+ }
+
+ ret = SPUREMD_SUCCESS;
}
- return atoms;
+ return ret;
+}
+
+
+/* Getter for atom charges
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * q: atom charges (allocated by caller)
+ */
+int get_atom_charges( const void * const handle, double * const q )
+{
+ int i, ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ for ( i = 0; i < spmd_handle->system->N; ++i )
+ {
+ q[i] = spmd_handle->system->atoms[i].q;
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
}
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index 980d3448..d84c1875 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -45,7 +45,10 @@ 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 get_atom_positions( const void * const, double * const,
+ double * const, double * const );
+
+int get_atom_charges( const void * const, double * const );
int set_output_enabled( const void * const, const int );
diff --git a/tools/driver.py b/tools/driver.py
index c00f9063..ea6e281c 100644
--- a/tools/driver.py
+++ b/tools/driver.py
@@ -193,6 +193,7 @@ class ReaxTiming(Structure):
class SimulationData(Structure):
_fields_ = [
+ ("sim_id", c_int),
("step", c_int),
("prev_step", c_int),
("time", c_double),
@@ -400,6 +401,10 @@ if __name__ == '__main__':
cleanup.argtypes = [c_void_p]
cleanup.restype = c_int
+ reset = lib.reset
+ reset.argtypes = [c_void_p, c_char_p, c_char_p, c_char_p]
+ reset.restype = c_int
+
get_atoms = lib.get_atoms
get_atoms.argtypes = [c_void_p]
get_atoms.restype = POINTER(ReaxAtom)
@@ -435,14 +440,14 @@ if __name__ == '__main__':
with conn:
conn.execute("INSERT INTO system_properties VALUES (?,?,?,?,?,?,?,?)",
- (0, data[0].step, data[0].E_Tot, data[0].E_Pot, data[0].E_Kin,
+ (data[0].sim_id, data[0].step, data[0].E_Tot, data[0].E_Pot, data[0].E_Kin,
data[0].therm.T, 0.0, data[0].iso_bar.P))
# MISSING: ID, system->box.volume
if record_potential:
with conn:
conn.execute("INSERT INTO potential VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?)",
- (0, data[0].step, data[0].E_BE, data[0].E_Ov + data[0].E_Un,
+ (data[0].sim_id, data[0].step, data[0].E_BE, data[0].E_Ov + data[0].E_Un,
data[0].E_Lp, data[0].E_Ang + data[0].E_Pen, data[0].E_Coa,
data[0].E_HB, data[0].E_Tor, data[0].E_Con, data[0].E_vdW,
data[0].E_Ele, data[0].E_Pol))
@@ -451,13 +456,13 @@ if __name__ == '__main__':
with conn:
for i in range(6540):
conn.execute("INSERT INTO trajectory VALUES (?,?,?,?,?,?,?)",
- (0, data[0].step, i, atoms[i].x[0], atoms[i].x[1],
+ (data[0].sim_id, data[0].step, i, atoms[i].x[0], atoms[i].x[1],
atoms[i].x[2], atoms[i].q))
if record_performance:
with conn:
conn.execute("INSERT INTO performance VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)",
- (0, data[0].step, data[0].timing.total, data[0].timing.nbrs, data[0].timing.init_forces,
+ (data[0].sim_id, data[0].step, data[0].timing.total, data[0].timing.nbrs, data[0].timing.init_forces,
data[0].timing.bonded, data[0].timing.nonb, data[0].timing.cm,
data[0].timing.cm_sort_mat_rows, data[0].timing.cm_solver_iters,
data[0].timing.cm_solver_pre_comp,
@@ -479,6 +484,22 @@ if __name__ == '__main__':
atoms = get_atoms(handle)
+ print()
+ print("{0:9}|{1:24}|{2:24}|{3:24}|{4:24}".format("Atom Num", "x-Position", "y-Position", "z-Position", "Charge"))
+ for i in range(10):
+ print("{0:9d} {1:24.15f} {2:24.15f} {3:24.15f} {4:24.15f}".format(
+ i + 1, atoms[i].x[0], atoms[i].x[1], atoms[i].x[2], atoms[i].q))
+
+ ret = reset(handle, b"data/benchmarks/silica/silica_6000.pdb",
+ b"data/benchmarks/silica/ffield-bio",
+ b"environ/control_silica")
+
+ print()
+ print("{0:24}|{1:24}|{2:24}".format("Total Energy", "Kinetic Energy", "Potential Energy"))
+ ret = simulate(handle)
+
+ atoms = get_atoms(handle)
+
print()
print("{0:9}|{1:24}|{2:24}|{3:24}|{4:24}".format("Atom Num", "x-Position", "y-Position", "z-Position", "Charge"))
for i in range(10):
--
GitLab