compile

-#!/bin/sh
+#! /bin/sh
 # Wrapper for compilers which do not understand '-c -o'.
 
 scriptversion=2016-01-11.22; # UTC

ltmain.sh

@@ -7220,6 +7220,11 @@ func_mode_link ()
 	arg=$func_stripname_result
 	;;
 
+      -Wl,--as-needed)
+	deplibs="$deplibs $arg"
+	continue
+	;;
+
       -Wl,*)
 	func_stripname '-Wl,' '' "$arg"
 	args=$func_stripname_result
@@ -7272,12 +7277,10 @@ func_mode_link ()
       # -tp=*                Portland pgcc target processor selection
       # --sysroot=*          for sysroot support
       # -O*, -g*, -flto*, -fwhopr*, -fuse-linker-plugin GCC link-time optimization
-      # -specs=*             GCC specs files
       # -stdlib=*            select c++ std lib with clang
       -64|-mips[0-9]|-r[0-9][0-9]*|-xarch=*|-xtarget=*|+DA*|+DD*|-q*|-m*| \
       -t[45]*|-txscale*|-p|-pg|--coverage|-fprofile-*|-F*|@*|-tp=*|--sysroot=*| \
-      -O*|-g*|-flto*|-fwhopr*|-fuse-linker-plugin|-fstack-protector*|-stdlib=*| \
-      -specs=*)
+      -O*|-g*|-flto*|-fwhopr*|-fuse-linker-plugin|-fstack-protector*|-stdlib=*)
         func_quote_for_eval "$arg"
 	arg=$func_quote_for_eval_result
         func_append compile_command " $arg"
@@ -7526,6 +7529,7 @@ func_mode_link ()
 
     case $linkmode in
     lib)
+	as_needed_flag=
 	passes="conv dlpreopen link"
 	for file in $dlfiles $dlprefiles; do
 	  case $file in
@@ -7537,6 +7541,7 @@ func_mode_link ()
 	done
 	;;
     prog)
+	as_needed_flag=
 	compile_deplibs=
 	finalize_deplibs=
 	alldeplibs=false
@@ -7606,6 +7611,15 @@ func_mode_link ()
 	lib=
 	found=false
 	case $deplib in
+	-Wl,--as-needed)
+	  if test prog,link = "$linkmode,$pass" ||
+	     test lib,link = "$linkmode,$pass"; then
+	    as_needed_flag="$deplib "
+	  else
+	    deplibs="$deplib $deplibs"
+	  fi
+	  continue
+	  ;;
 	-mt|-mthreads|-kthread|-Kthread|-pthread|-pthreads|--thread-safe \
         |-threads|-fopenmp|-openmp|-mp|-xopenmp|-omp|-qsmp=*)
 	  if test prog,link = "$linkmode,$pass"; then
@@ -10023,6 +10037,13 @@ EOF
 	  test "X$libobjs" = "X " && libobjs=
 	fi
 
+	# A bit hacky. I had wanted to add \$as_needed_flag to archive_cmds instead, but that
+	# comes from libtool.m4 which is part of the project being built. This should put it
+	# in the right place though.
+	if test lib,link = "$linkmode,$pass" && test -n "$as_needed_flag"; then
+	  libobjs=$as_needed_flag$libobjs
+	fi
+
 	save_ifs=$IFS; IFS='~'
 	for cmd in $cmds; do
 	  IFS=$sp$nl
@@ -10255,8 +10276,8 @@ EOF
       compile_deplibs=$new_libs
 
 
-      func_append compile_command " $compile_deplibs"
-      func_append finalize_command " $finalize_deplibs"
+      func_append compile_command " $as_needed_flag $compile_deplibs"
+      func_append finalize_command " $as_needed_flag $finalize_deplibs"
 
       if test -n "$rpath$xrpath"; then
 	# If the user specified any rpath flags, then add them.

m4/libtool.m4

@@ -2867,9 +2867,6 @@ linux* | k*bsd*-gnu | kopensolaris*-gnu | gnu*)
   # before this can be enabled.
   hardcode_into_libs=yes
 
-  # Add ABI-specific directories to the system library path.
-  sys_lib_dlsearch_path_spec="/lib64 /usr/lib64 /lib /usr/lib"
-
   # Ideally, we could use ldconfig to report *all* directores which are
   # searched for libraries, however this is still not possible.  Aside from not
   # being certain /sbin/ldconfig is available, command
@@ -2878,7 +2875,7 @@ linux* | k*bsd*-gnu | kopensolaris*-gnu | gnu*)
   # appending ld.so.conf contents (and includes) to the search path.
   if test -f /etc/ld.so.conf; then
     lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s 2>/dev/null", \[$]2)); skip = 1; } { if (!skip) print \[$]0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;/^[	 ]*hwcap[	 ]/d;s/[:,	]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;s/"//g;/^$/d' | tr '\n' ' '`
-    sys_lib_dlsearch_path_spec="$sys_lib_dlsearch_path_spec $lt_ld_extra"
+    sys_lib_dlsearch_path_spec="/lib /usr/lib $lt_ld_extra"
   fi
 
   # We used to test for /lib/ld.so.1 and disable shared libraries on

missing

-#!/bin/sh
+#! /bin/sh
 # Common wrapper for a few potentially missing GNU programs.
 
 scriptversion=2016-01-11.22; # UTC

sPuReMD/src/geo_tools.c

@@ -487,8 +487,13 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
     char element[6], charge[9];
     char chain_id;
     char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
-    char *endptr = NULL;
-    int i, atom_cnt, token_cnt, bgf_serial, ratom = 0;
+    char *endptr;
+    int i, atom_cnt, token_cnt, bgf_serial, ratom, crystx_found;
+
+    endptr = NULL;
+    ratom = 0;
+    crystx_found = FALSE;
+
 
     /* open biograf file */
     bgf = sfopen( bgf_file, "r" );
