diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 2ddd9eac7e6a3c57913a38dcb651eeb9dbe6678b..447e7cc326a7b10e4f0c60aca94bfe16829c30be 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -577,36 +577,9 @@ void Output_Results( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
{
- int i, type_i;
- real e_pol, q, f_update;
+ real f_update;
real t_elapsed = 0;
- /* Compute Polarization Energy */
- e_pol = 0.0;
-
-#ifdef _OPENMP
- #pragma omp parallel for default(none) private(q, type_i) shared(system) \
- reduction(+: e_pol) schedule(static)
-#endif
- for ( i = 0; i < system->N; i++ )
- {
- q = system->atoms[i].q;
- type_i = system->atoms[i].type;
-
- e_pol += ( system->reaxprm.sbp[ type_i ].chi * q +
- (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) ) *
- KCALpMOL_to_EV;
- }
-
- data->E_Pol = e_pol;
-
- data->E_Pot = data->E_BE + data->E_Ov + data->E_Un + data->E_Lp +
- data->E_Ang + data->E_Pen + data->E_Coa + data->E_HB +
- data->E_Tor + data->E_Con + data->E_vdW + data->E_Ele + data->E_Pol;
-
-
- data->E_Tot = data->E_Pot + E_CONV * data->E_Kin;
-
/* output energies if it is the time */
if ( out_control->energy_update_freq > 0 &&
data->step % out_control->energy_update_freq == 0 )
@@ -724,24 +697,6 @@ void Output_Results( reax_system *system, control_params *control,
//t_elapsed = Get_Timing_Info( t_start );
//fprintf(stdout, "append_frame took %.6f seconds\n", t_elapsed );
}
-
- if ( IS_NAN_REAL(data->E_Pol) )
- {
- fprintf( stderr, "[ERROR] NaN detected for polarization energy. Terminating...\n" );
- exit( NUMERIC_BREAKDOWN );
- }
-
- if ( IS_NAN_REAL(data->E_Pot) )
- {
- fprintf( stderr, "[ERROR] NaN detected for potential energy. Terminating...\n" );
- exit( NUMERIC_BREAKDOWN );
- }
-
- if ( IS_NAN_REAL(data->E_Tot) )
- {
- fprintf( stderr, "[ERROR] NaN detected for total energy. Terminating...\n" );
- exit( NUMERIC_BREAKDOWN );
- }
}
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index f18e71d521d5e0195db3957bb29bfb760d4b94bb..cb3e7418a86f6fec84ae147d7c19b91e974da5d0 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -74,6 +74,8 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
rvec_ScaledAdd( system->atoms[i].v, -1., cross );
}
}
+
+ Compute_Total_Energy( system, data );
}
@@ -178,8 +180,6 @@ void* setup( const char * const geo_file, const char * const ffield_file,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
- //TODO: if errors detected, set handle to NULL to indicate failure
-
return (void*) spmd_handle;
}
@@ -235,11 +235,18 @@ int simulate( const void * const handle )
Compute_Kinetic_Energy( spmd_handle->system, spmd_handle->data );
+ Compute_Total_Energy( spmd_handle->system, 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 );
}
+ if ( spmd_handle->callback != NULL )
+ {
+ spmd_handle->callback( spmd_handle->system->atoms, spmd_handle->data,
+ spmd_handle->lists );
+ }
++spmd_handle->data->step;
//}
@@ -262,11 +269,13 @@ int simulate( const void * const handle )
{
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+
Analysis( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
}
steps = spmd_handle->data->step - spmd_handle->data->prev_steps;
+
if ( steps && spmd_handle->out_control->restart_freq &&
steps % spmd_handle->out_control->restart_freq == 0 &&
spmd_handle->output_enabled == TRUE )
@@ -350,3 +359,21 @@ reax_atom* get_atoms( const void * const handle )
return atoms;
}
+
+
+int set_output_enabled( const void * const handle, const int enabled )
+{
+ int ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+ spmd_handle->output_enabled = enabled;
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index 3a78242bd768e2674c727ac882f74b69b2a5e336..679aeaf3bdd3a9a14a44099994100d059fbee19e 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -22,11 +22,11 @@
#ifndef __SPUREMD_H_
#define __SPUREMD_H_
-#define SPUREMD_SUCCESS (0)
-#define SPUREMD_FAILURE (-1)
+#include "mytypes.h"
-#include "mytypes.h"
+#define SPUREMD_SUCCESS (0)
+#define SPUREMD_FAILURE (-1)
#ifdef __cplusplus
@@ -44,6 +44,8 @@ int cleanup( const void * const );
reax_atom* get_atoms( const void * const );
+int set_output_enabled( const void * const, const int );
+
#ifdef __cplusplus
}
#endif
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index ad98d171927ac81f54a6ce2133a198ed9e6b8f3f..cff459fd897c8b46df06ff3d01ce93687365e37a 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -210,6 +210,57 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
}
+void Compute_Total_Energy( reax_system* system, simulation_data* data )
+{
+ int i, type_i;
+ real e_pol, q;
+
+ /* Compute Polarization Energy */
+ e_pol = 0.0;
+
+#ifdef _OPENMP
+ #pragma omp parallel for default(none) private(q, type_i) shared(system) \
+ reduction(+: e_pol) schedule(static)
+#endif
+ for ( i = 0; i < system->N; i++ )
+ {
+ q = system->atoms[i].q;
+ type_i = system->atoms[i].type;
+
+ e_pol += ( system->reaxprm.sbp[ type_i ].chi * q +
+ (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) ) *
+ KCALpMOL_to_EV;
+ }
+
+ data->E_Pol = e_pol;
+
+ data->E_Pot = data->E_BE + data->E_Ov + data->E_Un + data->E_Lp +
+ data->E_Ang + data->E_Pen + data->E_Coa + data->E_HB +
+ data->E_Tor + data->E_Con + data->E_vdW + data->E_Ele + data->E_Pol;
+
+
+ data->E_Tot = data->E_Pot + E_CONV * data->E_Kin;
+
+ if ( IS_NAN_REAL(data->E_Pol) )
+ {
+ fprintf( stderr, "[ERROR] NaN detected for polarization energy. Terminating...\n" );
+ exit( NUMERIC_BREAKDOWN );
+ }
+
+ if ( IS_NAN_REAL(data->E_Pot) )
+ {
+ fprintf( stderr, "[ERROR] NaN detected for potential energy. Terminating...\n" );
+ exit( NUMERIC_BREAKDOWN );
+ }
+
+ if ( IS_NAN_REAL(data->E_Tot) )
+ {
+ fprintf( stderr, "[ERROR] NaN detected for total energy. Terminating...\n" );
+ exit( NUMERIC_BREAKDOWN );
+ }
+}
+
+
/* IMPORTANT: This function assumes that current kinetic energy and
* the center of mass of the system is already computed before.
*
diff --git a/sPuReMD/src/system_props.h b/sPuReMD/src/system_props.h
index 996e6a452299d25d987c620d5769135caf71760f..de65bf666b619732f56dc2b47049a57b268d7aa2 100644
--- a/sPuReMD/src/system_props.h
+++ b/sPuReMD/src/system_props.h
@@ -33,6 +33,8 @@ void Compute_Center_of_Mass( reax_system*, simulation_data*, FILE* );
void Compute_Kinetic_Energy( reax_system*, simulation_data* );
+void Compute_Total_Energy( reax_system*, simulation_data* );
+
void Compute_Pressure_Isotropic( reax_system*, control_params*, simulation_data*, output_controls* );
void Compute_Pressure_Isotropic_Klein( reax_system*, simulation_data* );
diff --git a/tools/driver.py b/tools/driver.py
index 12a4a0882981ef73d8b252dbb360505e830d35f4..59a8c4a8ba7b6eda19fe08b8dcb9e61fd18e38fc 100644
--- a/tools/driver.py
+++ b/tools/driver.py
@@ -227,6 +227,7 @@ class SimulationData(Structure):
class ReaxAtom(Structure):
_fields_ = [
+ ("type", c_int),
("name", c_char * 8),
("x", c_double * 3),
("v", c_double * 3),
@@ -264,12 +265,18 @@ if __name__ == '__main__':
setup_callback = lib.setup_callback
setup_callback.restype = c_int
+ set_output_enabled = lib.set_output_enabled
+ set_output_enabled.argtypes = [c_void_p, c_int]
+ set_output_enabled.restype = c_int
+
handle = setup(b"data/benchmarks/water/water_6540.pdb",
b"data/benchmarks/water/ffield.water",
b"environ/param.gpu.water")
ret = setup_callback(handle, CALLBACKFUNC(get_simulation_step_results))
+ ret = set_output_enabled(handle, c_int(0))
+
print("{0:24}|{1:24}|{2:24}".format("Total Energy", "Kinetic Energy", "Potential Energy"))
ret = simulate(handle)