1 #ifdef TNG_BUILD_OPENMP_EXAMPLES
3 /* This code is part of the tng binary trajectory format.
5 * Written by Magnus Lundborg
6 * Copyright (c) 2012-2013, The GROMACS development team.
7 * Check out http://www.gromacs.org for more information.
10 * This program is free software; you can redistribute it and/or
11 * modify it under the terms of the Revised BSD License.
14 #include "tng/tng_io.h"
20 /* N.B. this code is for testing parallel reading of trajectory frame sets. The
21 * performance is not improved very much and is to a large extent limited by
22 * disk i/o. It can however be used as inspiration for writing parallel code
23 * using the TNG library. The code is NOT fully tested and may behave strangely. */
25 int main(int argc
, char **argv
)
27 tng_trajectory_t traj
, local_traj
= 0;
28 union data_values
***local_positions
= 0; // A 3-dimensional array to be populated
29 union data_values
**particle_pos
= 0;
30 int64_t n_particles
, n_values_per_frame
, n_frame_sets
, n_frames
;
31 int64_t n_frames_per_frame_set
, tot_n_frames
= 0;
34 int64_t particle
= 0, local_first_frame
, local_last_frame
;
35 char atom_name
[64], res_name
[64];
36 tng_trajectory_frame_set_t frame_set
= 0;
40 printf("No file specified\n");
42 printf("tng_parallel_read <tng_file> [particle number = %"PRId64
"]\n",
47 // A reference must be passed to allocate memory
48 if(tng_trajectory_init(&traj
) != TNG_SUCCESS
)
50 tng_trajectory_destroy(&traj
);
53 tng_input_file_set(traj
, argv
[1]);
55 tng_current_frame_set_get(traj
, &frame_set
);
57 // Read the file headers
58 tng_file_headers_read(traj
, TNG_USE_HASH
);
62 particle
= strtoll(argv
[2], 0, 10);
65 tng_num_frame_sets_get(traj
, &n_frame_sets
);
66 tng_num_frames_per_frame_set_get(traj
, &n_frames_per_frame_set
);
68 particle_pos
= malloc(sizeof(union data_values
*) * n_frame_sets
*
69 n_frames_per_frame_set
);
70 for(i
= n_frame_sets
* n_frames_per_frame_set
; i
--;)
72 /* Assume 3 values per frame even if it's not determined yet */
73 particle_pos
[i
] = malloc(sizeof(union data_values
) * 3);
76 printf("%"PRId64
" frame sets\n", n_frame_sets
);
78 if(tng_atom_name_of_particle_nr_get(traj
, particle
, atom_name
,
81 tng_residue_name_of_particle_nr_get(traj
, particle
, res_name
,
85 printf("Particle: %s (%s)\n", atom_name
, res_name
);
89 printf("Particle name not found\n");
94 #pragma omp parallel \
95 private (n_frames, n_particles, n_values_per_frame, \
96 local_first_frame, local_last_frame, j, fail) \
97 firstprivate (local_traj, local_positions, frame_set)\
98 shared(data_type, traj, n_frame_sets, particle_pos, particle, i, tot_n_frames)\
101 /* Each tng_trajectory_t keeps its own file pointers and i/o positions.
102 * Therefore there must be a copy for each thread. */
103 tng_trajectory_init_from_src(traj
, &local_traj
);
105 for(i
= 0; i
< n_frame_sets
; i
++)
107 if(tng_frame_set_nr_find(local_traj
, i
) != TNG_SUCCESS
)
109 printf("FAILED finding frame set %d!\n", i
);
113 if(tng_particle_data_get(local_traj
, TNG_TRAJ_POSITIONS
, &local_positions
,
114 &n_frames
, &n_particles
, &n_values_per_frame
,
115 &data_type
) != TNG_SUCCESS
)
117 printf("FAILED getting particle data\n");
123 tng_current_frame_set_get(local_traj
, &frame_set
);
124 tng_frame_set_frame_range_get(local_traj
, frame_set
, &local_first_frame
, &local_last_frame
);
125 // printf("Frame %"PRId64"-%"PRId64":\n", local_first_frame, local_last_frame);
126 // printf("%"PRId64" %"PRId64" %"PRId64"\n", n_frames, n_particles, n_values_per_frame);
127 tot_n_frames
+= n_frames
;
128 for(j
= 0; j
< n_frames
; j
++)
130 particle_pos
[local_first_frame
+ j
][0] = local_positions
[j
][particle
][0];
131 particle_pos
[local_first_frame
+ j
][1] = local_positions
[j
][particle
][1];
132 particle_pos
[local_first_frame
+ j
][2] = local_positions
[j
][particle
][2];
140 tng_particle_data_values_free(local_traj
, local_positions
, n_frames
, n_particles
,
141 n_values_per_frame
, data_type
);
143 tng_trajectory_destroy(&local_traj
);
148 for(j
= 0; j
< tot_n_frames
; j
++)
150 printf("\t%"PRId64
"\t%"PRId64
"\t%"PRId64
"\n", particle_pos
[j
][0].i
,
151 particle_pos
[j
][1].i
, particle_pos
[j
][2].i
);
155 for(j
= 0; j
< tot_n_frames
; j
++)
157 printf("\t%f\t%f\t%f\n", particle_pos
[j
][0].f
,
158 particle_pos
[j
][1].f
, particle_pos
[j
][2].f
);
161 case TNG_DOUBLE_DATA
:
162 for(j
= 0; j
< tot_n_frames
; j
++)
164 printf("\t%f\t%f\t%f\n", particle_pos
[j
][0].d
,
165 particle_pos
[j
][1].d
, particle_pos
[j
][2].d
);
172 /* Free more memory */
173 for(i
= n_frame_sets
* n_frames_per_frame_set
; i
--;)
175 free(particle_pos
[i
]);
179 tng_trajectory_destroy(&traj
);