Finalize last written frame set when destroying trajectory. Fix bug reading data interval.

2 files changed, 188 insertions, 64 deletions

diff --git a/src/lib/tng_io.c b/src/lib/tng_io.c
index c1960ea..9367ef5 100644
--- a/src/lib/tng_io.c
+++ b/src/lib/tng_io.c
@@ -138,6 +138,10 @@ struct tng_trajectory_frame_set {
     int64_t first_frame;
     /** The number of frames in this frame set */
     int64_t n_frames;
+    /** The number of written frames in this frame set (used when writing one
+     * frame at a time. */
+    int64_t n_written_frames;
+
     /** A list of the number of each molecule type - only used when using
      * variable number of atoms */
     int64_t *molecule_cnt_list;
@@ -550,7 +554,7 @@ static tng_function_status hash_match_verify(tng_gen_block_t block,
     md5_state_t md5_state;
     char hash[TNG_HASH_LEN];
 
-    if(strlen(block->hash) == 0)
+    if(strncmp(block->hash, "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0", 16) == 0)
     {
         *results = TNG_TRUE;
         return(TNG_FAILURE);
@@ -5670,6 +5674,101 @@ static tng_function_status tng_frame_set_pointers_update
     return(TNG_SUCCESS);
 }
 
+/** Finish writing the current frame set. Update the number of frames
+ * and the hashes of the frame set and all its data blocks (if hash_mode
+ * == TNG_USE_HASH).
+ * @param tng_data is a trajectory data container.
+ * @param hash_mode specifies whether to update the block md5 hash when
+ * updating the pointers.
+ * @return TNG_SUCCESS (0) if successful or TNG_CRITICAL (2) if a major
+ * error has occured.
+ */
+static tng_function_status tng_frame_set_finalize
+                (tng_trajectory_t tng_data, const tng_hash_mode hash_mode)
+{
+    tng_gen_block_t block;
+    tng_trajectory_frame_set_t frame_set;
+    FILE *temp = tng_data->input_file;
+    int64_t pos, output_file_pos, contents_start_pos,
+            output_file_len;
+
+    frame_set = &tng_data->current_trajectory_frame_set;
+
+    if(frame_set->n_written_frames == frame_set->n_frames)
+    {
+        return(TNG_SUCCESS);
+    }
+
+    if(tng_output_file_init(tng_data) != TNG_SUCCESS)
+    {
+        printf("Cannot initialise destination file. %s: %d\n",
+               __FILE__, __LINE__);
+        return(TNG_CRITICAL);
+    }
+
+    tng_block_init(&block);
+    output_file_pos = ftell(tng_data->output_file);
+
+    tng_data->input_file = tng_data->output_file;
+
+    pos = tng_data->current_trajectory_frame_set_output_file_pos;
+
+    fseek(tng_data->output_file, pos, SEEK_SET);
+
+    if(tng_block_header_read(tng_data, block) != TNG_SUCCESS)
+    {
+        printf("Cannot read frame set header. %s: %d\n",
+            __FILE__, __LINE__);
+        tng_data->input_file = temp;
+        tng_block_destroy(&block);
+        return(TNG_CRITICAL);
+    }
+
+    contents_start_pos = ftell(tng_data->output_file);
+
+    fseek(tng_data->output_file, sizeof(frame_set->first_frame), SEEK_CUR);
+    if(fwrite(&frame_set->n_written_frames, sizeof(frame_set->n_frames),
+              1, tng_data->output_file) != 1)
+    {
+        tng_data->input_file = temp;
+        tng_block_destroy(&block);
+        return(TNG_CRITICAL);
+    }
+
+
+    if(hash_mode == TNG_SKIP_HASH)
+    {
+        fseek(tng_data->output_file, output_file_pos, SEEK_SET);
+        return(TNG_SUCCESS);
+    }
+
+    tng_md5_hash_update(tng_data, block, pos, pos + block->header_contents_size);
+
+    fseek(tng_data->output_file, 0, SEEK_END);
+    output_file_len = ftell(tng_data->output_file);
+    pos = contents_start_pos + block->block_contents_size;
+    fseek(tng_data->output_file, pos,
+          SEEK_SET);
+
+    while(pos < output_file_len)
+    {
+        if(tng_block_header_read(tng_data, block) != TNG_SUCCESS)
+        {
+            printf("Cannot read block header. %s: %d\n", __FILE__, __LINE__);
+            tng_data->input_file = temp;
+            tng_block_destroy(&block);
+            return(TNG_CRITICAL);
+        }
+        tng_md5_hash_update(tng_data, block, pos,
+                            pos + block->header_contents_size);
+        pos += block->header_contents_size + block->block_contents_size;
+        fseek(tng_data->output_file, pos, SEEK_SET);
+    }
+
+    return(TNG_SUCCESS);
+}
+
+
 /** Sets the name of a file contents block
  * @param tng_data is a trajectory data container.
  * @param block is the block, of which to change names.
@@ -6868,6 +6967,8 @@ tng_function_status tng_trajectory_destroy(tng_trajectory_t *tng_data_p)
 
     if(tng_data->output_file)
     {
+        /* FIXME: Do not always write the hash */
+        tng_frame_set_finalize(tng_data, TNG_USE_HASH);
         fclose(tng_data->output_file);
         tng_data->output_file = 0;
     }
@@ -9224,6 +9325,7 @@ tng_function_status tng_frame_set_new(tng_trajectory_t tng_data,
 
     frame_set->first_frame = first_frame;
     frame_set->n_frames = n_frames;
+    frame_set->n_written_frames = 0;
 
     if(tng_data->first_trajectory_frame_set_output_file_pos == -1 ||
        tng_data->first_trajectory_frame_set_output_file_pos == 0)
@@ -9366,6 +9468,11 @@ tng_function_status tng_data_block_add(tng_trajectory_t tng_data,
     {
         orig = new_data;
 
+        if(n_frames > frame_set->n_written_frames)
+        {
+            frame_set->n_written_frames = n_frames;
+        }
+
         switch(datatype)
         {
         case TNG_FLOAT_DATA:
@@ -9567,6 +9674,11 @@ tng_function_status tng_particle_data_block_add(tng_trajectory_t tng_data,
     {
         orig = new_data;
 
+        if(n_frames > frame_set->n_written_frames)
+        {
+            frame_set->n_written_frames = n_frames;
+        }
+
         switch(datatype)
         {
         case TNG_FLOAT_DATA:
@@ -10035,6 +10147,12 @@ tng_function_status tng_frame_data_write(tng_trajectory_t tng_data,
 
     fflush(tng_data->output_file);
 
+    /* Update the number of written frames in the frame set. */
+    if(frame_nr - frame_set->first_frame + 1 > frame_set->n_written_frames)
+    {
+        frame_set->n_written_frames = frame_nr - frame_set->first_frame + 1;
+    }
+
     /* If the last frame has been written update the hash */
     if(hash_mode == TNG_USE_HASH && (frame_nr + data.stride_length -
        data.first_frame_with_data) >=
@@ -10554,6 +10672,12 @@ tng_function_status tng_frame_particle_data_write(tng_trajectory_t tng_data,
     }
     fflush(tng_data->output_file);
 
+    /* Update the number of written frames in the frame set. */
+    if(frame_nr - frame_set->first_frame + 1 > frame_set->n_written_frames)
+    {
+        frame_set->n_written_frames = frame_nr - frame_set->first_frame + 1;
+    }
+
     /* If the last frame has been written update the hash */
     if(hash_mode == TNG_USE_HASH && (frame_nr + data.stride_length -
        data.first_frame_with_data) >=
@@ -10906,13 +11030,6 @@ tng_function_status tng_data_interval_get(tng_trajectory_t tng_data,
     case TNG_CHAR_DATA:
         for(i=0; i<n_frames; i++)
         {
-            for(j=*n_values_per_frame; j--;)
-            {
-                len = strlen(data->values[current_frame_pos][j].c) + 1;
-                (*values)[i][j].c = malloc(len);
-                strncpy((*values)[i][j].c, data->values[current_frame_pos][j].c, len);
-            }
-            current_frame_pos++;
             if(current_frame_pos == frame_set->n_frames)
             {
                 stat = tng_frame_set_read_next(tng_data, hash_mode);
@@ -10922,16 +11039,18 @@ tng_function_status tng_data_interval_get(tng_trajectory_t tng_data,
                 }
                 current_frame_pos = 0;
             }
+            for(j=*n_values_per_frame; j--;)
+            {
+                len = strlen(data->values[current_frame_pos][j].c) + 1;
+                (*values)[i][j].c = malloc(len);
+                strncpy((*values)[i][j].c, data->values[current_frame_pos][j].c, len);
+            }
+            current_frame_pos++;
         }
         break;
     case TNG_INT_DATA:
         for(i=0; i<n_frames; i++)
         {
-            for(j=*n_values_per_frame; j--;)
-            {
-                (*values)[i][j].i = data->values[current_frame_pos][j].i;
-            }
-            current_frame_pos++;
             if(current_frame_pos == frame_set->n_frames)
             {
                 stat = tng_frame_set_read_next(tng_data, hash_mode);
@@ -10941,16 +11060,16 @@ tng_function_status tng_data_interval_get(tng_trajectory_t tng_data,
                 }
                 current_frame_pos = 0;
             }
+            for(j=*n_values_per_frame; j--;)
+            {
+                (*values)[i][j].i = data->values[current_frame_pos][j].i;
+            }
+            current_frame_pos++;
         }
         break;
     case TNG_FLOAT_DATA:
         for(i=0; i<n_frames; i++)
         {
-            for(j=*n_values_per_frame; j--;)
-            {
-                (*values)[i][j].f = data->values[current_frame_pos][j].f;
-            }
-            current_frame_pos++;
             if(current_frame_pos == frame_set->n_frames)
             {
                 stat = tng_frame_set_read_next(tng_data, hash_mode);
@@ -10960,17 +11079,17 @@ tng_function_status tng_data_interval_get(tng_trajectory_t tng_data,
                 }
                 current_frame_pos = 0;
             }
+            for(j=*n_values_per_frame; j--;)
+            {
+                (*values)[i][j].f = data->values[current_frame_pos][j].f;
+            }
+            current_frame_pos++;
         }
         break;
     case TNG_DOUBLE_DATA:
     default:
         for(i=0; i<n_frames; i++)
         {
-            for(j=*n_values_per_frame; j--;)
-            {
-                (*values)[i][j].d = data->values[current_frame_pos][j].d;
-            }
-            current_frame_pos++;
             if(current_frame_pos == frame_set->n_frames)
             {
                 stat = tng_frame_set_read_next(tng_data, hash_mode);
@@ -10980,6 +11099,11 @@ tng_function_status tng_data_interval_get(tng_trajectory_t tng_data,
                 }
                 current_frame_pos = 0;
             }
+            for(j=*n_values_per_frame; j--;)
+            {
+                (*values)[i][j].d = data->values[current_frame_pos][j].d;
+            }
+            current_frame_pos++;
         }
     }
 
@@ -11303,6 +11427,15 @@ tng_function_status tng_particle_data_interval_get(tng_trajectory_t tng_data,
     case TNG_CHAR_DATA:
         for(i=0; i<n_frames; i++)
         {
+            if(current_frame_pos == frame_set->n_frames)
+            {
+                stat = tng_frame_set_read_next(tng_data, hash_mode);
+                if(stat != TNG_SUCCESS)
+                {
+                    return(stat);
+                }
+                current_frame_pos = 0;
+            }
             for(j=*n_particles; j--;)
             {
                 tng_particle_mapping_get_real_particle(frame_set, j, &mapping);
@@ -11314,6 +11447,11 @@ tng_function_status tng_particle_data_interval_get(tng_trajectory_t tng_data,
                 }
             }
             current_frame_pos++;
+        }
+        break;
+    case TNG_INT_DATA:
+        for(i=0; i<n_frames; i++)
+        {
             if(current_frame_pos == frame_set->n_frames)
             {
                 stat = tng_frame_set_read_next(tng_data, hash_mode);
@@ -11323,11 +11461,6 @@ tng_function_status tng_particle_data_interval_get(tng_trajectory_t tng_data,
                 }
                 current_frame_pos = 0;
             }
-        }
-        break;
-    case TNG_INT_DATA:
-        for(i=0; i<n_frames; i++)
-        {
             for(j=*n_particles; j--;)
             {
                 tng_particle_mapping_get_real_particle(frame_set, j, &mapping);
@@ -11337,6 +11470,11 @@ tng_function_status tng_particle_data_interval_get(tng_trajectory_t tng_data,
                 }
             }
             current_frame_pos++;
+        }
+        break;
+    case TNG_FLOAT_DATA:
+        for(i=0; i<n_frames; i++)
+        {
             if(current_frame_pos == frame_set->n_frames)
             {
                 stat = tng_frame_set_read_next(tng_data, hash_mode);
@@ -11346,11 +11484,6 @@ tng_function_status tng_particle_data_interval_get(tng_trajectory_t tng_data,
                 }
                 current_frame_pos = 0;
             }
-        }
-        break;
-    case TNG_FLOAT_DATA:
-        for(i=0; i<n_frames; i++)
-        {
             for(j=*n_particles; j--;)
             {
                 tng_particle_mapping_get_real_particle(frame_set, j, &mapping);
@@ -11360,6 +11493,12 @@ tng_function_status tng_particle_data_interval_get(tng_trajectory_t tng_data,
                 }
             }
             current_frame_pos++;
+        }
+        break;
+    case TNG_DOUBLE_DATA:
+    default:
+        for(i=0; i<n_frames; i++)
+        {
             if(current_frame_pos == frame_set->n_frames)
             {
                 stat = tng_frame_set_read_next(tng_data, hash_mode);
@@ -11369,12 +11508,6 @@ tng_function_status tng_particle_data_interval_get(tng_trajectory_t tng_data,
                 }
                 current_frame_pos = 0;
             }
-        }
-        break;
-    case TNG_DOUBLE_DATA:
-    default:
-        for(i=0; i<n_frames; i++)
-        {
             for(j=*n_particles; j--;)
             {
                 tng_particle_mapping_get_real_particle(frame_set, j, &mapping);
@@ -11384,15 +11517,6 @@ tng_function_status tng_particle_data_interval_get(tng_trajectory_t tng_data,
                 }
             }
             current_frame_pos++;
-            if(current_frame_pos == frame_set->n_frames)
-            {
-                stat = tng_frame_set_read_next(tng_data, hash_mode);
-                if(stat != TNG_SUCCESS)
-                {
-                    return(stat);
-                }
-                current_frame_pos = 0;
-            }
         }
     }
 
diff --git a/src/tests/md_openmp.c b/src/tests/md_openmp.c
index a779159..5572746 100644
--- a/src/tests/md_openmp.c
+++ b/src/tests/md_openmp.c
@@ -7,14 +7,14 @@
 #include "tng_io_testing.h"
 
 int main ( int argc, char *argv[] );
-void compute ( int np, int nd, double pos[], double vel[], 
+void compute ( int np, int nd, double pos[], double vel[],
     double mass, double f[], double *pot, double *kin );
 double dist ( int nd, double r1[], double r2[], double dr[] );
-void initialize ( int np, int nd, double box[], int *seed, double pos[], 
+void initialize ( int np, int nd, double box[], int *seed, double pos[],
     double vel[], double acc[] );
 double r8_uniform_01 ( int *seed );
 void timestamp ( void );
-void update ( int np, int nd, double pos[], double vel[], double f[], 
+void update ( int np, int nd, double pos[], double vel[], double f[],
     double acc[], double mass, double dt );
 
 /******************************************************************************/
@@ -75,7 +75,7 @@ int main ( int argc, char *argv[] )
     int proc_num;
     int seed = 123456789;
     int step;
-    int step_num = 1000;
+    int step_num = 1050;
     int step_print;
     int step_print_index;
     int step_print_num;
@@ -117,7 +117,7 @@ int main ( int argc, char *argv[] )
     printf ( "  Number of processors available = %d\n", proc_num );
     printf ( "  Number of threads =              %d\n", omp_get_max_threads ( ) );
 
-    
+
     printf("\n");
     printf("  Initializing trajectory storage.\n");
     if(tng_trajectory_init(&traj) != TNG_SUCCESS)
@@ -147,7 +147,7 @@ int main ( int argc, char *argv[] )
     }
     tng_molecule_cnt_set(traj, molecule, np);
 
-    
+
 /*
     Set the dimensions of the box.
 */
@@ -190,7 +190,7 @@ int main ( int argc, char *argv[] )
     compute ( np, nd, pos, vel, mass, force, &potential, &kinetic );
 
     e0 = potential + kinetic;
-    
+
     /* Saving frequency */
     step_save = 5;
 
@@ -232,7 +232,7 @@ int main ( int argc, char *argv[] )
                 i, __FILE__, __LINE__);
         exit(1);
     }
-    
+
     /* Add empty data blocks */
     if(tng_particle_data_block_add(traj, TNG_TRAJ_POSITIONS,
                                 "POSITIONS",
@@ -269,7 +269,7 @@ int main ( int argc, char *argv[] )
         printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
         exit(1);
     }
-    
+
     /* There is no standard ID for potential energy. Pick one. The
        potential energy will not be saved every frame - it is sparsely
        saved. */
@@ -284,14 +284,14 @@ int main ( int argc, char *argv[] )
         printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
         exit(1);
     }
-    
+
     /* Write the frame set to disk */
     if(tng_frame_set_write(traj, TNG_USE_HASH) != TNG_SUCCESS)
     {
         printf("Error writing frame set. %s: %d\n", __FILE__, __LINE__);
         exit(1);
     }
-        
+
     wtime = omp_get_wtime ( );
 
     for ( step = 1; step <= step_num; step++ )
@@ -368,7 +368,7 @@ int main ( int argc, char *argv[] )
 }
 /******************************************************************************/
 
-void compute ( int np, int nd, double pos[], double vel[], 
+void compute ( int np, int nd, double pos[], double vel[],
     double mass, double f[], double *pot, double *kin )
 
 /******************************************************************************/
@@ -439,7 +439,7 @@ void compute ( int np, int nd, double pos[], double vel[],
 # pragma omp parallel \
     shared ( f, nd, np, pos, vel ) \
     private ( i, j, k, rij, d, d2 )
-    
+
 
 # pragma omp for reduction ( + : pe, ke )
     for ( k = 0; k < np; k++ )
@@ -487,7 +487,7 @@ void compute ( int np, int nd, double pos[], double vel[],
     }
 
     ke = ke * 0.5 * mass;
-    
+
     *pot = pe;
     *kin = ke;
 
@@ -542,7 +542,7 @@ double dist ( int nd, double r1[], double r2[], double dr[] )
 }
 /******************************************************************************/
 
-void initialize ( int np, int nd, double box[], int *seed, double pos[], 
+void initialize ( int np, int nd, double box[], int *seed, double pos[],
     double vel[], double acc[] )
 
 /******************************************************************************/
@@ -729,7 +729,7 @@ void timestamp ( void )
 }
 /******************************************************************************/
 
-void update ( int np, int nd, double pos[], double vel[], double f[], 
+void update ( int np, int nd, double pos[], double vel[], double f[],
     double acc[], double mass, double dt )
 
 /******************************************************************************/