sPuReMD: separate total energy computation and fix Python driver bug with...

sPuReMD: separate total energy computation and fix Python driver bug with structure definition for ctypes. Allow disabling file logging during simulation.

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 );
-    }
 }
 
 

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;
+}

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

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.
  *

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* );

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)