@@ -512,14 +517,16 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
         tokens[0][0] = 0;
         token_cnt = Tokenize( line, &tokens, MAX_TOKEN_LEN );
 
-        if ( !strncmp( tokens[0], "ATOM", MAX_TOKEN_LEN )
-                || !strncmp( tokens[0], "HETATM", MAX_TOKEN_LEN ) )
+        if ( !strncmp( tokens[0], "ATOM", 4 )
+                || !strncmp( tokens[0], "HETATM", 6 ) )
+
         {
             ++system->N;
         }
 
         line[0] = 0;
     }
+
     if ( ferror( bgf ) )
     {
         fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
@@ -664,6 +671,285 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
             Setup_Box( atof(s_a), atof(s_b), atof(s_c),
                     atof(s_alpha), atof(s_beta), atof(s_gamma),
                     &system->box );
+            crystx_found = TRUE;
+
+        }
+        else if ( !strncmp( tokens[0], "CONECT", 6 ) )
+        {
+            /* check number of restrictions */
+            Check_Input_Range( token_cnt - 2, 0, MAX_RESTRICT,
+                    "CONECT line exceeds max restrictions allowed.\n" );
+
+            /* read bond restrictions */
+            if ( is_Valid_Serial( workspace, bgf_serial = atoi(tokens[1]) ) )
+            {
+                ratom = workspace->map_serials[ bgf_serial ];
+            }
+
+            workspace->restricted[ ratom ] = token_cnt - 2;
+            for ( i = 2; i < token_cnt; ++i )
+            {
+                if ( is_Valid_Serial( workspace, bgf_serial = atoi(tokens[i]) ) )
+                {
+                    workspace->restricted_list[ratom][i - 2] =
+                        workspace->map_serials[bgf_serial];
+                }
+            }
+
+            /* fprintf( stderr, "restriction on %d:", ratom );
+            for( i = 0; i < workspace->restricted[ ratom ]; ++i )
+             fprintf( stderr, "  %d", workspace->restricted_list[ratom][i] );
+             fprintf( stderr, "\n" ); */
+        }
+
+        /* clear previous input line */
+        line[0] = 0;
+
+        for ( i = 0; i < token_cnt; ++i )
+        {
+            tokens[i][0] = 0;
+        }
+    }
+    if ( ferror( bgf ) )
+    {
+        fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
+        exit( INVALID_INPUT );
+    }
+    
+    if ( crystx_found == FALSE )
+    {
+        fprintf( stderr, "[ERROR] improperly formatted BGF file (no CRYSTX keyword found). Terminating...\n" );
+        exit( INVALID_INPUT );
+    }
+
+    sfree( line, "Read_BGF::line" );
+    sfree( backup, "Read_BGF::backup" );
+    for ( i = 0; i < MAX_TOKENS; i++ )
+    {
+        sfree( tokens[i], "Read_BGF::tokens[i]" );
+    }
+    sfree( tokens, "Read_BGF::tokens" );
+
+
+    sfclose( bgf, "Read_BGF::bgf" );
+}
+
+char * sgets(char *s, int n, const char **strp){
+    if(**strp == '\0')return NULL;
+    int i;
+    for(i=0;i<n-1;++i, ++(*strp)){
+        s[i] = **strp;
+        if(**strp == '\0')
+            break;
+        if(**strp == '\n'){
+            s[i+1]='\0';
+            ++(*strp);
+            break;
+        }
+    }
+    if(i==n-1)
+        s[i] = '\0';
+    return s;
+}
+
+void Read_BGF_from_string( char * bgf_str, reax_system* system, control_params *control,
+        simulation_data *data, static_storage *workspace )
+{
+    char **tokens;
+    char *line, *backup;
+    char descriptor[7], serial[6];
+    char atom_name[6], res_name[4], res_seq[6];
+    char s_x[11], s_y[11], s_z[11];
+    char occupancy[4], temp_factor[3];
+    char element[6], charge[9];
+    char chain_id;
+    char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
+    char *endptr;
+    int i, atom_cnt, token_cnt, bgf_serial, ratom, crystx_found;
+
+    endptr = NULL;
+    ratom = 0;
+    crystx_found = FALSE;
+
+
+
+    /* allocate memory for tokenizing biograf file lines */
+    line = smalloc( sizeof(char) * MAX_LINE, "Read_BGF::line" );
+    backup = smalloc( sizeof(char) * MAX_LINE, "Read_BGF::backup" );
+    tokens = smalloc( sizeof(char*) * MAX_TOKENS, "Read_BGF::tokens" );
+    for ( i = 0; i < MAX_TOKENS; i++ )
+    {
+        tokens[i] = smalloc( sizeof(char) * MAX_TOKEN_LEN,
+               "Read_BGF::tokens[i]" );
+    }
+
+    /* count number of atoms in the pdb file */
+    system->N = 0;
+    line[0] = 0;
+    const char** new_bgf_ptr = &bgf_str;
+    while(sgets(line, MAX_LINE, new_bgf_ptr))
+    {
+        tokens[0][0] = 0;
+        token_cnt = Tokenize( line, &tokens, MAX_TOKEN_LEN );
+
+        if ( !strncmp( tokens[0], "ATOM", 4 )
+                || !strncmp( tokens[0], "HETATM", 6 ) )
+
+        {
+            ++system->N;
+        }
+
+        line[0] = 0;
+    }
+    printf("Number of atoms: %d\n", system->N);
+    /*
+    if ( ferror( bgf ) )
+    {
+        fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
+        exit( INVALID_INPUT );
+    }
+    
+    sfclose( bgf, "Read_BGF::bgf" );
+    */
+    /* memory allocations for atoms, atom maps, bond restrictions */
+    system->atoms = scalloc( system->N, sizeof(reax_atom),
+            "Read_BGF::system->atoms" );
+
+    workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
+            "Read_BGF::workspace->map_serials" );
+    for ( i = 0; i < MAX_ATOM_ID; ++i )
+    {
+        workspace->map_serials[i] = -1;
+    }
+
+    workspace->orig_id = scalloc( system->N, sizeof(int),
+            "Read_BGF::workspace->orig_id" );
+    workspace->restricted  = scalloc( system->N, sizeof(int),
+            "Read_BGF::workspace->restricted" );
+    workspace->restricted_list = scalloc( system->N, sizeof(int*),
+            "Read_BGF::workspace->restricted_list" );
+    for ( i = 0; i < system->N; ++i )
+    {
+        workspace->restricted_list[i] = scalloc( MAX_RESTRICT, sizeof(int),
+                "Read_BGF::workspace->restricted_list[i]" );
+    }
+
+    /* start reading and processing bgf file 
+    bgf = sfopen( bgf_file, "r" );
+    */
+    atom_cnt = 0;
+    token_cnt = 0;
+    const char** new_bgf_ptr_2 = &bgf_str;
+    while(sgets(line, MAX_LINE, new_bgf_ptr_2))
+    {
+        /* read new line and tokenize it */
+        strncpy( backup, line, MAX_LINE - 1 );
+        backup[MAX_LINE - 1] = '\0';
+        token_cnt = Tokenize( line, &tokens, MAX_TOKEN_LEN );
+
+        /* process new line */
+        if ( !strncmp(tokens[0], "ATOM", 4) || !strncmp(tokens[0], "HETATM", 6) )
+        {
+            if ( !strncmp(tokens[0], "ATOM", 4) )
+            {
+                strncpy( descriptor, backup, sizeof(descriptor) - 1 );
+                descriptor[sizeof(descriptor) - 1] = '\0';
+                strncpy( serial, backup + 7, sizeof(serial) - 1 );
+                serial[sizeof(serial) - 1] = '\0';
+                strncpy( atom_name, backup + 13, sizeof(atom_name) - 1 );
+                atom_name[sizeof(atom_name) - 1] = '\0';
+                strncpy( res_name, backup + 19, sizeof(res_name) - 1 );
+                res_name[sizeof(res_name) - 1] = '\0';
+                chain_id = backup[23];
+                strncpy( res_seq, backup + 25, sizeof(res_seq) - 1 );
+                res_seq[sizeof(res_seq) - 1] = '\0';
+                strncpy( s_x, backup + 30, sizeof(s_x) - 1 );
+                s_x[sizeof(s_x) - 1] = '\0';
+                strncpy( s_y, backup + 40, sizeof(s_y) - 1 );
+                s_y[sizeof(s_y) - 1] = '\0';
+                strncpy( s_z, backup + 50, sizeof(s_x) - 1 );
+                s_z[sizeof(s_x) - 1] = '\0';
+                strncpy( element, backup + 61, sizeof(element) - 1 );
+                element[sizeof(element) - 1] = '\0';
+                strncpy( occupancy, backup + 66, sizeof(occupancy) - 1 );
+                occupancy[sizeof(occupancy) - 1] = '\0';
+                strncpy( temp_factor, backup + 69, sizeof(temp_factor) - 1 );
+                temp_factor[sizeof(temp_factor) - 1] = '\0';
+                strncpy( charge, backup + 72, sizeof(charge) - 1 );
+                charge[sizeof(charge) - 1] = '\0';
+            }
+            else if ( !strncmp(tokens[0], "HETATM", 6) )
+            {
+                /* bgf hetatm:
+                   (7x,i5,1x,a5,1x,a3,1x,a1,1x,a5,3f10.5,1x,a5,i3,i2,1x,f8.5) */
+                strncpy( descriptor, backup, sizeof(descriptor) - 1 );
+                descriptor[sizeof(descriptor) - 1] = '\0';
+                strncpy( serial, backup + 7, sizeof(serial) - 1 );
+                serial[sizeof(serial) - 1] = '\0';
+                strncpy( atom_name, backup + 13, sizeof(atom_name) - 1 );
+                atom_name[sizeof(atom_name) - 1] = '\0';
+                strncpy( res_name, backup + 19, sizeof(res_name) - 1 );
+                res_name[sizeof(res_name) - 1] = '\0';
+                chain_id = backup[23];
+                strncpy( res_seq, backup + 25, sizeof(res_seq) - 1 );
+                res_seq[sizeof(res_seq) - 1] = '\0';
+                strncpy( s_x, backup + 30, sizeof(s_x) - 1 );
+                s_x[sizeof(s_x) - 1] = '\0';
+                strncpy( s_y, backup + 40, sizeof(s_y) - 1 );
+                s_y[sizeof(s_y) - 1] = '\0';
+                strncpy( s_z, backup + 50, sizeof(s_z) - 1 );
+                s_z[sizeof(s_z) - 1] = '\0';
+                strncpy( element, backup + 61, sizeof(element) - 1 );
+                element[sizeof(element) - 1] = '\0';
+                strncpy( occupancy, backup + 66, sizeof(occupancy) - 1 );
+                occupancy[sizeof(occupancy) - 1] = '\0';
+                strncpy( temp_factor, backup + 69, sizeof(temp_factor) - 1 );
+                temp_factor[sizeof(temp_factor) - 1] = '\0';
+                strncpy( charge, backup + 72, sizeof(charge) - 1 );
+                charge[sizeof(charge) - 1] = '\0';
+            }
+
+            /* add to mapping */
+            bgf_serial = strtod( serial, &endptr );
+            Check_Input_Range( bgf_serial, 0, MAX_ATOM_ID, "[ERROR] Invalid bgf serial" );
+            workspace->map_serials[ bgf_serial ] = atom_cnt;
+            workspace->orig_id[ atom_cnt ] = bgf_serial;
+            // fprintf( stderr, "map %d --> %d\n", bgf_serial, atom_cnt );
+
+            /* copy atomic positions */
+            system->atoms[atom_cnt].x[0] = strtod( s_x, &endptr );
+            system->atoms[atom_cnt].x[1] = strtod( s_y, &endptr );
+            system->atoms[atom_cnt].x[2] = strtod( s_z, &endptr );
+
+            /* atom name and type */
+            strncpy( system->atoms[atom_cnt].name, atom_name,
+                    sizeof(system->atoms[atom_cnt].name) - 1 );
+            system->atoms[atom_cnt].name[sizeof(system->atoms[atom_cnt].name) - 1] = '\0';
+            Trim_Spaces( element, 6 );
+            system->atoms[atom_cnt].type =
+                Get_Atom_Type( &system->reax_param, element );
+
+            /* fprintf( stderr,
+            "a:%3d(%1d) c:%10.5f%10.5f%10.5f q:%10.5f occ:%s temp:%s seg_id:%s element:%s\n",
+             atom_cnt, system->atoms[ atom_cnt ].type,
+             system->atoms[ atom_cnt ].x[0],
+             system->atoms[ atom_cnt ].x[1], system->atoms[ atom_cnt ].x[2],
+             system->atoms[ atom_cnt ].q, occupancy, temp_factor,
+             seg_id, element ); */
+
+            atom_cnt++;
+        }
+        else if ( !strncmp( tokens[0], "CRYSTX", 6 ) )
+        {
+            sscanf( backup, BGF_CRYSTX_FORMAT, descriptor,
+                    s_a, s_b, s_c, s_alpha, s_beta, s_gamma );
+
+            /* Compute full volume tensor from the angles */
+            Setup_Box( atof(s_a), atof(s_b), atof(s_c),
+                    atof(s_alpha), atof(s_beta), atof(s_gamma),
+                    &system->box );
+                    
+            crystx_found = TRUE;
         }
         else if ( !strncmp( tokens[0], "CONECT", 6 ) )
         {
@@ -701,11 +987,29 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
             tokens[i][0] = 0;
         }
     }
