1 #ifdef TNG_BUILD_OPENMP_EXAMPLES
3 /* This code is part of the tng binary trajectory format.
7 * Written by Magnus Lundborg
8 * Copyright (c) 2012-2013, The GROMACS development team.
9 * Check out http://www.gromacs.org for more information.
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the Revised BSD License.
18 #include "../../include/tng_io.h"
21 /* N.B. this code is for testing parallel reading of trajectory frame sets. The
22 * performance is not improved very much and is to a large extent limited by
23 * disk i/o. It can however be used as inspiration for writing parallel code
24 * using the TNG library. The code is NOT fully tested and may behave strangely. */
26 int main(int argc, char **argv)
28 tng_trajectory_t traj, local_traj = 0;
29 union data_values ***local_positions = 0; // A 3-dimensional array to be populated
30 union data_values **particle_pos = 0;
31 int64_t n_particles, n_values_per_frame, n_frame_sets, n_frames;
32 int64_t n_frames_per_frame_set, tot_n_frames = 0;
35 int64_t particle = 0, local_first_frame, local_last_frame;
36 char atom_name[64], res_name[64];
37 tng_trajectory_frame_set_t frame_set = 0;
41 printf("No file specified\n");
43 printf("tng_parallel_read <tng_file> [particle number = %"PRId64"]\n",
48 // A reference must be passed to allocate memory
49 if(tng_trajectory_init(&traj) != TNG_SUCCESS)
51 tng_trajectory_destroy(&traj);
54 tng_input_file_set(traj, argv[1]);
56 tng_current_frame_set_get(traj, &frame_set);
58 // Read the file headers
59 tng_file_headers_read(traj, TNG_USE_HASH);
63 particle = strtoll(argv[2], 0, 10);
66 tng_num_frame_sets_get(traj, &n_frame_sets);
67 tng_num_frames_per_frame_set_get(traj, &n_frames_per_frame_set);
69 particle_pos = malloc(sizeof(union data_values *) * n_frame_sets *
70 n_frames_per_frame_set);
71 for(i = n_frame_sets * n_frames_per_frame_set; i--;)
73 /* Assume 3 values per frame even if it's not determined yet */
74 particle_pos[i] = malloc(sizeof(union data_values) * 3);
77 printf("%"PRId64" frame sets\n", n_frame_sets);
79 if(tng_atom_name_of_particle_nr_get(traj, particle, atom_name,
82 tng_residue_name_of_particle_nr_get(traj, particle, res_name,
86 printf("Particle: %s (%s)\n", atom_name, res_name);
90 printf("Particle name not found\n");
95 #pragma omp parallel \
96 private (n_frames, n_particles, n_values_per_frame, \
97 local_first_frame, local_last_frame, j, fail) \
98 firstprivate (local_traj, local_positions, frame_set)\
99 shared(data_type, traj, n_frame_sets, particle_pos, particle, i, tot_n_frames)\
102 /* Each tng_trajectory_t keeps its own file pointers and i/o positions.
103 * Therefore there must be a copy for each thread. */
104 tng_trajectory_init_from_src(traj, &local_traj);
106 for(i = 0; i < n_frame_sets; i++)
108 if(tng_frame_set_nr_find(local_traj, i) != TNG_SUCCESS)
110 printf("FAILED finding frame set %d!\n", i);
114 if(tng_particle_data_get(local_traj, TNG_TRAJ_POSITIONS, &local_positions,
115 &n_frames, &n_particles, &n_values_per_frame,
116 &data_type) != TNG_SUCCESS)
118 printf("FAILED getting particle data\n");
124 tng_current_frame_set_get(local_traj, &frame_set);
125 tng_frame_set_frame_range_get(local_traj, frame_set, &local_first_frame, &local_last_frame);
126 // printf("Frame %"PRId64"-%"PRId64":\n", local_first_frame, local_last_frame);
127 // printf("%"PRId64" %"PRId64" %"PRId64"\n", n_frames, n_particles, n_values_per_frame);
128 tot_n_frames += n_frames;
129 for(j = 0; j < n_frames; j++)
131 particle_pos[local_first_frame + j][0] = local_positions[j][particle][0];
132 particle_pos[local_first_frame + j][1] = local_positions[j][particle][1];
133 particle_pos[local_first_frame + j][2] = local_positions[j][particle][2];
141 tng_particle_data_values_free(local_traj, local_positions, n_frames, n_particles,
142 n_values_per_frame, data_type);
144 tng_trajectory_destroy(&local_traj);
149 for(j = 0; j < tot_n_frames; j++)
151 printf("\t%"PRId64"\t%"PRId64"\t%"PRId64"\n", particle_pos[j][0].i,
152 particle_pos[j][1].i, particle_pos[j][2].i);
156 for(j = 0; j < tot_n_frames; j++)
158 printf("\t%f\t%f\t%f\n", particle_pos[j][0].f,
159 particle_pos[j][1].f, particle_pos[j][2].f);
162 case TNG_DOUBLE_DATA:
163 for(j = 0; j < tot_n_frames; j++)
165 printf("\t%f\t%f\t%f\n", particle_pos[j][0].d,
166 particle_pos[j][1].d, particle_pos[j][2].d);
173 /* Free more memory */
174 for(i = n_frame_sets * n_frames_per_frame_set; i--;)
176 free(particle_pos[i]);
180 tng_trajectory_destroy(&traj);