diff options
| -rw-r--r-- | include/tng_io.h | 26 | ||||
| -rw-r--r-- | src/lib/tng_io.c | 165 | ||||
| -rw-r--r-- | src/tests/CMakeLists.txt | 14 | ||||
| -rw-r--r-- | src/tests/md_openmp.c | 4 | ||||
| -rw-r--r-- | src/tests/md_openmp_util.c | 707 | ||||
| -rw-r--r-- | src/tests/tng_io_read_pos_util.c | 113 |
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, ¤t_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, ¤t_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, ¤t_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, ¤t_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); +} |