From c64018078904bc649b061e23a79ba08445e45e7c Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Wed, 7 Feb 2018 14:47:47 -0500
Subject: [PATCH] sPuReMD: update library interface and driver. Add prototype
 driver in Python.

---
 sPuReMD/src/driver.c    |  27 ++++--
 sPuReMD/src/geo_tools.c |   6 +-
 sPuReMD/src/geo_tools.h |   6 +-
 sPuReMD/src/mytypes.h   |  20 ++++
 sPuReMD/src/restart.c   |   4 +-
 sPuReMD/src/restart.h   |   4 +-
 sPuReMD/src/spuremd.c   | 204 ++++++++++++++++++++++++++--------------
 sPuReMD/src/spuremd.h   |  13 ++-
 sPuReMD/src/tool_box.c  |   2 -
 tools/driver.py         | 142 ++++++++++++++++++++++++++++
 10 files changed, 328 insertions(+), 100 deletions(-)
 create mode 100644 tools/driver.py

diff --git a/sPuReMD/src/driver.c b/sPuReMD/src/driver.c
index 6d2ff6ad..8e48dc1b 100644
--- a/sPuReMD/src/driver.c
+++ b/sPuReMD/src/driver.c
@@ -20,22 +20,24 @@
   <http://www.gnu.org/licenses/>.
   ----------------------------------------------------------------------*/
 
-#include "mytypes.h"
+#include <stdio.h>
+#include <stdlib.h>
 
 #include "spuremd.h"
 
+#define INVALID_INPUT (-1)
+
 
 static void usage( char * argv[] )
 {
-    fprintf( stderr, "usage: ./%s geometry ffield control\n", argv[0] );
+    fprintf( stderr, "usage: ./%s geometry_file force_field_file control_file\n", argv[0] );
 }
 
 
 int main( int argc, char* argv[] )
 {
-    reax_system system;
-    control_params control;
-    simulation_data data;
+    void *handle;
+    int ret;
 
     if ( argc != 4 )
     {
@@ -43,11 +45,18 @@ int main( int argc, char* argv[] )
         exit( INVALID_INPUT );
     }
 
-    Setup( argv + 1, &system, &control, &data );
+    handle = setup( argv[1], argv[2], argv[3] );
+    ret = SPUREMD_FAILURE;
 
-    Run( &system, &control, &data, TRUE );
+    if ( handle != NULL )
+    {
+        ret = simulate( handle );
+    }
 
-    Cleanup( &system, &control, &data, TRUE );
+    if ( ret == SPUREMD_SUCCESS )
+    {
+        ret = cleanup( handle );
+    }
 
-    return 0;
+    return (ret == SPUREMD_SUCCESS) ? 0 : (-1);
 }
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 69675645..740a9648 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -55,7 +55,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
 }
 
 
