More work on high-level API. New example files.

Fixed many bugs in the high-level API.
Use the high-level API (where appropriate) in new example files.

6 files changed, 964 insertions, 65 deletions

diff --git a/include/tng_io.h b/include/tng_io.h
index 57ae8a4..1085d17 100644
--- a/include/tng_io.h
+++ b/include/tng_io.h
@@ -1729,6 +1729,7 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get
                  const int64_t end_frame_nr,
                  const tng_hash_mode hash_mode,
                  union data_values ***values,
+                 int64_t *stride_length,
                  int64_t *n_values_per_frame,
                  tng_data_type *type);
 
@@ -1875,6 +1876,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get
                  const tng_hash_mode hash_mode,
                  void **values,
                  int64_t *n_particles,
+                 int64_t *stride_length,
                  int64_t *n_values_per_frame,
                  tng_data_type *type);
 
@@ -1933,43 +1935,51 @@ tng_function_status DECLSPECDLLEXPORT tng_util_molecule_particles_set
 
 tng_function_status DECLSPECDLLEXPORT tng_util_pos_read
                 (tng_trajectory_t tng_data,
-                 float *positions);
+                 float **positions,
+                 int64_t *stride_length);
 
 tng_function_status DECLSPECDLLEXPORT tng_util_vel_read
                 (tng_trajectory_t tng_data,
-                 float *velocities);
+                 float **velocities,
+                 int64_t *stride_length);
 
 tng_function_status DECLSPECDLLEXPORT tng_util_force_read
                 (tng_trajectory_t tng_data,
-                 float *forces);
+                 float **forces,
+                 int64_t *stride_length);
 
 tng_function_status DECLSPECDLLEXPORT tng_util_box_shape_read
                 (tng_trajectory_t tng_data,
-                 float *box_shape);
+                 float **box_shape,
+                 int64_t *stride_length);
 
 tng_function_status DECLSPECDLLEXPORT tng_util_pos_read_range
                 (tng_trajectory_t tng_data,
                  const int64_t first_frame,
                  const int64_t last_frame,
-                 float *positions);
+                 float **positions,
+                 int64_t *stride_length);
 
 tng_function_status DECLSPECDLLEXPORT tng_util_vel_read_range
                 (tng_trajectory_t tng_data,
                  const int64_t first_frame,
                  const int64_t last_frame,
-                 float *velocities);
+                 float **velocities,
+                 int64_t *stride_length);
 
 tng_function_status DECLSPECDLLEXPORT tng_util_force_read_range
                 (tng_trajectory_t tng_data,
                  const int64_t first_frame,
                  const int64_t last_frame,
-                 float *forces);
+                 float **forces,
+                 int64_t *stride_length);
 
 tng_function_status DECLSPECDLLEXPORT tng_util_box_shape_read_range
                 (tng_trajectory_t tng_data,
                  const int64_t first_frame,
                  const int64_t last_frame,
-                 float *box_shape);
+                 float **box_shape,
+                 int64_t *stride_length);
 
 tng_function_status DECLSPECDLLEXPORT tng_util_pos_write_frequence_set
                 (tng_trajectory_t tng_data,
diff --git a/src/lib/tng_io.c b/src/lib/tng_io.c
index 4a76d25..3397fcd 100644
--- a/src/lib/tng_io.c
+++ b/src/lib/tng_io.c
@@ -147,6 +147,10 @@ struct tng_trajectory_frame_set {
     /** The number of written frames in this frame set (used when writing one
      * frame at a time). */
     int64_t n_written_frames;
+    /** The number of frames not yet written to file in this frame set
+     * (used from the utility functions to finish the writing properly. */
+    int64_t n_unwritten_frames;
+
 
     /** A list of the number of each molecule type - only used when using
      * variable number of atoms */
@@ -3814,10 +3818,25 @@ static tng_function_status tng_uncompress(tng_trajectory_t tng_data,
 
     if(type == TNG_FLOAT_DATA)
     {
+        f_dest = malloc(uncompressed_len);
+        if(!f_dest)
+        {
+            printf("Cannot allocate memory (%"PRId64" bytes). %s: %d\n",
+                uncompressed_len, __FILE__, __LINE__);
+            return(TNG_CRITICAL);
+        }
+        
         result = tng_compress_uncompress_float(start_pos, f_dest);
     }
     else
     {
+        d_dest = malloc(uncompressed_len);
+        if(!d_dest)
+        {
+            printf("Cannot allocate memory (%"PRId64" bytes). %s: %d\n",
+                uncompressed_len, __FILE__, __LINE__);
+            return(TNG_CRITICAL);
+        }
         result = tng_compress_uncompress(start_pos, d_dest);
     }
 
@@ -4322,7 +4341,8 @@ static tng_function_status tng_particle_data_read
         data->compression_multiplier = multiplier;
     }
 
-    if(block_type_flag == TNG_TRAJECTORY_BLOCK &&
+    if(/*block_type_flag == TNG_TRAJECTORY_BLOCK &&*/
+       tng_data->current_trajectory_frame_set_input_file_pos > 0 &&
        tng_data->var_num_atoms_flag)
     {
         tot_n_particles = frame_set->n_particles;
@@ -4349,8 +4369,6 @@ static tng_function_status tng_particle_data_read
                        __FILE__, __LINE__);
                 return(TNG_CRITICAL);
             }
-
-            printf("TNG compression not implemented yet.\n");
             break;
     #ifdef USE_ZLIB
         case TNG_GZIP_COMPRESSION:
@@ -4891,6 +4909,8 @@ static tng_function_status tng_particle_data_block_write
         memset(block->block_contents+offset, 0, block->block_contents_size - offset);
     }
 
+    frame_set->n_written_frames += frame_set->n_unwritten_frames;
+    frame_set->n_unwritten_frames = 0;
 
     if(frame_set->n_written_frames > 0)
     {
@@ -7491,6 +7511,7 @@ tng_function_status DECLSPECDLLEXPORT tng_trajectory_init(tng_trajectory_t *tng_
     frame_set->tr_data = 0;
 
     frame_set->n_written_frames = 0;
+    frame_set->n_unwritten_frames = 0;
 
     frame_set->next_frame_set_file_pos = -1;
     frame_set->prev_frame_set_file_pos = -1;
@@ -10045,6 +10066,8 @@ tng_function_status tng_frame_set_write(tng_trajectory_t tng_data,
 
     tng_block_destroy(&block);
 
+    frame_set->n_unwritten_frames = 0;
+
     return(stat);
 }
 
@@ -10218,6 +10241,7 @@ tng_function_status DECLSPECDLLEXPORT tng_frame_set_new
     frame_set->first_frame = first_frame;
     frame_set->n_frames = n_frames;
     frame_set->n_written_frames = 0;
+    frame_set->n_unwritten_frames = 0;
 
     if(tng_data->first_trajectory_frame_set_output_file_pos == -1 ||
        tng_data->first_trajectory_frame_set_output_file_pos == 0)
@@ -11892,7 +11916,7 @@ tng_function_status tng_data_vector_get(tng_trajectory_t tng_data,
 
     *values = temp;
 
-    memcpy(values, data->values, data_size);
+    memcpy(*values, data->values, data_size);
 
     return(TNG_SUCCESS);
 }
@@ -12090,11 +12114,12 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get
                  const int64_t end_frame_nr,
                  const tng_hash_mode hash_mode,
                  union data_values ***values,
+                 int64_t *stride_length,
                  int64_t *n_values_per_frame,
                  tng_data_type *type)
 {
     int64_t n_frames, tot_n_frames;
-    int64_t stride_length, current_frame_pos, data_size, frame_size;
+    int64_t current_frame_pos, data_size, frame_size;
     int size;
     tng_trajectory_frame_set_t frame_set;
     tng_function_status stat;
@@ -12110,7 +12135,7 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get
     frame_set = &tng_data->current_trajectory_frame_set;
 
     stat = tng_data_vector_get(tng_data, block_id, &current_values,
-                               &n_frames, &stride_length,
+                               &n_frames, stride_length,
                                n_values_per_frame, type);
 
     if(stat != TNG_SUCCESS)
@@ -12138,7 +12163,7 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get
         size = sizeof(double);
     }
 
-    data_size = (tot_n_frames / stride_length) * size *
+    data_size = (tot_n_frames / *stride_length) * size *
                 (*n_values_per_frame);
 
     temp = realloc(*values, data_size);
@@ -12158,7 +12183,7 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get
 
     memcpy(*values, current_values + current_frame_pos * frame_size,
            frame_size * (frame_set->n_frames - frame_set->first_frame -
-           current_frame_pos));
+           current_frame_pos) / *stride_length);
 
     current_frame_pos += frame_set->n_frames - frame_set->first_frame -
                          current_frame_pos;