+    
+    if ( crystx_found == FALSE )
+    {
+        fprintf( stderr, "[ERROR] improperly formatted BGF file (no CRYSTX keyword found). Terminating...\n" );
+        exit( INVALID_INPUT );
+    }
+
+    sfree( line, "Read_BGF::line" );
+    sfree( backup, "Read_BGF::backup" );
+    for ( i = 0; i < MAX_TOKENS; i++ )
+    {
+        sfree( tokens[i], "Read_BGF::tokens[i]" );
+    }
+    sfree( tokens, "Read_BGF::tokens" );
+
+    /*
     if ( ferror( bgf ) )
     {
         fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
         exit( INVALID_INPUT );
     }
 
+
     sfclose( bgf, "Read_BGF::bgf" );
+    */
 }

sPuReMD/src/geo_tools.h

@@ -125,6 +125,9 @@ void Read_PDB( const char * const, reax_system*, control_params*,
 void Read_BGF( const char * const, reax_system*, control_params*,
         simulation_data*, static_storage* );
 
+void Read_BGF_from_string( char *, reax_system*, control_params*,
+        simulation_data*, static_storage* );
+
 void Write_PDB( reax_system*, reax_list*, simulation_data*,
         control_params*, static_storage*, output_controls* );
 

sPuReMD/src/init_md.c

@@ -766,11 +766,6 @@ static void Init_Lists( reax_system *system, control_params *control,
      *   (non-bonded, hydrogen, 3body, etc.) */
     Allocate_Matrix( &workspace->H_sp, system->N_cm, max_nnz );
 
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "estimated storage - Htop: %d\n", Htop );
-    fprintf( stderr, "memory allocated: H = %ldMB\n",
-            Htop * sizeof(sparse_matrix_entry) / (1024 * 1024) );
-#endif
 
     workspace->num_H = 0;
     if ( control->hbond_cut > 0.0 )

sPuReMD/src/nonbonded.c

@@ -282,7 +282,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
                         }
 
 #if defined(DEBUG_FOCUS)
