Additional example using a very light md implementation using openmp. Data is saved during simulation.

7 files changed, 2361 insertions, 4 deletions

diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt
index b12bb42..dda277d 100644
--- a/src/lib/CMakeLists.txt
+++ b/src/lib/CMakeLists.txt
@@ -1,2 +1,7 @@
+if(NOT OPENMP_FOUND)
+  add_library(openmp SHARED openmp_stubs.c)
+  set(OpenMP_LIBS "openmp")
+endif()
+
 add_library(tng_io SHARED tng_io.c md5.c)
 
diff --git a/src/lib/omp.h b/src/lib/omp.h
new file mode 100644
index 0000000..179797f
--- /dev/null
+++ b/src/lib/omp.h
@@ -0,0 +1,63 @@
+#ifndef _OMP_H_DEF
+#define _OMP_H_DEF
+/*
+  Define the lock data types
+*/
+typedef void *omp_lock_t;
+typedef void *omp_nest_lock_t;
+/*
+  Define the schedule kinds.
+  Add vendor specific schedule constants at end.
+*/
+typedef enum omp_sched_t
+{
+  omp_sched_static = 1,
+  omp_sched_dynamic = 2,
+  omp_sched_guided = 3,
+  omp_sched_auto = 4
+} omp_sched_t;
+/*
+  Exported OpenMP functions
+*/
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+extern void   omp_destroy_lock ( omp_lock_t *lock );
+extern void   omp_destroy_nest_lock ( omp_nest_lock_t *lock );
+extern int    omp_get_active_level ( void );
+extern int    omp_get_ancestor_thread_num ( int level );
+extern int    omp_get_dynamic ( void );
+extern int    omp_get_level ( void );
+extern int    omp_get_max_active_levels ( void );
+extern int    omp_get_max_threads ( void );
+extern int    omp_get_nested ( void );
+extern int    omp_get_num_procs ( void );
+extern int    omp_get_num_threads ( void );
+extern void   omp_get_schedule ( omp_sched_t *kind, int *modifier );
+extern int    omp_get_team_size ( int level );
+extern int    omp_get_thread_limit ( void );
+extern int    omp_get_thread_num ( void );
+extern double omp_get_wtick ( void );
+extern double omp_get_wtime ( void );
+extern int    omp_in_final ( void );
+extern int    omp_in_parallel ( void );
+extern void   omp_init_lock ( omp_lock_t *lock );
+extern void   omp_init_nest_lock ( omp_nest_lock_t *lock );
+extern void   omp_set_dynamic ( int dynamic_threads );
+extern void   omp_set_lock ( omp_lock_t *lock );
+extern void   omp_set_max_active_levels ( int max_active_levels );
+extern void   omp_set_nest_lock ( omp_nest_lock_t *lock );
+extern void   omp_set_nested ( int nested );
+extern void   omp_set_num_threads ( int num_threads );
+extern void   omp_set_schedule ( omp_sched_t kind, int modifier );
+extern int    omp_test_lock ( omp_lock_t *lock );
+extern int    omp_test_nest_lock ( omp_nest_lock_t *lock );
+extern void   omp_unset_lock ( omp_lock_t *lock );
+extern void   omp_unset_nest_lock ( omp_nest_lock_t *lock );
+
+#ifdef __cplusplus
+}
+#endif
+#endif
\ No newline at end of file
diff --git a/src/lib/openmp_stubs.c b/src/lib/openmp_stubs.c
new file mode 100644
index 0000000..5ce31dc
--- /dev/null
+++ b/src/lib/openmp_stubs.c
@@ -0,0 +1,1402 @@
+# include <stdio.h>
+# include <stdlib.h>
+# include <time.h>
+
+# include "omp.h"
+
+struct __omp_lock
+{
+  int lock;
+};
+
+enum { UNLOCKED = -1, INIT, LOCKED };
+
+struct __omp_nest_lock
+{
+  short owner;
+  short count;
+};
+
+enum { NOOWNER = -1, MASTER = 0 };
+
+/******************************************************************************/
+
+void omp_destroy_lock ( omp_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_DESTROY_LOCK destroys a simple lock.
+
+  Discussion:
+
+    The routine is intended to return the state of the lock to the 
+    uninitialized state.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, omp_lock_t *ARG, the simple lock.
+*/
+{
+  struct __omp_lock *lock = ( struct __omp_lock * ) arg;
+  lock->lock = INIT;
+
+  return;
+}
+/******************************************************************************/
+
+void omp_destroy_nest_lock ( omp_nest_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_DESTROY_NEST_LOCK destroys a nestable lock.
+
+  Discussion:
+
+    The routine is intended to return the state of the lock to the 
+    uninitialized state.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, omp_nest_lock_t *ARG, the nestable lock.
+*/
+{
+  struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg;
+  nlock->owner = NOOWNER;
+  nlock->count = UNLOCKED;
+  return;
+}
+/******************************************************************************/
+
+int omp_get_active_level ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_ACTIVE_LEVEL returns the number of nested active parallel regions.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_ACTIVE_LEVEL, the number of nested active parallel 
+    regions enclosing the task that contains this call.
+*/
+{
+  return 0;
+}
+/******************************************************************************/
+
+int omp_get_ancestor_thread_num ( int level )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_ANCESTOR_THREAD_NUM returns the thread number of the ancestor of this thread.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input, int LEVEL, the nested level.
+
+    Output, int OMP_GET_ANCESTOR_THREAD_NUM, the thread number of the 
+    ancestor of this thread.
+*/
+{
+  if ( level == 0 )
+  {
+    return 0;
+  }
+  else
+  {
+    return -1;
+  }
+}
+/******************************************************************************/
+
+int omp_get_dynamic ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_DYNAMIC reports if dynamic adjustment of thread number is allowed.
+
+  Discussion:
+
+    The user can request dynamic thread adjustment by calling OMP_SET_DYNAMIC.
+
+    For this stub library, the value FALSE is always returned.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_DYNAMIC, is TRUE (1) if dynamic adjustment of thread
+    number has been enable, by default or by a user call, and FALSE (0)
+    otherwise.
+*/
+{
+  return 0;
+}
+/******************************************************************************/
+
+int omp_get_level ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_LEVEL returns the number of nested parallel regions enclosing this task.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_LEVEL, the number of nested parallel regions 
+    enclosing this task.
+*/
+{
+  return 0;
+}
+/******************************************************************************/
+
+int omp_get_max_active_levels ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+     OMP_GET_MAX_ACTIVE_LEVELS gets the maximum number of nested active parallel regions.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_MAX_ACTIVE_LEVELS gets the maximum number of 
+    nested active parallel regions.
+*/
+{
+  return 0;
+}
+/******************************************************************************/
+
+int omp_get_max_threads ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_MAX_THREADS returns the default number of threads.
+
+  Discussion:
+
+    If a parallel region is reached, and no number of threads has been
+    specified explicitly, there is a default number of threads that will
+    be used to form the new team.  That value is returned by this function.
+
+    For this stub library, the value 1 is always returned.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_MAX_THREADS, the default number of threads.
+*/
+{
+  return 1;
+}
+/******************************************************************************/
+
+int omp_get_nested ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_NESTED reports if nested parallelism has been enabled.
+
+  Discussion:
+
+    The user can request nested parallelism by calling OMP_SET_NESTED.
+
+    For this stub library, the value FALSE is always returned.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_NESTED, is TRUE (1) if nested parallelism has been
+    enable by default or by a user call, and FALSE (0) otherwise.
+*/
+{
+  return 0;
+}
+/******************************************************************************/
+
+int omp_get_num_procs ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_NUM_PROCS returns the number of processors available to the program.
+
+  Discussion:
+
+    For this stub library, the value 1 is always returned.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_NUM_PROCS, the number of processors available.
+*/
+{
+  return 1;
+}
+/******************************************************************************/
+
+int omp_get_num_threads ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_NUM_THREADS returns the number of threads in the current team.
+
+  Discussion:
+
+    For this stub library, the value 1 is always returned.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_NUM_THREADS, the number of threads in the 
+    current team.
+*/
+{
+  return 1;
+}
+/******************************************************************************/
+
+void omp_get_schedule ( omp_sched_t *kind, int *modifier )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+     OMP_GET_SCHEDULE returns information about the "runtime" schedule.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, omp_sched_t *KIND, may be
+    1, omp_sched_static,
+    2, omp_sched_dynamic,
+    3, omp_sched_guided,
+    4, omp_sched_auto.
+
+    Output, int *MODIFIER; this contains the "chunk_size" information for
+    static, dynamic, or guided schedules, and is ignored for the auto schedule.
+    If the chunk_size is less than 1, then the default value is used instead.
+*/
+{
+  *kind = omp_sched_static;
+  *modifier = 0;
+  return;
+}
+/******************************************************************************/
+
+int omp_get_team_size ( int level )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+     OMP_GET_TEAM_SIZE returns the size of the thread team for a given level.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input, int LEVEL, the nested level.
+
+    Output, int OMP_GET_TEAM_SIZE, the size of the thread team for 
+    this level.
+*/
+{
+  if ( level == 0 )
+  {
+    return 1;
+  }
+  else
+  {
+    return -1;
+  }
+}
+/******************************************************************************/
+
+int omp_get_thread_limit ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+     OMP_GET_THREAD_LIMIT returns the maximum number of OpenMP threads available.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_THREAD_LIMIT, the maximum number of OpenMP
+    threads available.
+*/
+{
+  return 1;
+}
+/******************************************************************************/
+
+int omp_get_thread_num ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_THREAD_NUM returns the thread number of a thread in a team.
+
+  Discussion:
+
+    Thread numbers start at 0.
+
+    If this function is not called from a parallel region, then only one
+    thread is executing, so the value returned is 0.
+
+    For this stub library, the value 0 is always returned.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_GET_THREAD_NUM, the thread number.
+*/
+{
+  return 0;
+}
+/******************************************************************************/
+
+double omp_get_wtick ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_WTICK returns the precision of the timer used by OMP_GET_WTIME.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, double OMP_GET_WTICK, the number of seconds between
+    successive "ticks" of the wall clock timer.
+*/
+{
+  double value;
+
+  value = 1.0 / ( double ) CLOCKS_PER_SEC;
+
+  return value;
+}
+/******************************************************************************/
+
+double omp_get_wtime ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_GET_WTIME returns elapsed wall clock time in seconds.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, double OMP_GET_WTIME, the current reading of the
+    wall clock timer.
+*/
+{
+  double value;
+
+  value = ( double ) clock ( ) / ( double ) CLOCKS_PER_SEC;
+
+  return value;
+}
+/******************************************************************************/
+
+int omp_in_final ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_IN_FINAL is true if the routine is executed in a final task region.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_IN_FINAL, is true if the routine is executed in a
+    final task region.
+*/
+{
+  return 1;
+}
+/******************************************************************************/
+
+int omp_in_parallel ( void )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_IN_PARALLEL returns TRUE if the call is made from a parallel region.
+
+  Discussion:
+
+    For this stub library, the value FALSE is always returned.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, int OMP_IN_PARALLEL, is "TRUE" (1) if the routine was called
+    from a parallel region and "FALSE" (0) otherwise.
+*/
+{
+  return 0;
+}
+/******************************************************************************/
+
+void omp_init_lock ( omp_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_INIT_LOCK initializes a simple lock.
+
+  Discussion:
+
+    This routine is intended to initialize the lock to the unlocked state.
+
+    For this stub library, the lock points to -1.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, omp_lock_t *ARG, the simple lock.
+*/
+{
+  struct __omp_lock *lock = ( struct __omp_lock * ) arg;
+  lock->lock = UNLOCKED;
+  return;
+}
+/******************************************************************************/
+
+void omp_init_nest_lock ( omp_nest_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_INIT_NEST_LOCK initializes a nestable lock.
+
+  Discussion:
+
+    This routine is intended to initialize the lock to the unlocked state.
+
+    For this stub library, the lock points to -1.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Output, omp_nest_lock_t *ARG, the lock.
+*/
+{
+  struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg;
+  nlock->owner = NOOWNER;
+  nlock->count = 0;
+  return;
+}
+/******************************************************************************/
+
+void omp_set_dynamic ( int dynamic_threads )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_SET_DYNAMIC enables dynamic adjustment of the number of threads.
+
+  Discussion:
+
+    If DYNAMIC_THREADS is TRUE, then the number of threads available for
+    execution in a parallel region may be dynamically adjusted.
+
+    For this stub library, the input value is ignored.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input, int DYNAMIC_THREADS, is TRUE if the user wishes to allow
+    dynamic adjustment of the number of threads available for execution
+    in any parallel region.
+*/
+{
+  return;
+}
+/******************************************************************************/
+
+void omp_set_lock ( omp_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_SET_LOCK sets a simple lock.
+
+  Discussion:
+
+    The lock must already have been initialized.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input/output, omp_lock_t *ARG, the simple lock.
+*/
+{
+  struct __omp_lock *lock = ( struct __omp_lock * ) arg;
+  if ( lock->lock == UNLOCKED )
+  {
+    lock->lock = LOCKED;
+  }
+  else if ( lock->lock == LOCKED )
+  {
+    fprintf ( stderr, "error: deadlock in using lock variable\n" );
+    exit ( 1 );
+  }
+  else
+  {
+    fprintf ( stderr, "error: lock not initialized\n" );
+    exit ( 1 );
+  }
+  return;
+}
+/******************************************************************************/
+
+void omp_set_max_active_levels ( int max_active_levels )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_SET_MAX_ACTIVE_LEVELS limits the number of nested active parallel regions.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input, int MAX_LEVELS, the maximum number of nested active parallel
+    regions.
+*/
+{
+  return;
+}
+/******************************************************************************/
+
+void omp_set_nest_lock ( omp_nest_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_SET_NEST_LOCK sets a nestable lock.
+
+  Discussion:
+
+    The lock must already have been initialized.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input/output, omp_nest_lock_t *ARG, the nestable lock.
+*/
+{
+  struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg;
+  if ( nlock->owner == MASTER && nlock->count >= 1 )
+  {
+    nlock->count++;
+  }
+  else if ( nlock->owner == NOOWNER && nlock->count == 0 )
+  {
+    nlock->owner = MASTER;
+    nlock->count = 1;
+  }
+  else
+  {
+    fprintf ( stderr, "error: lock corrupted or not initialized\n" );
+    exit ( 1 );
+  }
+  return;
+}
+/******************************************************************************/
+
+void omp_set_nested ( int nested )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_SET_NESTED controls the use of nested parallelism.
+
+  Discussion:
+
+    For this stub library, the input value is ignored.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input, int NESTED, is TRUE (1) if nested parallelism is to be enabled.
+*/
+{
+  return;
+}
+/******************************************************************************/
+
+void omp_set_num_threads ( int num_threads )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_SET_NUM_THREADS sets the number of threads.
+
+  Discussion:
+
+    This routine sets the number of threads to be used in all subsequent
+    parallel regions for which an explicit number of threads is not
+    specified.
+
+    For this stub library, the input value is ignored.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input, int NUM_THREADS, the number of threads to be used in all
+    subsequent parallel regions for which an explicit number of threads
+    is not specified.  0 < NUM_THREADS.
+*/
+{
+  return;
+}
+/******************************************************************************/
+
+void omp_set_schedule ( omp_sched_t kind, int modifier )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_SET_SCHEDULE chooses the schedule when "runtime" is the schedule kind.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input, omp_sched_t KIND, may be
+    1, omp_sched_static,
+    2, omp_sched_dynamic,
+    3, omp_sched_guided,
+    4, omp_sched_auto.
+
+    Input, int MODIFIER; this contains the "chunk_size" information for
+    static, dynamic, or guided schedules, and is ignored for the auto schedule.
+    If the chunk_size is less than 1, then the default value is used instead.
+*/
+{
+  return;
+}
+/******************************************************************************/
+
+int omp_test_lock ( omp_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_TEST_LOCK tests a simple lock.
+
+  Discussion:
+
+    Calling this routine with an uninitialized lock causes an error.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input/output, omp_lock_t *ARG, the simple lock.
+    If the lock was initialized but not set, it is set by this call.
+
+    Output, int OMP_TEST_LOCK:
+    TRUE (1), on input, the lock was initialized and not set;
+    FALSE (0), on input the lock was initialized, and set.
+*/
+{
+  struct __omp_lock *lock = ( struct __omp_lock * ) arg;
+  if ( lock->lock == UNLOCKED )
+  {
+    lock->lock = LOCKED;
+    return 1;
+  }
+  else if ( lock->lock == LOCKED )
+  {
+    return 0;
+  }
+  else
+  {
+    fprintf ( stderr, "error: lock not initialized\n" );
+    exit ( 1 );
+  }
+}
+/******************************************************************************/
+
+int omp_test_nest_lock ( omp_nest_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_TEST_NEST_LOCK tests a nestable lock.
+
+  Discussion:
+
+    Calling this routine with an uninitialized lock causes an error.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input/output, omp_nest_lock_t *ARG, the nestable lock.
+    If the lock was initialized but not set, it is set by this call.
+
+    Output, int OMP_TEST_NEST_LOCK, returns the new nesting count,
+    if the call was successful.  Otherwise, the value 0 is returned.
+*/
+{
+  struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg;
+  omp_set_nest_lock ( arg) ;
+  return nlock->count;
+}
+/******************************************************************************/
+
+void omp_unset_lock ( omp_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_UNSET_LOCK unsets a simple lock.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input/output, omp_lock_t *ARG, the simple lock.
+*/
+{
+  struct __omp_lock *lock = ( struct __omp_lock * ) arg;
+  if ( lock->lock == LOCKED )
+  {
+    lock->lock = UNLOCKED;
+  }
+  else if ( lock->lock == UNLOCKED )
+  {
+    fprintf ( stderr, "error: lock not set\n" );
+    exit ( 1 );
+  }
+  else
+  {
+    fprintf ( stderr, "error: lock not initialized\n" );
+    exit ( 1 );
+  }
+  return;
+}
+/******************************************************************************/
+
+void omp_unset_nest_lock ( omp_nest_lock_t *arg )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    OMP_UNSET_NEST_LOCK unsets a nestable lock.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    17 May 2012
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    OpenMP Application Program Interface,
+    Version 3.1,
+    July 2011.
+
+  Parameters:
+
+    Input/output, omp_nest_lock_t *ARG, the nestable lock.
+*/
+{
+  struct __omp_nest_lock *nlock = ( struct __omp_nest_lock * ) arg;
+  if ( nlock->owner == MASTER && nlock->count >= 1 )
+  {
+    nlock->count--;
+    if ( nlock->count == 0 )
+    {
+      nlock->owner = NOOWNER;
+    }
+  }
+  else if ( nlock->owner == NOOWNER && nlock->count == 0 )
+  {
+    fprintf ( stderr, "error: lock not set\n" );
+    exit ( 1 );
+  }
+  else
+  {
+    fprintf ( stderr, "error: lock corrupted or not initialized\n" );
+    exit ( 1 );
+  }
+  return;
+}
diff --git a/src/lib/tng_io.c b/src/lib/tng_io.c
index c97a3fb..12b7a4e 100644
--- a/src/lib/tng_io.c
+++ b/src/lib/tng_io.c
@@ -8275,7 +8275,7 @@ tng_function_status tng_frame_particle_data_write(tng_trajectory_t tng_data,
         memcpy(copy, values, val_n_particles * n_values_per_frame * size);
         for(i = 0; i < val_n_particles * n_values_per_frame; i++)
         {
-            tng_swap_byte_order_64(tng_data, &copy[i*size]);
+            tng_swap_byte_order_64(tng_data, (int64_t *)copy+i*size);
         }
         fwrite(copy, val_n_particles * n_values_per_frame, size,
                tng_data->output_file);
@@ -8288,7 +8288,7 @@ tng_function_status tng_frame_particle_data_write(tng_trajectory_t tng_data,
         memcpy(copy, values, val_n_particles * n_values_per_frame * size);
         for(i = 0; i < val_n_particles * n_values_per_frame; i++)
         {
-            tng_swap_byte_order_32(tng_data, &copy[i*size]);
+            tng_swap_byte_order_32(tng_data, (int32_t *)copy+i*size);
         }
         fwrite(copy, val_n_particles * n_values_per_frame, size,
                tng_data->output_file);
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index 6747eb6..92c580d 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -1,3 +1,7 @@
+if(NOT OPENMP_FOUND)
+  set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+endif()
+
 include_directories(${TRAJECTORY_SOURCE_DIR}/src/lib)
 link_directories(${TRAJECTORY_BINARY_DIR}/src/lib)
 
@@ -7,4 +11,7 @@ configure_file(tng_io_testing.h.in ${CMAKE_BINARY_DIR}/generated/tng_io_testing.
 include_directories(${CMAKE_BINARY_DIR}/generated/)
 
 add_executable(tng_testing tng_io_testing.c)
-target_link_libraries(tng_testing tng_io)
\ No newline at end of file
+target_link_libraries(tng_testing tng_io)
+
+add_executable(md_openmp md_openmp.c)
+target_link_libraries(md_openmp tng_io ${OpenMP_LIBS} m)
\ No newline at end of file
diff --git a/src/tests/md_openmp.c b/src/tests/md_openmp.c
new file mode 100644
index 0000000..0b5241a
--- /dev/null
+++ b/src/tests/md_openmp.c
@@ -0,0 +1,880 @@
+# 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, 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[], 
+    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[], 
+    double acc[], double mass, double 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.
+
+    Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+    Modified:
+
+        30 July 2009
+
+    Author:
+
+        Original FORTRAN77 version by Bill Magro.
+        C version by John Burkardt.
+        TNG trajectory output by Magnus Lundborg.
+
+    Parameters:
+
+        None
+*/
+{
+    double *acc;
+    double *box;
+    double *box_shape;
+    double dt = 0.0001;
+    double e0;
+    double *force;
+    int i;
+    double kinetic;
+    double mass = 1.0;
+    int nd = 3;
+    int np = 100;
+    double *pos;
+    double potential;
+    int proc_num;
+    int seed = 123456789;
+    int step;
+    int step_num = 1000;
+    int step_print;
+    int step_print_index;
+    int step_print_num;
+    int step_save;
+    double *vel;
+    double wtime;
+    tng_trajectory_t traj;
+    tng_molecule_t molecule;
+    tng_chain_t chain;
+    tng_residue_t residue;
+    tng_atom_t atom;
+    int64_t n_frames_per_frame_set;
+    int frames_saved_cnt = 0;
+    int frame_set_cnt = 0;
+    double **data;
+
+    timestamp ( );
+
+    proc_num = omp_get_num_procs ( );
+
+    acc = ( double * ) malloc ( nd * np * sizeof ( double ) );
+    box = ( double * ) malloc ( nd * sizeof ( double ) );
+    box_shape = (double *) malloc (9 * sizeof (double));
+    force = ( double * ) malloc ( nd * np * sizeof ( double ) );
+    pos = ( double * ) malloc ( nd * np * sizeof ( double ) );
+    vel = ( double * ) malloc ( nd * np * sizeof ( double ) );
+
+    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", omp_get_num_procs ( ) );
+    printf ( "  Number of threads =              %d\n", omp_get_max_threads ( ) );
+
+    
+    printf("\n");
+    printf("  Initializing trajectory storage.\n");
+    if(tng_trajectory_init(&traj) != TNG_SUCCESS)
+    {
+        tng_trajectory_destroy(&traj);
+        printf("  Cannot init trajectory.\n");
+        exit(1);
+    }
+#ifdef EXAMPLE_FILES_DIR
+    tng_output_file_set(traj, EXAMPLE_FILES_DIR "tng_md_out.tng");
+#else
+    tng_output_file_set(traj, "/tmp/tng_md_out.tng");
+#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_trajectory_destroy(&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];
+    }
+
+
+    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)
+    {
+        tng_trajectory_destroy(&traj);
+        printf("  Cannot write trajectory headers and box shape.\n");
+        exit(1);
+    }
+
+    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;
+/*
+    This is the main time stepping loop:
+        Compute forces and energies,
+        Update positions, velocities, accelerations.
+*/
+    printf ( "\n" );
+    printf ( "  At each step, 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_print = 0;
+    step_print_index = 0;
+    step_print_num = 10;
+    
+    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;
+
+    /* Saving frequency */
+    step_save = 5;
+
+    wtime = omp_get_wtime ( );
+
+    /* Create a frame set for writing data */
+    tng_num_frames_per_frame_set_get(traj, &n_frames_per_frame_set);
+    if(tng_frame_set_new(traj, 0,
+        n_frames_per_frame_set) != TNG_SUCCESS)
+    {
+        printf("Error creating frame set %d. %s: %d\n",
+                i, __FILE__, __LINE__);
+        exit(1);
+    }
+    frame_set_cnt++;
+
+    /* Setup an empty data array - this will be written to the frame set
+     * (since the block size needs to be correct in the file. After that it
+     * is filled with data. */
+    data = malloc(n_frames_per_frame_set * 3 * np * sizeof(double));
+
+    for(i = 0; i < n_frames_per_frame_set * 3 * np; i++)
+    {
+        data[i] = 0;
+    }
+
+    
+    /* Add data blocks */
+    if(tng_particle_data_block_add(traj, TNG_TRAJ_POSITIONS,
+                                "POSITIONS",
+                                TNG_DOUBLE_DATA,
+                                TNG_TRAJECTORY_BLOCK,
+                                n_frames_per_frame_set, 3,
+                                1, 0, np,
+                                TNG_UNCOMPRESSED,
+                                data) != TNG_SUCCESS)
+    {
+        printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+        free(data);
+        exit(1);
+    }
+    if(tng_particle_data_block_add(traj, TNG_TRAJ_VELOCITIES,
+                                "VELOCITIES",
+                                TNG_DOUBLE_DATA,
+                                TNG_TRAJECTORY_BLOCK,
+                                n_frames_per_frame_set, 3,
+                                1, 0, np,
+                                TNG_UNCOMPRESSED,
+                                data) != TNG_SUCCESS)
+    {
+        printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+        free(data);
+        exit(1);
+    }
+    if(tng_particle_data_block_add(traj, TNG_TRAJ_FORCES,
+                                "FORCES",
+                                TNG_DOUBLE_DATA,
+                                TNG_TRAJECTORY_BLOCK,
+                                n_frames_per_frame_set, 3,
+                                1, 0, np,
+                                TNG_UNCOMPRESSED,
+                                data) != TNG_SUCCESS)
+    {
+        printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+        free(data);
+        exit(1);
+    }
+    
+    /* Write the frame set to disk */
+    if(tng_frame_set_write(traj, TNG_SKIP_HASH) != TNG_SUCCESS)
+    {
+        printf("Error writing frame set. %s: %d\n", __FILE__, __LINE__);
+        free(data);
+        exit(1);
+    }
+
+    if(tng_frame_particle_data_write(traj, 0,
+                                    TNG_TRAJ_POSITIONS, 0, np,
+                                    pos, TNG_USE_HASH) != TNG_SUCCESS)
+    {
+        printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+        free(data);
+        exit(1);
+    }
+    if(tng_frame_particle_data_write(traj, 0,
+                                    TNG_TRAJ_VELOCITIES, 0, np,
+                                    vel, TNG_USE_HASH) != TNG_SUCCESS)
+    {
+        printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+        free(data);
+        exit(1);
+    }
+    if(tng_frame_particle_data_write(traj, 0,
+                                    TNG_TRAJ_FORCES, 0, np,
+                                    force, TNG_USE_HASH) != TNG_SUCCESS)
+    {
+        printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+        free(data);
+        exit(1);
+    }
+    frames_saved_cnt++;
+
+    
+//     printf("pos: %f, vel: %f, force: %f\n", pos[0], vel[0], force[0]);
+        
+    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(frames_saved_cnt % n_frames_per_frame_set == 0)
+            {
+                if(tng_frame_set_new(traj, frames_saved_cnt,
+                    n_frames_per_frame_set) != TNG_SUCCESS)
+                {
+                    printf("Error creating frame set %d. %s: %d\n",
+                            i, __FILE__, __LINE__);
+                    free(data);
+                    exit(1);
+                }
+                /* Add data blocks */
+                if(tng_particle_data_block_add(traj, TNG_TRAJ_POSITIONS,
+                                            "POSITIONS",
+                                            TNG_DOUBLE_DATA,
+                                            TNG_TRAJECTORY_BLOCK,
+                                            n_frames_per_frame_set, 3,
+                                            1, 0, np,
+                                            TNG_UNCOMPRESSED,
+                                            data) != TNG_SUCCESS)
+                {
+                    printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+                    free(data);
+                    exit(1);
+                }
+                if(tng_particle_data_block_add(traj, TNG_TRAJ_VELOCITIES,
+                                            "VELOCITIES",
+                                            TNG_DOUBLE_DATA,
+                                            TNG_TRAJECTORY_BLOCK,
+                                            n_frames_per_frame_set, 3,
+                                            1, 0, np,
+                                            TNG_UNCOMPRESSED,
+                                            data) != TNG_SUCCESS)
+                {
+                    printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+                    free(data);
+                    exit(1);
+                }
+                if(tng_particle_data_block_add(traj, TNG_TRAJ_FORCES,
+                                            "FORCES",
+                                            TNG_DOUBLE_DATA,
+                                            TNG_TRAJECTORY_BLOCK,
+                                            n_frames_per_frame_set, 3,
+                                            1, 0, np,
+                                            TNG_UNCOMPRESSED,
+                                            data) != TNG_SUCCESS)
+                {
+                    printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+                    free(data);
+                    exit(1);
+                }
+
+                /* Write the frame set to disk */
+                if(tng_frame_set_write(traj, TNG_SKIP_HASH) != TNG_SUCCESS)
+                {
+                    printf("Error writing frame set. %s: %d\n", __FILE__, __LINE__);
+                    free(data);
+                    exit(1);
+                }
+
+            }
+            if(tng_frame_particle_data_write(traj, frames_saved_cnt,
+                                            TNG_TRAJ_POSITIONS, 0, np,
+                                            pos, TNG_USE_HASH) != TNG_SUCCESS)
+            {
+                printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+                free(data);
+                exit(1);
+            }
+            if(tng_frame_particle_data_write(traj, frames_saved_cnt,
+                                            TNG_TRAJ_VELOCITIES, 0, np,
+                                            vel, TNG_USE_HASH) != TNG_SUCCESS)
+            {
+                printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+                free(data);
+                exit(1);
+            }
+            if(tng_frame_particle_data_write(traj, frames_saved_cnt,
+                                            TNG_TRAJ_FORCES, 0, np,
+                                            force, TNG_USE_HASH) != TNG_SUCCESS)
+            {
+                printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
+                free(data);
+                exit(1);
+            }
+            frames_saved_cnt++;
+
+        }
+        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.
+*/
+    printf ( "\n" );
+    printf ( "MD_OPENMP\n" );
+    printf ( "  Normal end of execution.\n" );
+
+    printf ( "\n" );
+    timestamp ( );
+
+    return 0;
+}
+/******************************************************************************/
+
+void compute ( int np, int nd, double pos[], double vel[], 
+    double mass, double f[], double *pot, double *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, double POS[ND*NP], the position of each particle.
+
+        Input, double VEL[ND*NP], the velocity of each particle.
+
+        Input, double MASS, the mass of each particle.
+
+        Output, double F[ND*NP], the forces.
+
+        Output, double *POT, the total potential energy.
+
+        Output, double *KIN, the total kinetic energy.
+*/
+{
+    double d;
+    double d2;
+    int i;
+    int j;
+    int k;
+    double ke;
+    double pe;
+    double PI2 = 3.141592653589793 / 2.0;
+    double 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;
+}
+/******************************************************************************/
+
+double dist ( int nd, double r1[], double r2[], double 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, double R1[ND], R2[ND], the positions of the particles.
+
+        Output, double DR[ND], the displacement vector.
+
+        Output, double D, the Euclidean norm of the displacement.
+*/
+{
+    double 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, double box[], int *seed, double pos[], 
+    double vel[], double 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, double BOX[ND], specifies the maximum position
+        of particles in each dimension.
+
+        Input, int *SEED, a seed for the random number generator.
+
+        Output, double POS[ND*NP], the position of each particle.
+
+        Output, double VEL[ND*NP], the velocity of each particle.
+
+        Output, double 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;
+}
+/******************************************************************************/
+
+double 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, double R8_UNIFORM_01, a new pseudorandom variate, strictly between
+        0 and 1.
+*/
+{
+    int k;
+    double r;
+
+    k = *seed / 127773;
+
+    *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
+
+    if ( *seed < 0 )
+    {
+        *seed = *seed + 2147483647;
+    }
+
+    r = ( double ) ( *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;
+    size_t len;
+    time_t now;
+
+    now = time ( NULL );
+    tm = localtime ( &now );
+
+    len = 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, double pos[], double vel[], double f[], 
+    double acc[], double mass, double 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, double POS[ND*NP], the position of each particle.
+
+        Input/output, double VEL[ND*NP], the velocity of each particle.
+
+        Input, double F[ND*NP], the force on each particle.
+
+        Input/output, double ACC[ND*NP], the acceleration of each particle.
+
+        Input, double MASS, the mass of each particle.
+
+        Input, double DT, the time step.
+*/
+{
+    int i;
+    int j;
+    double 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_testing.c b/src/tests/tng_io_testing.c
index 2e7f9ed..c45668f 100644
--- a/src/tests/tng_io_testing.c
+++ b/src/tests/tng_io_testing.c
@@ -1,7 +1,7 @@
 #include <inttypes.h>
 #include <stdlib.h>
 #include <string.h>
-#include "tng_io.h"
+#include <tng_io.h>
 #include "tng_io_testing.h"
 
 static tng_function_status tng_test_setup_molecules(tng_trajectory_t traj)