@@ -12172,7 +12197,7 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get
         }
 
         stat = tng_data_vector_get(tng_data, block_id, &current_values,
-                                   &n_frames, &stride_length,
+                                   &n_frames, stride_length,
                                    n_values_per_frame, type);
 
         if(stat != TNG_SUCCESS)
@@ -12186,10 +12211,10 @@ tng_function_status DECLSPECDLLEXPORT tng_data_vector_interval_get
             return(stat);
         }
 
-        memcpy(*values + (current_frame_pos - start_frame_nr) / stride_length,
+        memcpy(*values + (current_frame_pos - start_frame_nr) / *stride_length,
                current_values,
                frame_size * (frame_set->n_frames - frame_set->first_frame -
-               current_frame_pos) / stride_length);
+               current_frame_pos) / *stride_length);
 
         current_frame_pos += frame_set->n_frames;
     }
@@ -12401,15 +12426,6 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_get
 
     if(tng_particle_data_find(tng_data, block_id, &data) != TNG_SUCCESS)
     {
-        if(tng_data->current_trajectory_frame_set_input_file_pos > 0)
-        {
-            block_type_flag = TNG_TRAJECTORY_BLOCK;
-        }
-        else
-        {
-            block_type_flag = TNG_NON_TRAJECTORY_BLOCK;
-        }
-
         tng_block_init(&block);
         file_pos = ftell(tng_data->input_file);
         /* Read all blocks until next frame set block */
@@ -12444,7 +12460,6 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_get
             if(data->block_id == block_id)
             {
                 block_index = i;
-                block_type_flag = TNG_TRAJECTORY_BLOCK;
                 break;
             }
         }
@@ -12454,8 +12469,17 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_get
         }
     }
 
-    if(block_type_flag == TNG_TRAJECTORY_BLOCK &&
-       tng_data->var_num_atoms_flag)
+    if(tng_data->current_trajectory_frame_set_input_file_pos > 0)
+    {
+        block_type_flag = TNG_TRAJECTORY_BLOCK;
+    }
+    else
+    {
+        block_type_flag = TNG_NON_TRAJECTORY_BLOCK;
+    }
+
+   if(block_type_flag == TNG_TRAJECTORY_BLOCK &&
+      tng_data->var_num_atoms_flag)
     {
         *n_particles = frame_set->n_particles;
     }
