Bumped TNG to latest version.
[gromacs.git] / src / external / tng_io / src / tests / tng_io_testing.c
blob973ac4fb7b6268c0c5f2e9133f2133382fb9e7db
1 /* This code is part of the tng binary trajectory format.
3 * VERSION 1.0
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 #ifdef USE_STD_INTTYPES_H
15 #include <inttypes.h>
16 #endif
18 #include <stdlib.h>
19 #include <string.h>
20 #include "../../include/tng_io.h"
22 static tng_function_status tng_test_setup_molecules(tng_trajectory_t traj)
24 tng_molecule_t molecule;
25 tng_chain_t chain;
26 tng_residue_t residue;
27 tng_atom_t atom;
28 int64_t cnt;
30 tng_molecule_add(traj, "water", &molecule);
31 tng_molecule_chain_add(traj, molecule, "W", &chain);
32 tng_chain_residue_add(traj, chain, "WAT", &residue);
33 if(tng_residue_atom_add(traj, residue, "O", "O", &atom) == TNG_CRITICAL)
35 return(TNG_CRITICAL);
37 if(tng_residue_atom_add(traj, residue, "HO1", "H", &atom) == TNG_CRITICAL)
39 return(TNG_CRITICAL);
41 if(tng_residue_atom_add(traj, residue, "HO2", "H", &atom) == TNG_CRITICAL)
43 return(TNG_CRITICAL);
45 tng_molecule_cnt_set(traj, molecule, 200);
46 tng_molecule_cnt_get(traj, molecule, &cnt);
47 /* printf("Created %"PRId64" %s molecules.\n", cnt, molecule->name); */
49 /* traj->molecule_cnt_list[traj->n_molecules-1] = 5;
50 // tng_molecule_name_set(traj, &traj->molecules[1], "ligand");
51 // tng_molecule_name_set(traj, &traj->molecules[2], "water");
52 // tng_molecule_name_set(traj, &traj->molecules[3], "dummy");
53 // traj->molecules[0].id = 0;
54 // traj->molecules[1].id = 1;
55 // traj->molecules[2].id = 2;
56 // traj->molecules[3].id = 3;
58 // if(tng_add_atom_to_molecule(traj, &traj->molecules[0], "atom1", "type1") == TNG_CRITICAL)
59 // {
60 // return(TNG_CRITICAL);
61 // }
62 // if(tng_add_atom_to_molecule(traj, &traj->molecules[0], "atom2", "type1") == TNG_CRITICAL)
63 // {
64 // return(TNG_CRITICAL);
65 // }
66 // if(tng_add_atom_to_molecule(traj, &traj->molecules[0], "atom3", "type1") == TNG_CRITICAL)
67 // {
68 // return(TNG_CRITICAL);
69 // }
70 // if(tng_add_atom_to_molecule(traj, &traj->molecules[0], "atom4", "type2") == TNG_CRITICAL)
71 // {
72 // return(TNG_CRITICAL);
73 // }
74 // if(tng_add_atom_to_molecule(traj, &traj->molecules[0], "atom5", "type2") == TNG_CRITICAL)
75 // {
76 // return(TNG_CRITICAL);
77 // }
78 // if(tng_add_atom_to_molecule(traj, &traj->molecules[0], "atom6", "type2") == TNG_CRITICAL)
79 // {
80 // return(TNG_CRITICAL);
81 // }
82 // if(tng_add_atom_to_molecule(traj, &traj->molecules[0], "atom7", "type3") == TNG_CRITICAL)
83 // {
84 // return(TNG_CRITICAL);
85 // }
86 // if(tng_add_atom_to_molecule(traj, &traj->molecules[1], "C1", "C") == TNG_CRITICAL)
87 // {
88 // return(TNG_CRITICAL);
89 // }
90 // if(tng_add_atom_to_molecule(traj, &traj->molecules[1], "O1", "O") == TNG_CRITICAL)
91 // {
92 // return(TNG_CRITICAL);
93 // }
94 // if(tng_add_atom_to_molecule(traj, &traj->molecules[1], "H11", "H") == TNG_CRITICAL)
95 // {
96 // return(TNG_CRITICAL);
97 // }
98 // if(tng_add_atom_to_molecule(traj, &traj->molecules[1], "H12", "H") == TNG_CRITICAL)
99 // {
100 // return(TNG_CRITICAL);
101 // }
102 // if(tng_add_atom_to_molecule(traj, &traj->molecules[1], "H13", "H") == TNG_CRITICAL)
103 // {
104 // return(TNG_CRITICAL);
105 // }
107 return(TNG_SUCCESS);
110 static tng_function_status tng_test_read_and_write_file
111 (tng_trajectory_t traj)
113 tng_function_status stat;
115 stat = tng_file_headers_read(traj, TNG_USE_HASH);
116 if(stat == TNG_CRITICAL)
118 return(stat);
120 stat = tng_file_headers_write(traj, TNG_USE_HASH);
121 if(stat == TNG_CRITICAL)
123 return(stat);
126 while(stat == TNG_SUCCESS)
128 stat = tng_frame_set_read_next(traj, TNG_USE_HASH);
129 if(stat != TNG_SUCCESS)
131 return(stat);
133 stat = tng_frame_set_write(traj, TNG_USE_HASH);
136 return(stat);
139 static tng_function_status tng_test_write_and_read_traj(tng_trajectory_t *traj)
141 int i, j, k, nr, cnt;
142 float *data, *molpos, *charges;
143 int64_t mapping[300], n_particles, n_frames_per_frame_set, tot_n_mols;
144 // int64_t frame_nr;
145 double box_shape[9];
146 char atom_type[16], annotation[128];
147 tng_function_status stat = TNG_SUCCESS;
149 tng_medium_stride_length_set(*traj, 10);
150 tng_long_stride_length_set(*traj, 100);
152 /* Create molecules */
153 if(tng_test_setup_molecules(*traj) == TNG_CRITICAL)
155 return(TNG_CRITICAL);
158 /* Set the box shape */
159 box_shape[1] = box_shape[2] = box_shape[3] = box_shape[5] = box_shape[6] =
160 box_shape[7] = 0;
161 box_shape[0] = 150.0;
162 box_shape[4] = 145.5;
163 box_shape[8] = 155.5;
164 if(tng_data_block_add(*traj, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE", TNG_DOUBLE_DATA,
165 TNG_NON_TRAJECTORY_BLOCK, 1, 9, 1, TNG_UNCOMPRESSED,
166 box_shape) == TNG_CRITICAL)
168 tng_trajectory_destroy(traj);
169 printf("Cannot write trajectory box shape.\n");
170 exit(1);
173 /* Set partial charges (treat the water as TIP3P. */
174 tng_num_particles_get(*traj, &n_particles);
175 charges = malloc(sizeof(float) * n_particles);
176 for(i = 0; i < n_particles; i++)
178 stat = tng_atom_type_of_particle_nr_get(*traj, i, atom_type,
179 sizeof(atom_type));
180 if(stat == TNG_CRITICAL)
182 break;
184 if(atom_type[0] == 'O')
186 charges[i] = -0.834;
188 else if(atom_type[0] == 'H')
190 charges[i] = 0.417;
193 if(stat == TNG_CRITICAL)
195 free(charges);
196 printf("Failed setting partial charges.\n");
197 return(TNG_CRITICAL);
200 stat = tng_particle_data_block_add(*traj, TNG_TRAJ_PARTIAL_CHARGES, "PARTIAL CHARGES",
201 TNG_FLOAT_DATA, TNG_NON_TRAJECTORY_BLOCK,
202 1, 1, 1, 0, n_particles,
203 TNG_UNCOMPRESSED, charges);
204 free(charges);
205 if(stat != TNG_SUCCESS)
207 printf("Failed adding partial charges\n");
208 return(TNG_CRITICAL);
212 /* Generate a custom annotation data block */
213 strcpy(annotation, "This trajectory was generated from tng_io_testing. "
214 "It is not a real MD trajectory.");
215 if(tng_data_block_add(*traj, 10100, "DETAILS", TNG_CHAR_DATA,
216 TNG_NON_TRAJECTORY_BLOCK, 1, 1, 1, TNG_UNCOMPRESSED,
217 annotation) != TNG_SUCCESS)
219 printf("Failed adding details annotation data block.\n");
220 return(TNG_CRITICAL);
223 /* Write file headers (includes non trajectory data blocks */
224 if(tng_file_headers_write(*traj, TNG_SKIP_HASH) == TNG_CRITICAL)
226 printf("Cannot write file headers.\n");
230 tng_num_frames_per_frame_set_get(*traj, &n_frames_per_frame_set);
232 data = malloc(sizeof(float) * n_particles *
233 n_frames_per_frame_set * 3);
234 if(!data)
236 printf("Cannot allocate memory. %s: %d\n", __FILE__, __LINE__);
237 return(TNG_CRITICAL);
240 tng_num_molecules_get(*traj, &tot_n_mols);
241 molpos = malloc(sizeof(float) * tot_n_mols * 3);
243 /* Set initial coordinates */
244 for(i = 0; i < tot_n_mols; i++)
246 nr = i * 3;
247 /* Somewhat random coordinates (between 0 and 100),
248 * but not specifying a random seed */
249 molpos[nr] = 100.0 * rand() / (RAND_MAX + 1.0);
250 molpos[nr+1] = 100.0 * rand() / (RAND_MAX + 1.0);
251 molpos[nr+2] = 100.0 * rand() / (RAND_MAX + 1.0);
254 /* Generate 200 frame sets - each with 100 frames (by default) */
255 for(i = 0; i < 200; i++)
257 cnt = 0;
258 for(j = 0; j < n_frames_per_frame_set; j++)
260 for(k = 0; k < tot_n_mols; k++)
262 nr = k * 3;
263 /* Move -1 to 1 */
264 molpos[nr] += 2 * (rand() / (RAND_MAX + 1.0)) - 1;
265 molpos[nr+1] += 2 * (rand() / (RAND_MAX + 1.0)) - 1;
266 molpos[nr+2] += 2 * (rand() / (RAND_MAX + 1.0)) - 1;
268 data[cnt++] = molpos[nr];
269 data[cnt++] = molpos[nr + 1];
270 data[cnt++] = molpos[nr + 2];
271 data[cnt++] = molpos[nr] + 1;
272 data[cnt++] = molpos[nr + 1] + 1;
273 data[cnt++] = molpos[nr + 2] + 1;
274 data[cnt++] = molpos[nr] - 1;
275 data[cnt++] = molpos[nr + 1] - 1;
276 data[cnt++] = molpos[nr + 2] - 1;
279 if(tng_frame_set_new(*traj, i * n_frames_per_frame_set,
280 n_frames_per_frame_set) != TNG_SUCCESS)
282 printf("Error creating frame set %d. %s: %d\n",
283 i, __FILE__, __LINE__);
284 free(molpos);
285 free(data);
286 return(TNG_CRITICAL);
289 tng_frame_set_particle_mapping_free(*traj);
291 /* Setup particle mapping. Use 4 different mapping blocks with arbitrary
292 * mappings. */
293 for(k=0; k<150; k++)
295 mapping[k]=k;
297 if(tng_particle_mapping_add(*traj, 0, 150, mapping) != TNG_SUCCESS)
299 printf("Error creating particle mapping. %s: %d\n",
300 __FILE__, __LINE__);
301 free(molpos);
302 free(data);
303 return(TNG_CRITICAL);
305 for(k=0; k<150; k++)
307 mapping[k]=599-k;
309 if(tng_particle_mapping_add(*traj, 150, 150, mapping) != TNG_SUCCESS)
311 printf("Error creating particle mapping. %s: %d\n",
312 __FILE__, __LINE__);
313 free(molpos);
314 free(data);
315 return(TNG_CRITICAL);
317 for(k=0; k<150; k++)
319 mapping[k]=k+150;
321 if(tng_particle_mapping_add(*traj, 300, 150, mapping) != TNG_SUCCESS)
323 printf("Error creating particle mapping. %s: %d\n",
324 __FILE__, __LINE__);
325 free(molpos);
326 free(data);
327 return(TNG_CRITICAL);
329 for(k=0; k<150; k++)
331 mapping[k]=449-k;
333 if(tng_particle_mapping_add(*traj, 450, 150, mapping) != TNG_SUCCESS)
335 printf("Error creating particle mapping. %s: %d\n",
336 __FILE__, __LINE__);
337 free(molpos);
338 free(data);
339 return(TNG_CRITICAL);
342 /* Add the positions in a data block */
343 if(tng_particle_data_block_add(*traj, TNG_TRAJ_POSITIONS,
344 "POSITIONS",
345 TNG_FLOAT_DATA,
346 TNG_TRAJECTORY_BLOCK,
347 n_frames_per_frame_set, 3,
348 1, 0, n_particles,
349 /* TNG_UNCOMPRESSED, */
350 TNG_GZIP_COMPRESSION,
351 data) != TNG_SUCCESS)
353 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
354 free(molpos);
355 free(data);
356 return(TNG_CRITICAL);
358 /* Write the frame set */
359 if(tng_frame_set_write(*traj, TNG_SKIP_HASH) != TNG_SUCCESS)
361 printf("Error writing frame set. %s: %d\n", __FILE__, __LINE__);
362 free(molpos);
363 free(data);
364 return(TNG_CRITICAL);
368 /* Write two more frame sets one frame at a time */
370 /* Make a new frame set - if always using the same mapping blocks
371 * it is not necessary to explicitly add a new frame set - it will
372 * be added automatically when adding data for a frame */
373 /* if(tng_frame_set_new(*traj, i * n_frames_per_frame_set,
374 n_frames_per_frame_set) != TNG_SUCCESS)
376 printf("Error creating frame set %d. %s: %d\n",
377 i, __FILE__, __LINE__);
378 free(molpos);
379 free(data);
380 return(TNG_CRITICAL);
383 frame_nr = i * n_frames_per_frame_set;
385 for(k=0; k<300; k++)
387 mapping[k]=k;
389 *//* Just use two particle mapping blocks in this frame set *//*
390 if(tng_particle_mapping_add(*traj, 0, 300, mapping) != TNG_SUCCESS)
392 printf("Error creating particle mapping. %s: %d\n",
393 __FILE__, __LINE__);
394 free(molpos);
395 free(data);
396 return(TNG_CRITICAL);
398 for(k=0; k<300; k++)
400 mapping[k]=599-k;
402 if(tng_particle_mapping_add(*traj, 300, 300, mapping) != TNG_SUCCESS)
404 printf("Error creating particle mapping. %s: %d\n",
405 __FILE__, __LINE__);
406 free(molpos);
407 free(data);
408 return(TNG_CRITICAL);
411 *//* Add the data block to the current frame set *//*
412 if(tng_particle_data_block_add(*traj, TNG_TRAJ_POSITIONS,
413 "POSITIONS",
414 TNG_FLOAT_DATA,
415 TNG_TRAJECTORY_BLOCK,
416 n_frames_per_frame_set, 3,
417 1, 0, n_particles,
418 TNG_UNCOMPRESSED,
419 0) != TNG_SUCCESS)
421 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
422 free(molpos);
423 free(data);
424 return(TNG_CRITICAL);
427 *//* Write the frame set to disk *//*
428 if(tng_frame_set_write(*traj, TNG_SKIP_HASH) != TNG_SUCCESS)
430 printf("Error writing frame set. %s: %d\n", __FILE__, __LINE__);
431 free(molpos);
432 free(data);
433 return(TNG_CRITICAL);
436 *//* Write particle data to disk - one frame at a time *//*
437 for(i = 0; i < n_frames_per_frame_set * 2; i++)
439 for(j = 0; j < 2; j++)
441 cnt = 0;
442 for(k = 0; k < tot_n_mols/2; k++)
444 nr = k * 3;
445 *//* Move -1 to 1 *//*
446 molpos[nr] += 2 * (rand() / (RAND_MAX + 1.0)) - 1;
447 molpos[nr+1] += 2 * (rand() / (RAND_MAX + 1.0)) - 1;
448 molpos[nr+2] += 2 * (rand() / (RAND_MAX + 1.0)) - 1;
450 data[cnt++] = molpos[nr];
451 data[cnt++] = molpos[nr + 1];
452 data[cnt++] = molpos[nr + 2];
453 data[cnt++] = molpos[nr] + 1;
454 data[cnt++] = molpos[nr + 1] + 1;
455 data[cnt++] = molpos[nr + 2] + 1;
456 data[cnt++] = molpos[nr] - 1;
457 data[cnt++] = molpos[nr + 1] - 1;
458 data[cnt++] = molpos[nr + 2] - 1;
460 if(tng_frame_particle_data_write(*traj, frame_nr + i,
461 TNG_TRAJ_POSITIONS, j * 300, 300,
462 data, TNG_SKIP_HASH) != TNG_SUCCESS)
464 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
465 free(molpos);
466 free(data);
467 return(TNG_CRITICAL);
472 free(molpos);
473 free(data);
475 tng_trajectory_destroy(traj);
476 tng_trajectory_init(traj);
477 tng_input_file_set(*traj, TNG_EXAMPLE_FILES_DIR "tng_test.tng");
479 stat = tng_file_headers_read(*traj, TNG_SKIP_HASH);
481 while(stat == TNG_SUCCESS)
483 stat = tng_frame_set_read_next(*traj, TNG_SKIP_HASH);
486 return(stat);
489 /* This test relies on knowing that the box shape is stored as double */
490 tng_function_status tng_test_get_box_data(tng_trajectory_t traj)
492 int64_t n_frames, n_values_per_frame;
493 union data_values **values = 0;
494 char type;
496 if(tng_data_get(traj, TNG_TRAJ_BOX_SHAPE, &values, &n_frames,
497 &n_values_per_frame, &type) != TNG_SUCCESS)
499 printf("Failed getting box shape. %s: %d\n", __FILE__, __LINE__);
500 return(TNG_CRITICAL);
504 // int64_t i, j;
505 // printf("Box shape:");
506 // for(i=0; i<n_frames; i++)
507 // {
508 // for(j=0; j<n_values_per_frame; j++)
509 // {
510 // printf("\t%f", (values[i][j]).d);
511 // }
512 // printf("\n");
513 // }
515 tng_data_values_free(traj, values, n_frames, n_values_per_frame, type);
517 return(TNG_SUCCESS);
520 /* This test relies on knowing that the positions are stored as float
521 * and that the data is not sparse (i.e. as many frames in the data
522 * as in the frame set */
523 tng_function_status tng_test_get_positions_data(tng_trajectory_t traj)
525 int64_t n_frames, n_particles, n_values_per_frame;
526 union data_values ***values = 0;
527 char type;
529 if(tng_particle_data_get(traj, TNG_TRAJ_POSITIONS, &values, &n_frames,
530 &n_particles, &n_values_per_frame, &type) !=
531 TNG_SUCCESS)
533 printf("Failed getting particle positions. %s: %d\n", __FILE__, __LINE__);
534 return(TNG_CRITICAL);
538 // int64_t i, j, k;
539 // struct tng_trajectory_frame_set *frame_set =
540 // &traj->current_trajectory_frame_set;
541 // for(i = 0; i<n_frames; i++)
542 // {
543 // printf("Frame %"PRId64"\n", frame_set->first_frame + i);
544 // for(j = 0; j<n_particles; j++)
545 // {
546 // printf("Particle %"PRId64":", j);
547 // for(k=0; k<n_values_per_frame; k++)
548 // {
549 // printf("\t%f", (values[i][j][k]).f);
550 // }
551 // printf("\n");
552 // }
553 // }
555 tng_particle_data_values_free(traj, values, n_frames, n_particles,
556 n_values_per_frame, type);
558 values = 0;
560 tng_particle_data_interval_get(traj, TNG_TRAJ_POSITIONS, 11000, 11499,
561 TNG_SKIP_HASH, &values, &n_particles,
562 &n_values_per_frame, &type);
564 /* Here the particle positions can be printed */
566 tng_particle_data_values_free(traj, values, 500, n_particles,
567 n_values_per_frame, type);
569 return(TNG_SUCCESS);
572 int main()
574 tng_trajectory_t traj;
575 tng_function_status stat;
576 char time_str[TNG_MAX_DATE_STR_LEN];
578 if(tng_trajectory_init(&traj) != TNG_SUCCESS)
580 tng_trajectory_destroy(&traj);
581 printf("Test Init trajectory:\t\t\t\tFailed. %s: %d.\n",
582 __FILE__, __LINE__);
583 exit(1);
585 printf("Test Init trajectory:\t\t\t\tSucceeded.\n");
587 tng_time_get_str(traj, time_str);
589 printf("Creation time: %s\n", time_str);
591 tng_input_file_set(traj, TNG_EXAMPLE_FILES_DIR "tng_example.tng");
592 tng_output_file_set(traj, TNG_EXAMPLE_FILES_DIR "tng_example_out.tng");
595 if(tng_test_read_and_write_file(traj) == TNG_CRITICAL)
597 printf("Test Read and write file:\t\t\tFailed. %s: %d\n",
598 __FILE__, __LINE__);
600 else
602 printf("Test Read and write file:\t\t\tSucceeded.\n");
605 if(tng_test_get_box_data(traj) != TNG_SUCCESS)
607 printf("Test Get data:\t\t\t\t\tFailed. %s: %d\n",
608 __FILE__, __LINE__);
610 else
612 printf("Test Get data:\t\t\t\t\tSucceeded.\n");
615 if(tng_trajectory_destroy(&traj) == TNG_CRITICAL ||
616 tng_trajectory_init(&traj) == TNG_CRITICAL)
618 printf("Test Destroy and init trajectory:\t\tFailed. %s: %d\n",
619 __FILE__, __LINE__);
621 else
623 printf("Test Destroy and init trajectory:\t\tSucceeded.\n");
627 tng_output_file_set(traj, TNG_EXAMPLE_FILES_DIR "tng_test.tng");
629 if(tng_test_write_and_read_traj(&traj) == TNG_CRITICAL)
631 printf("Test Write and read file:\t\t\tFailed. %s: %d\n",
632 __FILE__, __LINE__);
634 else
636 printf("Test Write and read file:\t\t\tSucceeded.\n");
639 if(tng_test_get_positions_data(traj) != TNG_SUCCESS)
641 printf("Test Get particle data:\t\t\t\tFailed. %s: %d\n",
642 __FILE__, __LINE__);
644 else
646 printf("Test Get particle data:\t\t\t\tSucceeded.\n");
649 if(tng_trajectory_destroy(&traj) == TNG_CRITICAL)
651 printf("Test Destroy trajectory:\t\t\tFailed. %s: %d.\n",
652 __FILE__, __LINE__);
653 exit(1);
655 else
657 printf("Test Destroy trajectory:\t\t\tSucceeded.\n");
661 stat = tng_util_trajectory_open(TNG_EXAMPLE_FILES_DIR "tng_test.tng", 'r', &traj);
663 if(stat != TNG_SUCCESS)
665 printf("Test Utility function open:\t\t\tFailed. %s: %d.\n",
666 __FILE__, __LINE__);
667 exit(1);
669 else
671 printf("Test Utility function open:\t\t\tSucceeded.\n");
674 stat = tng_util_trajectory_close(&traj);
675 if(stat != TNG_SUCCESS)
677 printf("Test Utility function close:\t\t\tFailed. %s: %d.\n",
678 __FILE__, __LINE__);
679 exit(1);
681 else
683 printf("Test Utility function close:\t\t\tSucceeded.\n");
686 printf("Tests finished\n");
688 exit(0);