-                        fprintf( stderr, "nonbonded(%d,%d): rel_box (%f %f %f)",
+                        fprintf( stderr, "nonbonded(%d,%d): rel_box (%d %d %d)",
                                 i, j, nbr_pj->rel_box[0],
                                 nbr_pj->rel_box[1], nbr_pj->rel_box[2] );
 

sPuReMD/src/spuremd.c

@@ -22,6 +22,10 @@
 #include "spuremd.h"
 
 #include "analyze.h"
+#if defined(DEBUG_FOCUS)
+  #include "box.h"
+#endif
+
 #include "control.h"
 #include "ffield.h"
 #include "forces.h"
@@ -73,7 +77,6 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
     }
 }
 
-
 static void Read_System( const char * const geo_file,
         const char * const ffield_file,
         const char * const control_file,
@@ -128,10 +131,40 @@ static void Read_System( const char * const geo_file,
 
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "input files have been read...\n" );
-    Print_Box( &(system->box), stderr );
+    Print_Box( &system->box, stderr );
 #endif
 }
 
+static void Read_System2(char * geo_str,
+        const char * const ffield_file,
+        const char * const control_file,
+        reax_system * const system,
+        control_params * const control,
+        simulation_data * const data,
+        static_storage * const workspace,
+        output_controls * const out_control )
+{
+    FILE *ffield, *ctrl;
+
+    ffield = sfopen( ffield_file, "r" );
+    ctrl = sfopen( control_file, "r" );
+
+    /* ffield file */
+    Read_Force_Field( ffield, &system->reax_param );
+
+    /* control file */
+    Read_Control_File( ctrl, system, control, out_control );
+
+    Read_BGF_from_string( geo_str, system, control, data, workspace );
+
+    sfclose( ffield, "Read_System::ffield" );
+    sfclose( ctrl, "Read_System::ctrl" );
+
+#if defined(DEBUG_FOCUS)
+    fprintf( stderr, "input files have been read...\n" );
+    Print_Box( &(system->box), stderr );
+#endif
+}
 
 void* setup( const char * const geo_file, const char * const ffield_file,
         const char * const control_file )