@@ -12501,7 +12525,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_get
 
     if(frame_set->n_mapping_blocks <= 0)
     {
-        memcpy(values, data->values, data_size);
+        memcpy(*values, data->values, data_size);
     }
     else
     {
@@ -12511,7 +12535,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_get
             for(j = *n_particles; j--;)
             {
                 tng_particle_mapping_get_real_particle(frame_set, j, &mapping);
-                memcpy(values + size * (i * i_step + mapping *
+                memcpy(*values + size * (i * i_step + mapping *
                        (*n_values_per_frame)),
                        data->values + size *
                        (i * i_step + j * (*n_values_per_frame)),
@@ -12756,11 +12780,12 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get
                  const tng_hash_mode hash_mode,
                  void **values,
                  int64_t *n_particles,
+                 int64_t *stride_length,
                  int64_t *n_values_per_frame,
                  tng_data_type *type)
 {
     int64_t n_frames, tot_n_frames;
-    int64_t stride_length, current_frame_pos, data_size, frame_size;
+    int64_t current_frame_pos, last_frame_pos, data_size, frame_size;
     int size;
     tng_trajectory_frame_set_t frame_set;
     tng_function_status stat;
@@ -12775,7 +12800,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get
     frame_set = &tng_data->current_trajectory_frame_set;
 
     stat = tng_particle_data_vector_get(tng_data, block_id, &current_values,
-                                        &n_frames, &stride_length, n_particles,
+                                        &n_frames, stride_length, n_particles,
                                         n_values_per_frame, type);
 
     if(stat != TNG_SUCCESS)
@@ -12787,7 +12812,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get
         return(stat);
     }
 
-    tot_n_frames = end_frame_nr - start_frame_nr + 1;
+    tot_n_frames = end_frame_nr - start_frame_nr + 1 *(*stride_length);
     switch(*type)
     {
     case TNG_CHAR_DATA:
@@ -12803,7 +12828,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get
         size = sizeof(double);
     }
 
-    data_size = (tot_n_frames / stride_length) * size * (*n_particles) *
+    data_size = (tot_n_frames / *stride_length) * size * (*n_particles) *
                 (*n_values_per_frame);
 
     temp = realloc(*values, data_size);
@@ -12821,11 +12846,15 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get
 
     frame_size = size * (*n_particles) * (*n_values_per_frame);
 
-    memcpy(*values, current_values + current_frame_pos * frame_size,
-           frame_size * (frame_set->n_frames - frame_set->first_frame -
-           current_frame_pos));
+    last_frame_pos = tng_min(frame_set->first_frame + n_frames -
+                             *stride_length,
+                             end_frame_nr - frame_set->first_frame);
 
-    current_frame_pos += frame_set->n_frames - frame_set->first_frame -
+    memcpy(*values, current_values + current_frame_pos * frame_size /
+           *stride_length,
+           frame_size * (last_frame_pos / *stride_length + 1));
+
+    current_frame_pos += n_frames - frame_set->first_frame -
                          current_frame_pos;
 
     while(current_frame_pos <= end_frame_nr)
@@ -12837,7 +12866,7 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get
         }
 
         stat = tng_particle_data_vector_get(tng_data, block_id, &current_values,
-                                            &n_frames, &stride_length, n_particles,
+                                            &n_frames, stride_length, n_particles,
                                             n_values_per_frame, type);
 
         if(stat != TNG_SUCCESS)
@@ -12851,12 +12880,16 @@ tng_function_status DECLSPECDLLEXPORT tng_particle_data_vector_interval_get
             return(stat);
         }
 
-        memcpy(*values + (current_frame_pos - start_frame_nr) / stride_length,
+        last_frame_pos = tng_min(frame_set->first_frame + n_frames -
+                                 *stride_length,
+                                 end_frame_nr - frame_set->first_frame);
+        
+        memcpy(*values + (current_frame_pos - start_frame_nr) *frame_size /
+               *stride_length,
                current_values,
-               frame_size * (frame_set->n_frames - frame_set->first_frame -
-               current_frame_pos) / stride_length);
+               frame_size * (last_frame_pos / *stride_length + 1));
 
-        current_frame_pos += frame_set->n_frames;
+        current_frame_pos += n_frames;
     }
 
     free(current_values);
@@ -12919,6 +12952,15 @@ tng_function_status DECLSPECDLLEXPORT tng_util_trajectory_open
 tng_function_status DECLSPECDLLEXPORT tng_util_trajectory_close
                 (tng_trajectory_t *tng_data_p)
 {
+    tng_trajectory_frame_set_t frame_set = &(*tng_data_p)->
+                                           current_trajectory_frame_set;
+
+    if(frame_set->n_unwritten_frames > 0)
+    {
+        frame_set->n_frames = frame_set->n_unwritten_frames;
+        tng_frame_set_write(*tng_data_p, TNG_USE_HASH);
+    }
+    
     return(tng_trajectory_destroy(tng_data_p));
 }
 
@@ -13049,7 +13091,7 @@ tng_function_status DECLSPECDLLEXPORT tng_util_molecule_particles_set
 
 tng_function_status DECLSPECDLLEXPORT tng_util_pos_read
                 (tng_trajectory_t tng_data,
-                 float *positions)
+                 float **positions, int64_t *stride_length)
 {
     int64_t n_frames, n_particles, n_values_per_frame;
     tng_data_type type;
@@ -13063,8 +13105,9 @@ tng_function_status DECLSPECDLLEXPORT tng_util_pos_read
 
     stat = tng_particle_data_vector_interval_get(tng_data, TNG_TRAJ_POSITIONS,
                                                  0, n_frames - 1, TNG_USE_HASH,
-                                                 (void **)(&positions),
+                                                 (void **)positions,
                                                  &n_particles,
+                                                 stride_length,
                                                  &n_values_per_frame,
                                                  &type);
 
@@ -13073,7 +13116,7 @@ tng_function_status DECLSPECDLLEXPORT tng_util_pos_read
 
 tng_function_status DECLSPECDLLEXPORT tng_util_vel_read
                 (tng_trajectory_t tng_data,
-                 float *velocities)
+                 float **velocities, int64_t *stride_length)
 {
     int64_t n_frames, n_particles, n_values_per_frame;
     tng_data_type type;
@@ -13087,8 +13130,9 @@ tng_function_status DECLSPECDLLEXPORT tng_util_vel_read
 
     stat = tng_particle_data_vector_interval_get(tng_data, TNG_TRAJ_VELOCITIES,
                                                  0, n_frames - 1, TNG_USE_HASH,
-                                                 (void **)(&velocities),
+                                                 (void **)velocities,
                                                  &n_particles,
+                                                 stride_length,
                                                  &n_values_per_frame,
                                                  &type);
 
@@ -13097,7 +13141,7 @@ tng_function_status DECLSPECDLLEXPORT tng_util_vel_read
 
 tng_function_status DECLSPECDLLEXPORT tng_util_force_read
                 (tng_trajectory_t tng_data,
-                 float *forces)
+                 float **forces, int64_t *stride_length)
 {
     int64_t n_frames, n_particles, n_values_per_frame;
     tng_data_type type;
@@ -13111,8 +13155,9 @@ tng_function_status DECLSPECDLLEXPORT tng_util_force_read
 
     stat = tng_particle_data_vector_interval_get(tng_data, TNG_TRAJ_FORCES,
                                                  0, n_frames - 1, TNG_USE_HASH,
-                                                 (void **)(&forces),
+                                                 (void **)forces,
                                                  &n_particles,
+                                                 stride_length,
                                                  &n_values_per_frame,
                                                  &type);
 
@@ -13123,7 +13168,8 @@ tng_function_status DECLSPECDLLEXPORT tng_util_pos_read_range
                 (tng_trajectory_t tng_data,
                  const int64_t first_frame,
                  const int64_t last_frame,
-                 float *positions)
+                 float **positions,
+                 int64_t *stride_length)
 {
     int64_t n_particles, n_values_per_frame;
     tng_data_type type;
@@ -13132,8 +13178,9 @@ tng_function_status DECLSPECDLLEXPORT tng_util_pos_read_range
     stat = tng_particle_data_vector_interval_get(tng_data, TNG_TRAJ_POSITIONS,
                                                  first_frame, last_frame,
                                                  TNG_USE_HASH,
-                                                 (void **)(&positions),
+                                                 (void **)positions,
                                                  &n_particles,
+                                                 stride_length,
                                                  &n_values_per_frame,
                                                  &type);
 
@@ -13144,7 +13191,8 @@ tng_function_status DECLSPECDLLEXPORT tng_util_vel_read_range
                 (tng_trajectory_t tng_data,
                  const int64_t first_frame,
                  const int64_t last_frame,
-                 float *velocities)
+                 float **velocities,
+                 int64_t *stride_length)
 {
     int64_t n_particles, n_values_per_frame;
     tng_data_type type;
@@ -13153,8 +13201,9 @@ tng_function_status DECLSPECDLLEXPORT tng_util_vel_read_range
     stat = tng_particle_data_vector_interval_get(tng_data, TNG_TRAJ_VELOCITIES,
                                                  first_frame, last_frame,
                                                  TNG_USE_HASH,
-                                                 (void **)(&velocities),
+                                                 (void **)velocities,
                                                  &n_particles,
+                                                 stride_length,
                                                  &n_values_per_frame,
                                                  &type);
 
@@ -13165,7 +13214,8 @@ tng_function_status DECLSPECDLLEXPORT tng_util_force_read_range
                 (tng_trajectory_t tng_data,
                  const int64_t first_frame,
                  const int64_t last_frame,
-                 float *forces)
+                 float **forces,
+                 int64_t *stride_length)
 {
     int64_t n_particles, n_values_per_frame;
     tng_data_type type;
@@ -13174,8 +13224,9 @@ tng_function_status DECLSPECDLLEXPORT tng_util_force_read_range
     stat = tng_particle_data_vector_interval_get(tng_data, TNG_TRAJ_FORCES,
                                                  first_frame, last_frame,
                                                  TNG_USE_HASH,
-                                                 (void **)(&forces),
+                                                 (void **)forces,
                                                  &n_particles,
+                                                 stride_length,
                                                  &n_values_per_frame,
                                                  &type);
 
@@ -13193,7 +13244,7 @@ tng_function_status DECLSPECDLLEXPORT tng_util_pos_write
     int64_t n_particles, n_frames = 10000, stride_length = 100, frame_pos;
     tng_function_status stat;
 
-    if(!frame_set)
+    if(!frame_set || tng_data->n_trajectory_frame_sets <= 0)
     {
         stat = tng_frame_set_new(tng_data, 0, n_frames);
         if(stat != TNG_SUCCESS)
@@ -13203,7 +13254,7 @@ tng_function_status DECLSPECDLLEXPORT tng_util_pos_write
             return(stat);
         }
     }