-char Read_Geo( char* geo_file, reax_system* system, control_params *control,
+char Read_Geo( const char * const geo_file, reax_system* system, control_params *control,
         simulation_data *data, static_storage *workspace )
 {
 
@@ -221,7 +221,7 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
 }
 
 
-char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
+char Read_PDB( const char * const pdb_file, reax_system* system, control_params *control,
         simulation_data *data, static_storage *workspace )
 {
 
@@ -538,7 +538,7 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
 }
 
 
-char Read_BGF( char* bgf_file, reax_system* system, control_params *control,
+char Read_BGF( const char * const bgf_file, reax_system* system, control_params *control,
         simulation_data *data, static_storage *workspace )
 {
     FILE *bgf;
diff --git a/sPuReMD/src/geo_tools.h b/sPuReMD/src/geo_tools.h
index 0339a166..5e70721d 100644
--- a/sPuReMD/src/geo_tools.h
+++ b/sPuReMD/src/geo_tools.h
@@ -29,7 +29,7 @@
 // CUSTOM ATOM: serial element name x y z
 #define CUSTOM_ATOM_FORMAT " %d %s %s %lf %lf %lf"
 
-char Read_Geo( char*, reax_system*, control_params*,
+char Read_Geo( const char * const, reax_system*, control_params*,
         simulation_data*, static_storage* );
 
 /* PDB format :
@@ -117,10 +117,10 @@ COLUMNS       DATA TYPE       FIELD         DEFINITION
 
 #define BGF_CRYSTX_FORMAT "%8s%11s%11s%11s%11s%11s%11s"
 
-char Read_PDB( char*, reax_system*, control_params*,
+char Read_PDB( const char * const, reax_system*, control_params*,
         simulation_data*, static_storage* );
 
-char Read_BGF( char*, reax_system*, control_params*,
+char Read_BGF( const char * const, reax_system*, control_params*,
         simulation_data*, static_storage* );
 
 char Write_PDB( reax_system*, reax_list*, simulation_data*,
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 023542b9..38d1e35c 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -1122,6 +1122,26 @@ typedef struct
 } LR_lookup_table;
 
 
+/* Handle for working with an instance of the sPuReMD library */
+typedef struct
+{
+    /* System info. struct pointer */
+    reax_system *system;
+    /* System struct pointer */
+    control_params *control;
+    /* Control parameters struct pointer */
+    simulation_data *data;
+    /* Internal workspace struct pointer */
+    static_storage *workspace;
+    /* Reax interaction list struct pointer */
+    reax_list *lists;
+    /* Output controls struct pointer */
+    output_controls *out_control;
+    /* TRUE if file I/O for simulation output enabled, FALSE otherwise */
+    int output_enabled;
+} spuremd_handle;
+
+
 /* Function pointer definitions */
 typedef void (*interaction_function)(reax_system*, control_params*,
         simulation_data*, static_storage*, reax_list**, output_controls*);
diff --git a/sPuReMD/src/restart.c b/sPuReMD/src/restart.c
index 82474a8b..23b00726 100644
--- a/sPuReMD/src/restart.c
+++ b/sPuReMD/src/restart.c
@@ -68,7 +68,7 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
 }
 
 
-void Read_Binary_Restart( char *fname, reax_system *system,
+void Read_Binary_Restart( const char * const fname, reax_system *system,
                           control_params *control, simulation_data *data,
                           static_storage *workspace )
 {
@@ -185,7 +185,7 @@ void Write_ASCII_Restart( reax_system *system, control_params *control,
 }
 
 
-void Read_ASCII_Restart( char *fname, reax_system *system,
+void Read_ASCII_Restart( const char * const fname, reax_system *system,
                          control_params *control, simulation_data *data,
                          static_storage *workspace )
 {
diff --git a/sPuReMD/src/restart.h b/sPuReMD/src/restart.h
index d0ea4ca8..363bb581 100644
--- a/sPuReMD/src/restart.h
+++ b/sPuReMD/src/restart.h
@@ -56,9 +56,9 @@ typedef struct
 void Write_Restart( reax_system*, control_params*,
                     simulation_data*, static_storage*, output_controls* );
 
-void Read_Binary_Restart( char*, reax_system*, control_params*,
+void Read_Binary_Restart( const char * const, reax_system*, control_params*,
                           simulation_data*, static_storage* );
-void Read_ASCII_Restart( char*, reax_system*, control_params*,
+void Read_ASCII_Restart( const char * const, reax_system*, control_params*,
                          simulation_data*, static_storage* );
 
 #endif
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 82958abe..5eca7537 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -19,6 +19,8 @@
   <http://www.gnu.org/licenses/>.
   ----------------------------------------------------------------------*/
 
+#include "spuremd.h"
+
 #include "mytypes.h"
 
 #include "analyze.h"
@@ -36,11 +38,6 @@
 #include "vector.h"
 
 
-static static_storage workspace;
-static reax_list *lists;
-static output_controls out_control;
-
-
 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 )
@@ -82,9 +79,9 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
 }
 
 
-void static Read_System( char * const geo_file,
-        char * const ffield_file,
-        char * const control_file,
+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,
@@ -151,94 +148,157 @@ void static Read_System( char * const geo_file,
 }
 
 
-int Setup( char ** args, reax_system * const system, control_params * const control,
-        simulation_data * const data )
+void* setup( const char * const geo_file, const char * const ffield_file,
+        const char * const control_file )
 {
-    lists = (reax_list*) smalloc( sizeof(reax_list) * LIST_N,
-           "Setup::lists" );
-
-    Read_System( args[0], args[1], args[2], system, control,
-            data, &workspace, &out_control );
-
-    return SUCCESS;
+    spuremd_handle *spmd_handle;
+
+    /* top-level allocation */
+    spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
+            "setup::spmd_handle" );
+
+    /* second-level allocations */
+    spmd_handle->system = (reax_system*) smalloc( sizeof(reax_system),
+           "Setup::spmd_handle->system" );
+    spmd_handle->control = (control_params*) smalloc( sizeof(control_params),
+           "Setup::spmd_handle->control" );
+    spmd_handle->data = (simulation_data*) smalloc( sizeof(simulation_data),
+           "Setup::spmd_handle->data" );
+    spmd_handle->workspace = (static_storage*) smalloc( sizeof(static_storage),
+           "Setup::spmd_handle->workspace" );
+    spmd_handle->lists = (reax_list*) smalloc( sizeof(reax_list) * LIST_N,
+           "Setup::spmd_handle->lists" );
+    spmd_handle->out_control = (output_controls*) smalloc( sizeof(output_controls),
+           "Setup::spmd_handle->out_control" );
+
+    spmd_handle->output_enabled = TRUE;
+
+    /* parse geometry file */
+    Read_System( geo_file, ffield_file, control_file,
+            spmd_handle->system, spmd_handle->control,
+            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;
 }
 
 
-int Run( reax_system * const system, control_params * const control, simulation_data * const data,
-      const int output_enabled )
+int simulate( const void * const handle )
 {
-    int steps;
+    int steps, ret;
     evolve_function Evolve;
+    spuremd_handle *spmd_handle;
 
-    Initialize( system, control, data, &workspace, &lists,
-            &out_control, &Evolve, output_enabled );
-
-    /* compute f_0 */
-    //if( control.restart == 0 ) {
-    Reset( system, control, data, &workspace, &lists );
-    Generate_Neighbor_Lists( system, control, data, &workspace, 
-            &lists, &out_control );
+    ret = SPUREMD_FAILURE;
 
-    //fprintf( stderr, "total: %.2f secs\n", data.timing.nbrs);
-    Compute_Forces( system, control, data, &workspace, &lists, &out_control );
-    Compute_Kinetic_Energy( system, data );
-    if ( output_enabled == TRUE )
-    {
-        Output_Results( system, control, data, &workspace, &lists, &out_control );
-    }
-    ++data->step;
-    //}
-    //
-    
-    for ( ; data->step <= control->nsteps; data->step++ )
+    if ( handle != NULL )
     {
-        if ( control->T_mode )
+        spmd_handle = (spuremd_handle*) 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 );
+
+        /* compute f_0 */
+        //if( control.restart == 0 ) {
+        Reset( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+                spmd_handle->workspace, &spmd_handle->lists );
+        Generate_Neighbor_Lists( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+                spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+
+        //fprintf( stderr, "total: %.2f secs\n", data.timing.nbrs);
+        Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+                spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+        Compute_Kinetic_Energy( spmd_handle->system, spmd_handle->data );
+        if ( spmd_handle->output_enabled == TRUE )
         {
-            Temperature_Control( control, data, &out_control );
+            Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+                    spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
         }
-
-        Evolve( system, control, data, &workspace, &lists, &out_control );
-        Post_Evolve( system, control, data, &workspace, &lists, &out_control );
-
-        if ( output_enabled == TRUE )
+        ++spmd_handle->data->step;
+        //}
+        
+        for ( ; spmd_handle->data->step <= spmd_handle->control->nsteps; spmd_handle->data->step++ )
         {
-            Output_Results( system, control, data, &workspace, &lists, &out_control );
-            Analysis( system, control, data, &workspace, &lists, &out_control );
+            if ( spmd_handle->control->T_mode )
+            {
+                Temperature_Control( spmd_handle->control, spmd_handle->data,
+                        spmd_handle->out_control );
+            }
+
+            Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+                    spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+            Post_Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+                    spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
+
+            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 );
+                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 )
+            {
+                Write_Restart( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+                        spmd_handle->workspace, spmd_handle->out_control );
+            }
         }
 
-        steps = data->step - data->prev_steps;
-        if ( steps && out_control.restart_freq &&
-                steps % out_control.restart_freq == 0 &&
-                output_enabled == TRUE )
+        if ( spmd_handle->out_control->write_steps > 0 && spmd_handle->output_enabled == TRUE )
         {
-            Write_Restart( system, control, data, &workspace, &out_control );
+            fclose( spmd_handle->out_control->trj );
+            Write_PDB( spmd_handle->system, &(spmd_handle->lists[BONDS]), spmd_handle->data,
+                    spmd_handle->control, spmd_handle->workspace, spmd_handle->out_control );
         }
-    }
 
-    if ( out_control.write_steps > 0 && output_enabled == TRUE )
-    {
-        fclose( out_control.trj );
-        Write_PDB( system, &(lists[BONDS]), data, control, &workspace, &out_control );
-    }
+        spmd_handle->data->timing.end = Get_Time( );
+        spmd_handle->data->timing.elapsed = Get_Timing_Info( spmd_handle->data->timing.start );
+        if ( spmd_handle->output_enabled == TRUE )
+        {
+            fprintf( spmd_handle->out_control->log, "total: %.2f secs\n", spmd_handle->data->timing.elapsed );
+        }
 
-    data->timing.end = Get_Time( );
-    data->timing.elapsed = Get_Timing_Info( data->timing.start );
-    if ( output_enabled == TRUE )
-    {
-        fprintf( out_control.log, "total: %.2f secs\n", data->timing.elapsed );
+        ret = SPUREMD_SUCCESS;
     }
 
-    return SUCCESS;
+    return ret;
 }
 
 
-int Cleanup( reax_system * const system, control_params * const control,
-        simulation_data * const data, const int output_enabled )
+int cleanup( const void * const handle )
 {
-    Finalize( system, control, data, &workspace, &lists, &out_control,
-           output_enabled );
+    int ret;
+    spuremd_handle *spmd_handle;
+
+    ret = SPUREMD_FAILURE;
 
-    sfree( lists, "main::lists" );
+    if ( handle != NULL )
+    {
+        spmd_handle = (spuremd_handle*) 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 );
+
+        sfree( spmd_handle->out_control, "cleanup::spmd_handle->out_control" );
+        sfree( spmd_handle->lists, "cleanup::spmd_handle->lists" );
+        sfree( spmd_handle->workspace, "cleanup::spmd_handle->workspace" );
+        sfree( spmd_handle->data, "cleanup::spmd_handle->data" );
+        sfree( spmd_handle->control, "cleanup::spmd_handle->control" );
+        sfree( spmd_handle->system, "cleanup::spmd_handle->system" );
+
+        sfree( spmd_handle, "cleanup::spmd_handle" );
+
+        ret = SPUREMD_SUCCESS;
+    }
 
-    return SUCCESS;
+    return ret;
 }
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index 7cad1a69..7ed56288 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -22,17 +22,16 @@
 #ifndef __SPUREMD_H_
 #define __SPUREMD_H_
 
-#include "mytypes.h"
+#define SPUREMD_SUCCESS (0)
+#define SPUREMD_FAILURE (-1)
 
 
-int Setup( char **, reax_system * const, control_params * const,
-        simulation_data * const );
+void* setup( const char * const, const char * const,
+        const char * const );
 
-int Run( reax_system * const, control_params * const,
-        simulation_data * const, const int );
+int simulate( const void * const );
 
-int Cleanup( reax_system * const, control_params * const,
-        simulation_data * const, const int );
+int cleanup( const void * const );
 
 
 #endif
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 6b59d194..3853c270 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -561,6 +561,4 @@ void sfree( void *ptr, const char *name )
 #endif
 
     free( ptr );
-
-    ptr = NULL;
 }
diff --git a/tools/driver.py b/tools/driver.py
new file mode 100644
index 00000000..0144e62c
--- /dev/null
+++ b/tools/driver.py
@@ -0,0 +1,142 @@
+#!/bin/python3
+
+from ctypes import *
+
+
+class Thermostat(Structure):
+    _fields_ = [
+            ("T", c_double),
+            ("xi", c_double),
+            ("v_xi", c_double),
+            ("v_xi_old", c_double),
+            ("G_xi", c_double),
+            ]
+
+
+class IsotropicBarostat(Structure):
+    _fields_ = [
+            ("P", c_double),
+            ("eps", c_double),
+            ("v_eps", c_double),
+            ("v_eps_old", c_double),
+            ("a_eps", c_double),
+            ]
+
+
+class FlexibleBarostat(Structure):
+    _fields_ = [
+            ("P", c_double * 9),
+            ("P_scalar", c_double),
+            ("eps", c_double),
+            ("v_eps", c_double),
+            ("v_eps_old", c_double),
+            ("a_eps", c_double),
+            ("h0", c_double * 9),
+            ("v_g0", c_double * 9),
+            ("v_g0_old", c_double * 9),
+            ("a_g0", c_double * 9),
+            ]
+
+
+class ReaxTiming(Structure):
+    _fields_ = [
+            ("start", c_double),
+            ("end", c_double),
+            ("elapsed", c_double),
+            ("total", c_double),
+            ("nbrs", c_double),
+            ("init_forces", c_double),
+            ("bonded", c_double),
+            ("nonb", c_double),
+            ("cm", c_double),
+            ("cm_sort_mat_rows", c_double),
+            ("cm_solver_pre_comp", c_double),
+            ("cm_solver_pre_app", c_double),
+            ("cm_solver_iters", c_int),
+            ("cm_solver_spmv", c_double),
+            ("cm_solver_vector_ops", c_double),
+            ("cm_solver_orthog", c_double),
+            ("cm_solver_tri_solve", c_double),
+            ]
+
+
+class SimulationData(Structure):
+    _fields_ = [
+            ("step", c_int),
+            ("prev_step", c_int),
+            ("time", c_double),
+            ("M", c_double),
+            ("inv_M", c_double),
+            ("xcm", c_double * 3),
+            ("vcm", c_double * 3),
+            ("fcm", c_double * 3),
+            ("amcm", c_double * 3),
+            ("avcm", c_double * 3),
+            ("etran_cm", c_double),
+            ("erot_cm", c_double),
+            ("kinetic", c_double * 9),
+            ("virial", c_double * 9),
+            ("E_Tot", c_double),
+            ("E_Kin", c_double),
+            ("E_Pot", c_double),
+            ("E_BE", c_double),
+            ("E_Ov", c_double),
+            ("E_Un", c_double),
+            ("E_Lp", c_double),
+            ("E_Ang", c_double),
+            ("E_Pen", c_double),
+            ("E_Coa", c_double),
+            ("E_HB", c_double),
+            ("E_Tor", c_double),
+            ("E_Con", c_double),
+            ("E_vdW", c_double),
+            ("E_Ele", c_double),
+            ("E_Pol", c_double),
+            ("N_f", c_double),
+            ("t_scale", c_double * 3),
+            ("p_scale", c_double * 9),
+            ("therm", Thermostat),
+            ("iso_bar", IsotropicBarostat),
+            ("flex_bar", FlexibleBarostat),
+            ("inv_W", c_double),
+            ("int_press", c_double * 3),
+            ("ext_press", c_double * 3),
+            ("kin_press", c_double),
+            ("tot_press", c_double * 3),
+            ("tot_press", c_double * 3),
+            ("timing", ReaxTiming),
+            ]
+
+
+class ReaxAtom(Structure):
+    _fields_ = [
+            ("name", c_char * 8),
+            ("x", c_double * 3),
+            ("v", c_double * 3),
+            ("f", c_double * 3),
+            ("q", c_double),
+            ]
+
+
+if __name__ == '__main__':
+    lib = cdll.LoadLibrary("libspuremd.so.1")
+
+    setup = lib.setup
+    setup.argtypes = [c_char_p, c_char_p, c_char_p]
+    setup.restype = c_void_p
+
+    simulate = lib.simulate
+    simulate.argtypes = [c_void_p]
+    simulate.restype = c_int
+
+    cleanup = lib.cleanup
+    cleanup.argtypes = [c_void_p]
+    cleanup.restype = c_int
+
+    handle = setup(b"data/benchmarks/water/water_6540.pdb",
+            b"data/benchmarks/water/ffield.water",
+            b"environ/param.gpu.water")
+
+    simulate(handle)
+
+    cleanup(handle)
-- 
GitLab