@@ -175,6 +208,48 @@ void* setup( const char * const geo_file, const char * const ffield_file,
     return (void*) spmd_handle;
 }
 
+void* setup2( char * geo_str, const char * const ffield_file,
+        const char * const control_file )
+{
+    int i;
+    spuremd_handle *spmd_handle;
+
+    /* top-level allocation */
+    spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
+            "setup::spmd_handle" );
+
+    /* second-level allocations */
+    spmd_handle->system = smalloc( sizeof(reax_system),
+           "Setup::spmd_handle->system" );
+    spmd_handle->control = smalloc( sizeof(control_params),
+           "Setup::spmd_handle->control" );
+    spmd_handle->data = smalloc( sizeof(simulation_data),
+           "Setup::spmd_handle->data" );
+    spmd_handle->workspace = smalloc( sizeof(static_storage),
+           "Setup::spmd_handle->workspace" );
+    spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
+           "Setup::spmd_handle->lists" );
+    for ( i = 0; i < LIST_N; ++i )
+    {
+        spmd_handle->lists[i] = smalloc( sizeof(reax_list),
+                "Setup::spmd_handle->lists[i]" );
+//        spmd_handle->lists[i]->allocated = FALSE;
+    }
+    spmd_handle->out_control = smalloc( sizeof(output_controls),
+           "Setup::spmd_handle->out_control" );
+
+    spmd_handle->output_enabled = TRUE;
+    spmd_handle->callback = NULL;
+
+    /* parse geometry file */
+    Read_System2( geo_str, ffield_file, control_file,
+            spmd_handle->system, spmd_handle->control,
+            spmd_handle->data, spmd_handle->workspace,
+            spmd_handle->out_control );
+
+    return (void*) spmd_handle;
+}
+
 
 int setup_callback( const void * const handle, const callback_function callback  )
 {

sPuReMD/src/spuremd.h

@@ -27,7 +27,7 @@
 
 #define SPUREMD_SUCCESS (0)
 #define SPUREMD_FAILURE (-1)
-
+#define DEBUG_FOCUS (1)
 
 #ifdef __cplusplus
 extern "C"  {
@@ -36,6 +36,9 @@ extern "C"  {
 void* setup( const char * const, const char * const,
         const char * const );
 
+void* setup2( char *, const char * const,
+        const char * const );
+
 int setup_callback( const void * const, const callback_function );
 
 int simulate( const void * const );

tools/chosen_run

+BIOGRF 200
+DESCRP si3                                                         
+REMARK                                                             
+RUTYPE NORMAL RUN
+#              At1 At2   R12    Force1  Force2  dR12/dIter(MD) Start (MD) End (MD)
+BOND RESTRAINT    1   2  2.1100 2500.00  1.0000  0.0000000       0       0
+FORMAT ATOM   (a6,1x,i5,1x,a5,1x,a3,1x,a1,1x,a5,3f10.5,1x,a5,i3,i2,1x,f8.5)
+HETATM     1 Si                 40.00024  39.99963  41.06185    Si  1 1  0.21301
+HETATM     2 Si                 39.99975  40.00040  38.93811    Si  1 1  0.21301
+HETATM     3 H                  41.30813  40.00130  38.47248    H   1 1 -0.07100
+HETATM     4 H                  38.69272  39.99889  41.52902    H   1 1 -0.07100
+HETATM     5 H                  39.34629  41.13476  38.47406    H   1 1 -0.07100
+HETATM     6 H                  39.34396  38.86680  38.47073    H   1 1 -0.07100
+HETATM     7 H                  40.65381  41.13161  41.52612    H   1 1 -0.07100
+HETATM     8 H                  40.65560  38.86582  41.52858    H   1 1 -0.07100
+FORMAT CONECT (a6,12i6)
+CONECT     1     2     4     7     8
+CONECT     2     1     3     5     6
+CONECT     3     2
+CONECT     4     1
+CONECT     5     2
+CONECT     6     2
+CONECT     7     1
+CONECT     8     1
+UNIT ENERGY   kcal
+ENERGY        -530.795812
\ No newline at end of file

tools/driver.py

@@ -3,8 +3,10 @@
 from ctypes import c_int, c_double, c_char, c_char_p, c_void_p, \
         Structure, Union, POINTER, CFUNCTYPE, cdll
 import sqlite3 as sq3
+import os
 from os import path
-
+import read_geo
+import read_train
 
 class BondOrderData(Structure):
     _fields_ = [
@@ -248,6 +250,11 @@ class ReaxAtom(Structure):
             ("q", c_double),
             ]
 
+def print_file(filename, my_str):
+    fptr = open(filename, 'w')
+    fptr.write(my_str)
+    fptr.close()
+
 
 def create_db(name='spuremd.db'):
     conn = sq3.connect(name)
@@ -385,9 +392,11 @@ def create_db(name='spuremd.db'):
 
 
 if __name__ == '__main__':
-    lib = cdll.LoadLibrary("libspuremd.so.1")
+    lib = cdll.LoadLibrary("/mnt/home/kaymakme/SensitivityAnalysis/PuReMD/tools/libspuremd.so.1.0.0")
 
     setup = lib.setup
+    setup2 = lib.setup2
+
     setup.argtypes = [c_char_p, c_char_p, c_char_p]
     setup.restype = c_void_p
 
@@ -413,13 +422,14 @@ if __name__ == '__main__':
     set_output_enabled = lib.set_output_enabled
     set_output_enabled.argtypes = [c_void_p, c_int]
     set_output_enabled.restype = c_int
-
+    '''
     db_file = "spuremd.db"
 
     if not path.isfile(db_file):
         create_db(db_file)
 
     conn = sq3.connect(db_file)
+    '''
     record_potential = True
     record_trajectory = False
     record_performance = True
@@ -431,7 +441,7 @@ if __name__ == '__main__':
         if data[0].step == 0:
             #TODO: insert data into simulation table
             pass
-
+        '''
         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,
@@ -463,19 +473,29 @@ if __name__ == '__main__':
                             data[0].timing.cm_solver_pre_app, data[0].timing.cm_solver_spmv,
                             data[0].timing.cm_solver_vector_ops, data[0].timing.cm_solver_orthog,
                             data[0].timing.cm_solver_tri_solve))
-
-
-    handle = setup(b"data/benchmarks/water/water_6540.pdb",
-            b"data/benchmarks/water/ffield.water",
-            b"environ/control_water")
-
+        '''
+    geo_file = '/mnt/home/kaymakme/SensitivityAnalysis/SiOH/geo'
+    runs = read_geo.parse_geo_file(geo_file)
+    print(runs[5][:-1])
+    filename = '/mnt/home/kaymakme/SensitivityAnalysis/PuReMD/tools/chosen_run'
+    #print_file(filename, runs[2])
+    # 
+    handle = setup(b'/mnt/home/kaymakme/SensitivityAnalysis/PuReMD/data/water-systems/water-6/water_6540.pdb',
+            b"/mnt/home/kaymakme/SensitivityAnalysis/SiOH/ffield",
+            b"/mnt/home/kaymakme/SensitivityAnalysis/PuReMD/environ/control_water")
+    '''
+    handle = setup2(runs[5][:-1].encode(),
+        b"/mnt/home/kaymakme/SensitivityAnalysis/SiOH/ffield",
+        b"/mnt/home/kaymakme/SensitivityAnalysis/PuReMD/environ/control_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)
-
+    print('aa')
     atoms = get_atoms(handle)
 
     print()
@@ -484,5 +504,5 @@ if __name__ == '__main__':
         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))
 
-    conn.close()
+    #conn.close()
     cleanup(handle)

tools/libspuremd.so.1

+libspuremd.so.1.0.0
\ No newline at end of file

tools/read_geo.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Jun 18 16:29:59 2019
+
+@author: cagri
+"""
+import  os
+
+def parse_geo_file(geo_file):
+    if not os.path.exists(geo_file):
+        return
+    list_subruns = []
+    f = open(geo_file,'r')
+    run_str = ''
+    for line in f:
+        if len(line.strip()) < 1:
+            list_subruns.append(run_str)
+            run_str = ''
+        else:
+            run_str = run_str + line
+    return list_subruns
+            

tools/read_train.py

+#!/usr/bin/env python
+import  os
+
+'''
+Simple parser for training data
+TODO: right now, it is strict about the spaces, it could be solved with
+a tokenizer: tokenize the input then parse the tokens
+'''
+def clear_comments(input_str):
+    '''
+    clear the comments from a line (only # type of comments for now)
+    example:
+        input:   chexane  0.01   1              0.0    # Eucledian distance
+        output:  chexane  0.01   1              0.0
+    '''
+    com_start_ind = input_str.find('#')
+    return input_str[:com_start_ind]
+
+def parse_trainset(trainset_file='trainset.in'):
+    parsed_cells = initial_parse_trainset(trainset_file)
+    return process_parsed_data(parsed_cells)
+    
+
+def initial_parse_trainset(trainset_file='trainset.in'):
+    '''
+    reads all the lines *between* the section indicators in
+    trainset.in and stores them in a cells dictionary.
+
+    this is so I can append new lines to existing sections, and later
+    write them back out. Each line is stripped.
+    '''
+
+    cells = {}
+    cells['CHARGE'] = []
+    cells['GEOMETRY'] = []
+    cells['FORCES'] = []
+    cells['CELL PARAMETERS'] = []
+    cells['ENERGY'] = []
+    cells['HEATFO'] = []
+    if not os.path.exists(trainset_file):
+        return cells
+
+    CHARGE_FLAG, GEO_FLAG, FORCE_FLAG, CELL_PARAM_FLAG, ENERGY_FLAG, \
+        HEAT_FLAG  = False, False, False, False, False, False
+        
+    f = open(trainset_file,'r')
+    for line in f:
+        line = clear_comments(line).strip()
+        if len(line) < 1:
+            continue        
+        
+        if line.startswith('CHARGE'):
+                CHARGE_FLAG = True
+                continue
+        elif line.startswith('ENDCHARGE'):
+                CHARGE_FLAG = False
+                
+        elif line.startswith('GEOMETRY'):
+                GEO_FLAG = True
+                continue
+        elif line.startswith('ENDGEOMETRY'):
+                GEO_FLAG = False  
+                
+        elif line.startswith('FORCES'):
+                FORCE_FLAG = True
+                continue
+        elif line.startswith('ENDFORCES'):
+                FORCE_FLAG = False
+          
+        elif line.startswith('CELL PARAMETERS'):
+                CELL_PARAM_FLAG = True
+                continue
+        elif line.startswith('ENDCELL PARAMETERS'):
+                CELL_PARAM_FLAG = False
+                
+        elif line.startswith('ENERGY'):
+                ENERGY_FLAG = True
+                continue
+        elif line.startswith('ENDENERGY'):
+                ENERGY_FLAG = False 
+                
+
+        elif line.startswith('HEATFO'):
+                HEAT_FLAG = True
+                continue
+        elif line.startswith('ENDHEATFO'):
+                HEAT_FLAG = False
+
+
+            
+        elif CHARGE_FLAG:
+            cells['CHARGE'].append(line)
+        elif GEO_FLAG:
+            cells['GEOMETRY'].append(line)        
+        elif FORCE_FLAG:
+            cells['FORCES'].append(line)
+        elif CELL_PARAM_FLAG:
+            cells['CELL PARAMETERS'].append(line)
+        elif ENERGY_FLAG:
+            cells['ENERGY'].append(line)
+        elif HEAT_FLAG:
+            cells['HEATFO'].append(line)
+        else:
+            print('The input file is not in the correct format!')
+    f.close()
+    return cells
+
+
+def process_parsed_data(cells):
+    '''
+     process the parsed data and create the required lists and dictionaries
+     corresponding to the parsed data and print the required error messages
+    '''
+    processed_data = {}
+    processed_data['CHARGE'] = []
+    processed_data['GEOMETRY'] = []
+    processed_data['FORCES'] = []
+    processed_data['CELL PARAMETERS'] = []
+    processed_data['ENERGY'] = []
+    processed_data['HEATFO'] = []
+    # process lines for charge
+    for line in cells['CHARGE']: 
+        # format: Key Acc Atom Ref
+        words = line.split()
+        if len(words) != 4:
+            print('something wrong with charge data')
+        charge_input = {}
+        charge_input['key'] = words[0]
+        charge_input['acc'] = words[1]
+        charge_input['atom'] = words[2]
+        charge_input['ref'] = words[3]
+        processed_data['CHARGE'].append(charge_input)
+        
+    # process lines for geo
+    for line in cells['GEOMETRY']: 
+        # format: Key Acc [Atom1 [Atom2 [Atom3 [Atom4]]] Ref]
+        words = line.split()
+        if len(words) < 3 or len(words) > 7:
+            print('something wrong with geometry data')
+        geo_input = {}
+        geo_input['key'] = words[0]
+        geo_input['acc'] = words[1]
+        if len(words) > 3: # can be multiple atoms or not at all
+            geo_input['atoms'] = words[2:-1] 
+        geo_input['ref'] = words[-1]
+        processed_data['GEOMETRY'].append(geo_input)
+
+    # process lines for forces
+    for line in cells['FORCES']: 
+        # format: Key Acc Type Ref
+        words = line.split()
+        if len(words) != 4:
+            print('something wrong with force data')
+        force_input = {}
+        force_input['key'] = words[0]
+        force_input['acc'] = words[1]
+        force_input['atom'] = words[2]
+        force_input['ref'] = words[3]
+        processed_data['FORCES'].append(force_input)
+        
+    # process lines for CELL PARAMETERS
+    for line in cells['CELL PARAMETERS']: 
+        # format: Key Acc Type Ref
+        words = line.split()
+        if len(words) != 4:
+            print('something wrong with CELL PARAMETERS data')
+        force_input = {}
+        force_input['key'] = words[0]
+        force_input['acc'] = words[1]
+        force_input['type'] = words[2]
+        force_input['ref'] = words[3]
+        processed_data['CELL PARAMETERS'].append(force_input)
+        
+    # process lines for ENERGY
+    for line in cells['ENERGY']: 
+        # format: Acc [+-] Key1/n1 ... [+-] Key5/n5 Ref
+        words = line.split()
+        '''
+        if len(words) > 2 and len(words) % 2 == 0:
+            print('something wrong with ENERGY data')
+        '''
+        energy_input = {}
+        energy_input['acc'] = words[0]
+        energy_input['parts'] = []
+        item = {}
+        for word in words[1:-1]:
+            if word in ['+','-']:
+                item = {}
+                item['op'] = word
+            elif word.find('/') > 0:
+                word_sub_parts = word.split('/')
+                item['key'] = word_sub_parts[0]
+                item['n']  = word_sub_parts[1]
+                energy_input['parts'].append(item)
+            # if there is a space like in "Key1 /n1"
+            else:
+                if word.find('/') == 0:
+                    item['n'] = word[1:]
+                    energy_input['parts'].append(item)
+                else:
+                    item['key'] = word
+                    
+                
+                
+        energy_input['ref'] = words[-1]
+        processed_data['ENERGY'].append(energy_input)
+        
+    # process lines for HEAT
+    for line in cells['HEATFO']: 
+        # format: Key Acc Type Ref
+        words = line.split()
+        if len(words) != 3:
+            print('something wrong with HEAT data')
+        heat_input = {}
+        heat_input['key'] = words[0]
+        heat_input['acc'] = words[1]
+        heat_input['ref'] = words[2]
+        processed_data['HEATFO'].append(heat_input)
+        
+    return processed_data
+        

tools/trainset.py

+#!/usr/bin/env python
+'''
+trainset.in contains input training data for reax.
+
+Here is the basic format of these sections.
+I am not sure what a Heat of formation in reax is. those appear to be captured in the ENERGY section.
+
+HEATFO
+descriptor weight reference_value
+ENDHEATFO
+
+CELL PARAMETERS
+#descriptor weight type reference value
+fcc0 1 a 4.0
+fcc0 1 b 4.0
+fcc0 1 c 4.0
+fcc0 1 alpha 90
+fcc0 1 beta 90
+fcc0 1 gamma 90
+# comment
+fccm1 1 a 3.9
+fccm1 1 b 3.9
+fccm1 1 c 3.9
+fccm1 1 alpha 90
+fccm1 1 beta 90
+fccm1 1 gamma 90
+#
+fccp1 1 a 4.1
+fccp1 1 b 4.1
+fccp1 1 c 4.1
+fccp1 1 alpha 90
+fccp1 1 beta 90
+fccp1 1 gamma 90
+END CELL PARAMETERS
+
+ENERGY
+weight + descriptor1 /a +- descriptor2 /b +- descriptor3 /c reference_value
+ENDENERGY
+'''
+from Cheetah.Template import Template
+from Scientific.IO.FortranFormat import FortranFormat, FortranLine
+from Scientific.Geometry import Vector
+import  math, pickle, os
+from reax import *
+
+def parse_trainset(trainset_file='trainset.in'):
+    '''
+    reads all the lines *between* the section indicators in
+    trainset.in and stores them in a cells dictionary.
+
+    this is so I can append new lines to existing sections, and later
+    write them back out. Each line is stripped.
+    '''
+
+    cells = {}
+    cells['ENERGY'] = []
+    cells['CELL PARAMETERS'] = []
+    cells['HEATFO'] = []
+    if not os.path.exists(trainset_file):
+        return cells
+
+    EFLAG, HFLAG, CFLAG = False, False, False
+    f = open(trainset_file,'r')
+    for line in f:
+        if line.startswith('ENERGY'):
+                EFLAG = True
+                continue
+        elif line.startswith('ENDENERGY'):
+                EFLAG = False
+
+        elif line.startswith('HEATFO'):
+                HFLAG = True
+                continue
+        elif line.startswith('ENDHEATFO'):
+                HFLAG = False
+
+        elif line.startswith('CELL PARAMETERS'):
+                CFLAG = True
+                continue
+        elif line.startswith('END CELL PARAMETERS'):
+                CFLAG = False
+
+        if EFLAG:
+            cells['ENERGY'].append(line.strip())
+        elif CFLAG:
+            cells['CELL PARAMETERS'].append(line.strip())
+        elif HFLAG:
+            cells['HEATFO'].append(line.strip())
+    f.close()
+    return cells
+
+def write_trainset(trainset_file='trainset.in', cells=None):
+    '''
+    writes all entries contained in the cells dictionary to the
+    trainset file.
+    '''
+
+    f = open(trainset_file, 'w')
+    for key in cells:
+        if cells[key] == []:
+            continue
+        f.write('%s\n' % key)
+        for line in cells[key]:
+            f.write(line+'\n')
+        f.write('END%s\n' % key)
+    f.close()
+
+#######################################################
+# functions to generate trainset entries
+
+def cell_parameters(atoms, weight=1.0, cells=None):
+    '''
+    adds entries to CELL PARAMETERS for the unit cell geometry of the
+    atoms object.
+    '''
+
+    description = reax_hash(atoms)
+
+    # get all the parameters needed for a cell parameters block
+    uc = atoms.get_cell()
+    a,b,c = [Vector(x) for x in uc]
+    A = a.length()
+    B = b.length()
+    C = c.length()
+
+    rad2deg = 360./(2*math.pi)
+    alpha = b.angle(c)*rad2deg
+    beta  = a.angle(c)*rad2deg
+    gamma = a.angle(b)*rad2deg
+
+    cells['CELL PARAMETERS'].append('# created by reax.trainset')
+    # this is an org-link to the geo file
+    #cells['CELL PARAMETERS'].append('# file:geo::%s' % description)
+    for key,val in [('a',A),
+                    ('b',B),
+                    ('c',C),
+                    ('alpha',alpha),
+                    ('beta', beta),
+                    ('gamma', gamma)]:
+        s = '%(description)s %(weight)f %(key)s %(val)s' % locals()
+        if s not in cells['CELL PARAMETERS']:
+            cells['CELL PARAMETERS'].append(s)
+
+    return cells
+
+def energy(atoms, reference_energies,
+           comment=None,
+           weight=1.0, cells=None):
+    '''generate lines for formation energies
+
+    atoms is an ase.Atoms object
+    reference_energies is a dictionary:
+    chemical symbol: (description, energy)
+
+    cells is a dictionary that you are appending these entries too.
+    '''
+    description = reax_hash(atoms)
+    # add org-link
+    #cells['ENERGY'].append('# file:geo::%s ' % description)
+
+    e = atoms.get_potential_energy()
+    EE = e
+
+    c = {}
+    for atom in atoms:
+        if atom.symbol in c:
+            c[atom.symbol] += 1
+        else:
+            c[atom.symbol] = 1
+
+    x = {}
+    for key in c:
+        x[key] = c[key]/float(len(atoms)) # mole fractions
+
+    # compute formation energy
+    for atom in atoms:
+        e -= reference_energies[atom.symbol][1]
+
+    # construct the string for reax to compute formation energy
+    s = '%1.2f + %s /1' % (weight, description)
+    for entry in reference_energies:
+        desc, ref_energy = reference_energies[entry]
+        xi = x.get(entry,None)
+        if xi:
+            s += ' - %s / %f ' % (desc, 1./(float(xi)*len(atoms)))
+            EE -= ref_energy / (1./(float(xi)*len(atoms)))
+    s += '  %f' % (e*23.061)
+
+    cells['ENERGY'].append(s)
+    return cells
+
+def eos(atoms1, atoms2, weight=1.0, cells=None):
+    '''
+    add equation of state data to the trainset.in
+    '''
+    d1 = reax_hash(atoms1)
+    d2 = reax_hash(atoms2)
+    e = atoms1.get_potential_energy() - atoms2.get_potential_energy()
+    s = '%1.2f + %s /1 - %s /1 %f' % (weight,
+                                      d1, d2,
+                                      e*23.061)
+    cells['ENERGY'].append(s)
+    return cells