-    if(frame_set->first_frame + n_frames < frame_nr)
+    if(frame_nr >= frame_set->first_frame + n_frames)
     {
         stat = tng_frame_set_write(tng_data, TNG_USE_HASH);
         if(stat != TNG_SUCCESS)
@@ -13213,7 +13264,7 @@ tng_function_status DECLSPECDLLEXPORT tng_util_pos_write
             return(stat);
         }
         stat = tng_frame_set_new(tng_data, frame_set->first_frame +
-                                 n_frames + 1, n_frames);
+                                 n_frames, n_frames);
         if(stat != TNG_SUCCESS)
         {
             printf("Cannot create frame set.  %s: %d\n", __FILE__,
@@ -13242,12 +13293,16 @@ tng_function_status DECLSPECDLLEXPORT tng_util_pos_write
                    __LINE__);
             return(stat);
         }
+        data = &frame_set->tr_particle_data[frame_set->
+                                            n_particle_data_blocks - 1];
         stat = tng_allocate_particle_data_mem(tng_data, data, n_frames,
                                               stride_length, n_particles,
                                               3);
     }
-    memcpy(data->values + sizeof(float) * frame_pos, positions,
-           sizeof(float) * 3 * n_particles);
+    memcpy(data->values + sizeof(float) * frame_pos * n_particles * 3,
+           positions, sizeof(float) * n_particles * 3);
+
+    frame_set->n_unwritten_frames += stride_length;
 
     return(TNG_SUCCESS);
 }
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index bee72c9..aef21ee 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -33,6 +33,13 @@ if(BUILD_TNG_EXAMPLES)
             target_link_libraries(md_openmp m)
         endif()
         set_property(TARGET md_openmp PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples)
+
+        add_executable(md_openmp_util md_openmp_util.c)
+        target_link_libraries(md_openmp_util tng_io ${OpenMP_LIBS})
+        if(UNIX)
+            target_link_libraries(md_openmp_util m)
+        endif()
+        set_property(TARGET md_openmp_util PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples)
     endif()
 
     add_executable(tng_io_read_pos tng_io_read_pos.c)
@@ -42,6 +49,13 @@ if(BUILD_TNG_EXAMPLES)
     endif()
     set_property(TARGET tng_io_read_pos PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples)
 
+    add_executable(tng_io_read_pos_util tng_io_read_pos_util.c)
+    target_link_libraries(tng_io_read_pos_util tng_io)
+    if(UNIX)
+        target_link_libraries(tng_io_read_pos_util m)
+    endif()
+    set_property(TARGET tng_io_read_pos_util PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples)
+
     add_executable(tng_parallel_read tng_parallel_read.c)
     target_link_libraries(tng_parallel_read tng_io)
     if(UNIX)
diff --git a/src/tests/md_openmp.c b/src/tests/md_openmp.c
index 58f7c6d..108beb6 100644
--- a/src/tests/md_openmp.c
+++ b/src/tests/md_openmp.c
@@ -240,7 +240,7 @@ int main ( int argc, char *argv[] )
                                 TNG_TRAJECTORY_BLOCK,
                                 n_frames_per_frame_set, 3,
                                 1, 0, np,
-                                TNG_UNCOMPRESSED,
+                                TNG_TNG_COMPRESSION,
                                 0) != TNG_SUCCESS)
     {
         printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
@@ -252,7 +252,7 @@ int main ( int argc, char *argv[] )
                                 TNG_TRAJECTORY_BLOCK,
                                 n_frames_per_frame_set, 3,
                                 1, 0, np,
-                                TNG_UNCOMPRESSED,
+                                TNG_TNG_COMPRESSION,
                                 0) != TNG_SUCCESS)
     {
         printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
diff --git a/src/tests/md_openmp_util.c b/src/tests/md_openmp_util.c
new file mode 100644
index 0000000..6a2f377
--- /dev/null
+++ b/src/tests/md_openmp_util.c
@@ -0,0 +1,707 @@
+# include <stdlib.h>
+# include <stdio.h>
+# include <time.h>
+# include <math.h>
+# include <omp.h>
+# include <tng_io.h>
+#include "tng_io_testing.h"
+
+int main ( int argc, char *argv[] );
+void compute ( int np, int nd, float pos[], float vel[],
+    float mass, float f[], float *pot, float *kin );
+float dist ( int nd, float r1[], float r2[], float dr[] );
+void initialize ( int np, int nd, float box[], int *seed, float pos[],
+    float vel[], float acc[] );
+float r8_uniform_01 ( int *seed );
+void timestamp ( void );
+void update ( int np, int nd, float pos[], float vel[], float f[],
+    float acc[], float mass, float dt );
+
+/******************************************************************************/
+
+int main ( int argc, char *argv[] )
+
+/******************************************************************************/
+/*
+    Purpose:
+
+        MAIN is the main program for MD_OPENMP.
+
+    Discussion:
+
+        MD implements a simple molecular dynamics simulation.
+
+        The program uses Open MP directives to allow parallel computation.
+
+        The velocity Verlet time integration scheme is used.
+
+        The particles interact with a central pair potential.
+
+        Output of the program is saved in the TNG format, which is why this
+        code is included in the TNG API release. The high-level API of the
+        TNG API is used where appropriate.
+
+    Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+    Modified:
+
+        8 Jan 2013
+
+    Author:
+
+        Original FORTRAN77 version by Bill Magro.
+        C version by John Burkardt.
+        TNG trajectory output by Magnus Lundborg.
+
+    Parameters:
+
+        None
+*/
+{
+    float *acc;
+    float *box;
+    float *box_shape;
+    float dt = 0.0002;
+    float e0;
+    float *force;
+    int i;
+    float kinetic;
+    float mass = 1.0;
+    int nd = 3;
+    int np = 50;
+    float *pos;
+    float potential;
+    int proc_num;
+    int seed = 123456789;
+    int step;
+    int step_num = 50000;
+    int step_print;
+    int step_print_index;
+    int step_print_num;
+    int step_save;
+    float *vel;
+    float wtime;
+    tng_trajectory_t traj;
+    tng_molecule_t molecule;
+    tng_chain_t chain;
+    tng_residue_t residue;
+    tng_atom_t atom;
+
+    timestamp ( );
+
+    proc_num = omp_get_num_procs ( );
+
+    acc = ( float * ) malloc ( nd * np * sizeof ( float ) );
+    box = ( float * ) malloc ( nd * sizeof ( float ) );
+    box_shape = (float *) malloc (9 * sizeof (float));
+    force = ( float * ) malloc ( nd * np * sizeof ( float ) );
+    pos = ( float * ) malloc ( nd * np * sizeof ( float ) );
+    vel = ( float * ) malloc ( nd * np * sizeof ( float ) );
+
+    printf ( "\n" );
+    printf ( "MD_OPENMP\n" );
+    printf ( "  C/OpenMP version\n" );
+    printf ( "\n" );
+    printf ( "  A molecular dynamics program.\n" );
+
+    printf ( "\n" );
+    printf ( "  NP, the number of particles in the simulation is %d\n", np );
+    printf ( "  STEP_NUM, the number of time steps, is %d\n", step_num );
+    printf ( "  DT, the size of each time step, is %f\n", dt );
+
+    printf ( "\n" );
+    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");
+#ifdef EXAMPLE_FILES_DIR
+    tng_util_trajectory_open(EXAMPLE_FILES_DIR "tng_md_out.tng", 'w', &traj);
+#else
+    tng_util_trajectory_open("/tmp/tng_md_out.tng", 'w', &traj);
+#endif
+
+
+
+    /* Set molecules data */
+    printf("  Creating molecules in trajectory.\n");
+    tng_molecule_add(traj, "water", &molecule);
+    tng_molecule_chain_add(traj, molecule, "W", &chain);
+    tng_chain_residue_add(traj, chain, "WAT", &residue);
+    if(tng_residue_atom_add(traj, residue, "O", "O", &atom) == TNG_CRITICAL)
+    {
+        tng_util_trajectory_close(&traj);
+        printf("  Cannot create molecules.\n");
+        exit(1);
+    }
+    tng_molecule_cnt_set(traj, molecule, np);
+
+
+/*
+    Set the dimensions of the box.
+*/
+    for(i = 0; i < 9; i++)
+    {
+        box_shape[i] = 0.0;
+    }
+    for ( i = 0; i < nd; i++ )
+    {
+        box[i] = 10.0;
+        box_shape[i*nd + i] = box[i];
+    }
+
+
+    /* Add the box shape data block and write the file headers */
+    if(tng_data_block_add(traj, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE", TNG_DOUBLE_DATA,
+                       TNG_NON_TRAJECTORY_BLOCK, 1, 9, 1, TNG_UNCOMPRESSED,
+                       box_shape) == TNG_CRITICAL ||
+                       tng_file_headers_write(traj, TNG_USE_HASH) == TNG_CRITICAL)
+    {
+        free(box_shape);
+        tng_util_trajectory_close(&traj);
+        printf("  Cannot write trajectory headers and box shape.\n");
+        exit(1);
+    }
+    free(box_shape);
+
+    printf ( "\n" );
+    printf ( "  Initializing positions, velocities, and accelerations.\n" );
+/*
+    Set initial positions, velocities, and accelerations.
+*/
+    initialize ( np, nd, box, &seed, pos, vel, acc );
+/*
+    Compute the forces and energies.
+*/
+    printf ( "\n" );
+    printf ( "  Computing initial forces and energies.\n" );
+
+    compute ( np, nd, pos, vel, mass, force, &potential, &kinetic );
+
+    e0 = potential + kinetic;
+
+    /* Saving frequency */
+    step_save = 100;
+
+    step_print = 0;
+    step_print_index = 0;
+    step_print_num = 10;
+
+/*
+    This is the main time stepping loop:
+        Compute forces and energies,
+        Update positions, velocities, accelerations.
+*/
+    printf("  Every %d steps particle positions, velocities and forces are\n",
+           step_save);
+    printf("  saved to a TNG trajectory file.\n");
+    printf ( "\n" );
+    printf ( "  At certain step intervals, we report the potential and kinetic energies.\n" );
+    printf ( "  The sum of these energies should be a constant.\n" );
+    printf ( "  As an accuracy check, we also print the relative error\n" );
+    printf ( "  in the total energy.\n" );
+    printf ( "\n" );
+    printf ( "      Step      Potential       Kinetic        (P+K-E0)/E0\n" );
+    printf ( "                Energy P        Energy K       Relative Energy Error\n" );
+    printf ( "\n" );
+
+    step = 0;
+    printf ( "  %8d  %14f  %14f  %14e\n",
+        step, potential, kinetic, ( potential + kinetic - e0 ) / e0 );
+    step_print_index++;
+    step_print = ( step_print_index * step_num ) / step_print_num;
+
+    if(tng_util_pos_write(traj, 0, pos) != TNG_SUCCESS)
+    {
+        printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+        exit(1);
+    }
+
+    wtime = omp_get_wtime ( );
+
+    for ( step = 1; step <= step_num; step++ )
+    {
+        compute ( np, nd, pos, vel, mass, force, &potential, &kinetic );
+
+        if ( step == step_print )
+        {
+            printf ( "  %8d  %14f  %14f  %14e\n", step, potential, kinetic,
+             ( potential + kinetic - e0 ) / e0 );
+            step_print_index++;
+            step_print = ( step_print_index * step_num ) / step_print_num;
+        }
+        if(step % step_save == 0)
+        {
+            if(tng_util_pos_write(traj, step, pos) != TNG_SUCCESS)
+            {
+                printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+                exit(1);
+            }
+        }
+        update ( np, nd, pos, vel, force, acc, mass, dt );
+    }
+    wtime = omp_get_wtime ( ) - wtime;
+
+    printf ( "\n" );
+    printf ( "  Elapsed time for main computation:\n" );
+    printf ( "  %f seconds.\n", wtime );
+
+    free ( acc );
+    free ( box );
+    free ( force );
+    free ( pos );
+    free ( vel );
+/*
+    Terminate.
+*/
+    tng_util_trajectory_close(&traj);
+
+    printf ( "\n" );
+    printf ( "MD_OPENMP\n" );
+    printf ( "  Normal end of execution.\n" );
+
+    printf ( "\n" );
+    timestamp ( );
+
+    return 0;
+}
+/******************************************************************************/
+
+void compute ( int np, int nd, float pos[], float vel[],
+    float mass, float f[], float *pot, float *kin )
+
+/******************************************************************************/
+/*
+    Purpose:
+
+        COMPUTE computes the forces and energies.
+
+    Discussion:
+
+        The computation of forces and energies is fully parallel.
+
+        The potential function V(X) is a harmonic well which smoothly
+        saturates to a maximum value at PI/2:
+
+            v(x) = ( sin ( min ( x, PI2 ) ) )**2
+
+        The derivative of the potential is:
+
+            dv(x) = 2.0 * sin ( min ( x, PI2 ) ) * cos ( min ( x, PI2 ) )
+                        = sin ( 2.0 * min ( x, PI2 ) )
+
+    Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+    Modified:
+
+        21 November 2007
+
+    Author:
+
+        Original FORTRAN77 version by Bill Magro.
+        C version by John Burkardt.
+
+    Parameters:
+
+        Input, int NP, the number of particles.
+
+        Input, int ND, the number of spatial dimensions.
+
+        Input, float POS[ND*NP], the position of each particle.
+
+        Input, float VEL[ND*NP], the velocity of each particle.
+
+        Input, float MASS, the mass of each particle.
+
+        Output, float F[ND*NP], the forces.
+
+        Output, float *POT, the total potential energy.
+
+        Output, float *KIN, the total kinetic energy.
+*/
+{
+    float d;
+    float d2;
+    int i;
+    int j;
+    int k;
+    float ke;
+    float pe;
+    float PI2 = 3.141592653589793 / 2.0;
+    float rij[3];
+
+    pe = 0.0;
+    ke = 0.0;
+
+# 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++ )
+    {
+/*
+    Compute the potential energy and forces.
+*/
+        for ( i = 0; i < nd; i++ )
+        {
+            f[i+k*nd] = 0.0;
+        }
+
+        for ( j = 0; j < np; j++ )
+        {
+            if ( k != j )
+            {
+                d = dist ( nd, pos+k*nd, pos+j*nd, rij );
+/*
+    Attribute half of the potential energy to particle J.
+*/
+                if ( d < PI2 )
+                {
+                    d2 = d;
+                }
+                else
+                {
+                    d2 = PI2;
+                }
+
+                pe = pe + 0.5 * pow ( sin ( d2 ), 2 );
+
+                for ( i = 0; i < nd; i++ )
+                {
+                    f[i+k*nd] = f[i+k*nd] - rij[i] * sin ( 2.0 * d2 ) / d;
+                }
+            }
+        }
+/*
+    Compute the kinetic energy.
+*/
+        for ( i = 0; i < nd; i++ )
+        {
+            ke = ke + vel[i+k*nd] * vel[i+k*nd];
+        }
+    }
+
+    ke = ke * 0.5 * mass;
+
+    *pot = pe;
+    *kin = ke;
+
+    return;
+}
+/******************************************************************************/
+
+float dist ( int nd, float r1[], float r2[], float dr[] )
+
+/******************************************************************************/
+/*
+    Purpose:
+
+        DIST computes the displacement (and its norm) between two particles.
+
+    Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+    Modified:
+
+        21 November 2007
+
+    Author:
+
+        Original FORTRAN77 version by Bill Magro.
+        C version by John Burkardt.
+
+    Parameters:
+
+        Input, int ND, the number of spatial dimensions.
+
+        Input, float R1[ND], R2[ND], the positions of the particles.
+
+        Output, float DR[ND], the displacement vector.
+
+        Output, float D, the Euclidean norm of the displacement.
+*/
+{
+    float d;
+    int i;
+
+    d = 0.0;
+    for ( i = 0; i < nd; i++ )
+    {
+        dr[i] = r1[i] - r2[i];
+        d = d + dr[i] * dr[i];
+    }
+    d = sqrt ( d );
+
+    return d;
+}
+/******************************************************************************/
+
+void initialize ( int np, int nd, float box[], int *seed, float pos[],
+    float vel[], float acc[] )
+
+/******************************************************************************/
+/*
+    Purpose:
+
+        INITIALIZE initializes the positions, velocities, and accelerations.
+
+    Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+    Modified:
+
+        21 November 2007
+
+    Author:
+
+        Original FORTRAN77 version by Bill Magro.
+        C version by John Burkardt.
+
+    Parameters:
+
+        Input, int NP, the number of particles.
+
+        Input, int ND, the number of spatial dimensions.
+
+        Input, float BOX[ND], specifies the maximum position
+        of particles in each dimension.
+
+        Input, int *SEED, a seed for the random number generator.
+
+        Output, float POS[ND*NP], the position of each particle.
+
+        Output, float VEL[ND*NP], the velocity of each particle.
+
+        Output, float ACC[ND*NP], the acceleration of each particle.
+*/
+{
+    int i;
+    int j;
+/*
+    Give the particles random positions within the box.
+*/
+    for ( i = 0; i < nd; i++ )
+    {
+        for ( j = 0; j < np; j++ )
+        {
+            pos[i+j*nd] = box[i] * r8_uniform_01 ( seed );
+        }
+    }
+
+    for ( j = 0; j < np; j++ )
+    {
+        for ( i = 0; i < nd; i++ )
+        {
+            vel[i+j*nd] = 0.0;
+        }
+    }
+    for ( j = 0; j < np; j++ )
+    {
+        for ( i = 0; i < nd; i++ )
+        {
+            acc[i+j*nd] = 0.0;
+        }
+    }
+    return;
+}
+/******************************************************************************/
+
+float r8_uniform_01 ( int *seed )
+
+/******************************************************************************/
+/*
+    Purpose:
+
+        R8_UNIFORM_01 is a unit pseudorandom R8.
+
+    Discussion:
+
+        This routine implements the recursion
+
+            seed = 16807 * seed mod ( 2**31 - 1 )
+            unif = seed / ( 2**31 - 1 )
+
+        The integer arithmetic never requires more than 32 bits,
+        including a sign bit.
+
+    Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+    Modified:
+
+        11 August 2004
+
+    Author:
+
+        John Burkardt
+
+    Reference:
+
+        Paul Bratley, Bennett Fox, Linus Schrage,
+        A Guide to Simulation,
+        Springer Verlag, pages 201-202, 1983.
+
+        Bennett Fox,
+        Algorithm 647:
+        Implementation and Relative Efficiency of Quasirandom
+        Sequence Generators,
+        ACM Transactions on Mathematical Software,
+        Volume 12, Number 4, pages 362-376, 1986.
+
+    Parameters:
+
+        Input/output, int *SEED, a seed for the random number generator.
+
+        Output, float R8_UNIFORM_01, a new pseudorandom variate, strictly between
+        0 and 1.
+*/
+{
+    int k;
+    float r;
+
+    k = *seed / 127773;
+
+    *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
+
+    if ( *seed < 0 )
+    {
+        *seed = *seed + 2147483647;
+    }
+
+    r = ( float ) ( *seed ) * 4.656612875E-10;
+
+    return r;
+}
+/******************************************************************************/
+
+void timestamp ( void )
+
+/******************************************************************************/
+/*
+    Purpose:
+
+        TIMESTAMP prints the current YMDHMS date as a time stamp.
+
+    Example:
+
+        31 May 2001 09:45:54 AM
+
+    Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+    Modified:
+
+        24 September 2003
+
+    Author:
+
+        John Burkardt
+
+    Parameters:
+
+        None
+*/
+{
+# define TIME_SIZE 40
+
+    static char time_buffer[TIME_SIZE];
+    const struct tm *tm;
+    time_t now;
+
+    now = time ( NULL );
+    tm = localtime ( &now );
+
+    strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
+
+    printf ( "%s\n", time_buffer );
+
+    return;
+# undef TIME_SIZE
+}
+/******************************************************************************/
+
+void update ( int np, int nd, float pos[], float vel[], float f[],
+    float acc[], float mass, float dt )
+
+/******************************************************************************/
+/*
+    Purpose:
+
+        UPDATE updates positions, velocities and accelerations.
+
+    Discussion:
+
+        The time integration is fully parallel.
+
+        A velocity Verlet algorithm is used for the updating.
+
+        x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt * dt
+        v(t+dt) = v(t) + 0.5 * ( a(t) + a(t+dt) ) * dt
+        a(t+dt) = f(t) / m
+
+    Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+    Modified:
+
+        17 April 2009
+
+    Author:
+
+        Original FORTRAN77 version by Bill Magro.
+        C version by John Burkardt.
+
+    Parameters:
+
+        Input, int NP, the number of particles.
+
+        Input, int ND, the number of spatial dimensions.
+
+        Input/output, float POS[ND*NP], the position of each particle.
+
+        Input/output, float VEL[ND*NP], the velocity of each particle.
+
+        Input, float F[ND*NP], the force on each particle.
+
+        Input/output, float ACC[ND*NP], the acceleration of each particle.
+
+        Input, float MASS, the mass of each particle.
+
+        Input, float DT, the time step.
+*/
+{
+    int i;
+    int j;
+    float rmass;
+
+    rmass = 1.0 / mass;
+
+# pragma omp parallel \
+    shared ( acc, dt, f, nd, np, pos, rmass, vel ) \
+    private ( i, j )
+
+# pragma omp for
+    for ( j = 0; j < np; j++ )
+    {
+        for ( i = 0; i < nd; i++ )
+        {
+            pos[i+j*nd] = pos[i+j*nd] + vel[i+j*nd] * dt + 0.5 * acc[i+j*nd] * dt * dt;
+            vel[i+j*nd] = vel[i+j*nd] + 0.5 * dt * ( f[i+j*nd] * rmass + acc[i+j*nd] );
+            acc[i+j*nd] = f[i+j*nd] * rmass;
+        }
+    }
+
+    return;
+}
diff --git a/src/tests/tng_io_read_pos_util.c b/src/tests/tng_io_read_pos_util.c
new file mode 100644
index 0000000..12bf2ed
--- /dev/null
+++ b/src/tests/tng_io_read_pos_util.c
@@ -0,0 +1,113 @@
+/* This code is part of the tng binary trajectory format.
+ *
+ *  The high-level API of the TNG API is used where appropriate.
+ *
+ *                      VERSION 1.0
+ *
+ * Written by Magnus Lundborg
+ * Copyright (c) 2012, The GROMACS development team.
+ * check out http://www.gromacs.org for more information.
+ *
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <tng_io.h>
+
+int main(int argc, char **argv)
+{
+    tng_trajectory_t traj;
+    float *positions = 0; // A 1-dimensional array
+                                      // to be populated
+    int64_t n_particles, n_frames, tot_n_frames, stride_length, i, j;
+    // Set a default frame range
+    int64_t first_frame = 0, last_frame = 5000;
+    int k;
+
+    if(argc <= 1)
+    {
+        printf("No file specified\n");
+        printf("Usage:\n");
+        printf("tng_io_read_pos <tng_file> "
+               "[first_frame = %"PRId64"] [last_frame = %"PRId64"]\n",
+               first_frame, last_frame);
+        exit(1);
+    }
+
+    // A reference must be passed to allocate memory
+    tng_util_trajectory_open(argv[1], 'r', &traj);
+    
+    if(argc >= 3)
+    {
+        first_frame = strtoll(argv[2], 0, 10);
+        if(argc >= 4)
+        {
+            last_frame = strtoll(argv[3], 0, 10);
+        }
+    }
+
+    if(tng_num_frames_get(traj, &tot_n_frames) != TNG_SUCCESS)
+    {
+        printf("Cannot determine the number of frames in the file\n");
+        tng_util_trajectory_close(&traj);
+        exit(1);
+    }
+
+    if(tng_num_particles_get(traj, &n_particles) != TNG_SUCCESS)
+    {
+        printf("Cannot determine the number of particles in the file\n");
+        tng_util_trajectory_close(&traj);
+        exit(1);
+    }
+    
+    printf("%"PRId64" frames in file\n", tot_n_frames);
+
+    if(last_frame > tot_n_frames - 1)
+    {
+        last_frame = tot_n_frames - 1;
+    }
+
+    n_frames = last_frame - first_frame + 1;
+
+
+    // Get the positions of all particles in the requested frame range.
+    // The positions are stored in the positions array.
+    // N.B. No proper error checks.
+    if(tng_util_pos_read_range(traj, 0, last_frame, &positions, &stride_length)
+       == TNG_SUCCESS)
+    {
+        // Print the positions of the wanted particle (zero based)
+        for(i=0; i < n_frames; i += stride_length)
+        {
+            printf("\nFrame %"PRId64":\n", first_frame + i);
+            for(j=0; j < n_particles; j++)
+            {
+                printf("Atom nr: %"PRId64"", j);
+                for(k=0; k < 3; k++)
+                {
+                    printf("\t%f", positions[i/stride_length*n_particles*
+                                             3+j*3+k]);
+                }
+                printf("\n");
+            }
+        }
+    }
+    else
+    {
+        printf("Cannot read positions\n");
+    }
+
+    // Free memory
+    if(positions)
+    {
+        free(positions);
+    }
+    tng_util_trajectory_close(&traj);
+
+    return(0);
+}