1 #ifdef TNG_BUILD_OPENMP_EXAMPLES
2
3 #include "tng/tng_io.h"
4 #include <stdlib.h>
5 #include <stdio.h>
6 #include <time.h>
7 #include <math.h>
8 #include <omp.h>
9
10 int main ();
11 void compute ( int np, int nd, float pos[], float vel[],
12 float mass, float f[], float *pot, float *kin );
13 float dist ( int nd, float r1[], float r2[], float dr[] );
14 void initialize ( int np, int nd, float box[], int *seed, float pos[],
15 float vel[], float acc[] );
16 float r8_uniform_01 ( int *seed );
17 void timestamp ( void );
18 void update ( int np, int nd, float pos[], float vel[], float f[],
19 float acc[], float mass, float dt );
20
21 /******************************************************************************/
22
main()23 int main ()
24
25 /******************************************************************************/
26 /*
27 Purpose:
28
29 MAIN is the main program for MD_OPENMP.
30
31 Discussion:
32
33 MD implements a simple molecular dynamics simulation.
34
35 The program uses Open MP directives to allow parallel computation.
36
37 The velocity Verlet time integration scheme is used.
38
39 The particles interact with a central pair potential.
40
41 Output of the program is saved in the TNG format, which is why this
42 code is included in the TNG API release. The high-level API of the
43 TNG API is used where appropriate.
44
45 Licensing:
46
47 This code is distributed under the GNU LGPL license.
48
49 Modified:
50
51 8 Jan 2013
52
53 Author:
54
55 Original FORTRAN77 version by Bill Magro.
56 C version by John Burkardt.
57 TNG trajectory output by Magnus Lundborg.
58
59 Parameters:
60
61 None
62 */
63 {
64 float *acc;
65 float *box;
66 float *box_shape;
67 float dt = 0.0002;
68 float e0;
69 float *force;
70 int i;
71 float kinetic;
72 float mass = 1.0;
73 int nd = 3;
74 int np = 50;
75 float *pos;
76 float potential;
77 int proc_num;
78 int seed = 123456789;
79 int step;
80 int step_num = 50000;
81 int step_print;
82 int step_print_index;
83 int step_print_num;
84 int step_save;
85 float *vel;
86 float wtime;
87 tng_trajectory_t traj;
88 tng_molecule_t molecule;
89 tng_chain_t chain;
90 tng_residue_t residue;
91 tng_atom_t atom;
92
93 timestamp ( );
94
95 proc_num = omp_get_num_procs ( );
96
97 acc = ( float * ) malloc ( nd * np * sizeof ( float ) );
98 box = ( float * ) malloc ( nd * sizeof ( float ) );
99 box_shape = (float *) malloc (9 * sizeof (float));
100 force = ( float * ) malloc ( nd * np * sizeof ( float ) );
101 pos = ( float * ) malloc ( nd * np * sizeof ( float ) );
102 vel = ( float * ) malloc ( nd * np * sizeof ( float ) );
103
104 printf ( "\n" );
105 printf ( "MD_OPENMP\n" );
106 printf ( " C/OpenMP version\n" );
107 printf ( "\n" );
108 printf ( " A molecular dynamics program.\n" );
109
110 printf ( "\n" );
111 printf ( " NP, the number of particles in the simulation is %d\n", np );
112 printf ( " STEP_NUM, the number of time steps, is %d\n", step_num );
113 printf ( " DT, the size of each time step, is %f\n", dt );
114
115 printf ( "\n" );
116 printf ( " Number of processors available = %d\n", proc_num );
117 printf ( " Number of threads = %d\n", omp_get_max_threads ( ) );
118
119
120 printf("\n");
121 printf(" Initializing trajectory storage.\n");
122 /* Initialize the TNG trajectory */
123 tng_util_trajectory_open(TNG_EXAMPLE_FILES_DIR "tng_md_out.tng", 'w', &traj);
124
125
126
127 /* Set molecules data */
128 /* N.B. This is still not done using utility functions. The low-level API
129 * is used. */
130 printf(" Creating molecules in trajectory.\n");
131 tng_molecule_add(traj, "water", &molecule);
132 tng_molecule_chain_add(traj, molecule, "W", &chain);
133 tng_chain_residue_add(traj, chain, "WAT", &residue);
134 if(tng_residue_atom_add(traj, residue, "O", "O", &atom) == TNG_CRITICAL)
135 {
136 tng_util_trajectory_close(&traj);
137 printf(" Cannot create molecules.\n");
138 exit(1);
139 }
140 tng_molecule_cnt_set(traj, molecule, np);
141
142
143 /*
144 Set the dimensions of the box.
145 */
146 for(i = 0; i < 9; i++)
147 {
148 box_shape[i] = 0.0;
149 }
150 for ( i = 0; i < nd; i++ )
151 {
152 box[i] = 10.0;
153 /* box_shape stores 9 values according to the TNG specs */
154 box_shape[i*nd + i] = box[i];
155 }
156
157
158 printf ( "\n" );
159 printf ( " Initializing positions, velocities, and accelerations.\n" );
160 /*
161 Set initial positions, velocities, and accelerations.
162 */
163 initialize ( np, nd, box, &seed, pos, vel, acc );
164 /*
165 Compute the forces and energies.
166 */
167 printf ( "\n" );
168 printf ( " Computing initial forces and energies.\n" );
169
170 compute ( np, nd, pos, vel, mass, force, &potential, &kinetic );
171
172 e0 = potential + kinetic;
173
174 /* Saving frequency */
175 step_save = 400;
176
177 step_print = 0;
178 step_print_index = 0;
179 step_print_num = 10;
180
181 /*
182 This is the main time stepping loop:
183 Compute forces and energies,
184 Update positions, velocities, accelerations.
185 */
186 printf(" Every %d steps box shape, particle positions, velocities and forces are\n",
187 step_save);
188 printf(" saved to a TNG trajectory file.\n");
189 printf ( "\n" );
190 printf ( " At certain step intervals, we report the potential and kinetic energies.\n" );
191 printf ( " The sum of these energies should be a constant.\n" );
192 printf ( " As an accuracy check, we also print the relative error\n" );
193 printf ( " in the total energy.\n" );
194 printf ( "\n" );
195 printf ( " Step Potential Kinetic (P+K-E0)/E0\n" );
196 printf ( " Energy P Energy K Relative Energy Error\n" );
197 printf ( "\n" );
198
199 step = 0;
200 printf ( " %8d %14f %14f %14e\n",
201 step, potential, kinetic, ( potential + kinetic - e0 ) / e0 );
202 step_print_index++;
203 step_print = ( step_print_index * step_num ) / step_print_num;
204
205 /* Set the output frequency of box shape, positions, velocities and forces */
206 if(tng_util_box_shape_write_frequency_set(traj, step_save) != TNG_SUCCESS)
207 {
208 printf("Error setting writing frequency data. %s: %d\n",
209 __FILE__, __LINE__);
210 exit(1);
211 }
212 if(tng_util_pos_write_frequency_set(traj, step_save) != TNG_SUCCESS)
213 {
214 printf("Error setting writing frequency data. %s: %d\n",
215 __FILE__, __LINE__);
216 exit(1);
217 }
218 if(tng_util_vel_write_frequency_set(traj, step_save) != TNG_SUCCESS)
219 {
220 printf("Error setting writing frequency data. %s: %d\n",
221 __FILE__, __LINE__);
222 exit(1);
223 }
224 if(tng_util_force_write_frequency_set(traj, step_save) != TNG_SUCCESS)
225 {
226 printf("Error setting writing frequency data. %s: %d\n",
227 __FILE__, __LINE__);
228 exit(1);
229 }
230
231 /* Write the first frame of box shape, positions, velocities and forces */
232 if(tng_util_box_shape_write(traj, 0, box_shape) != TNG_SUCCESS)
233 {
234 printf("Error writing box shape. %s: %d\n",
235 __FILE__, __LINE__);
236 exit(1);
237 }
238 if(tng_util_pos_write(traj, 0, pos) != TNG_SUCCESS)
239 {
240 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
241 exit(1);
242 }
243 if(tng_util_vel_write(traj, 0, vel) != TNG_SUCCESS)
244 {
245 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
246 exit(1);
247 }
248 if(tng_util_force_write(traj, 0, force) != TNG_SUCCESS)
249 {
250 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
251 exit(1);
252 }
253
254 wtime = omp_get_wtime ( );
255
256 for ( step = 1; step < step_num; step++ )
257 {
258 compute ( np, nd, pos, vel, mass, force, &potential, &kinetic );
259
260 if ( step == step_print )
261 {
262 printf ( " %8d %14f %14f %14e\n", step, potential, kinetic,
263 ( potential + kinetic - e0 ) / e0 );
264 step_print_index++;
265 step_print = ( step_print_index * step_num ) / step_print_num;
266 }
267 if(step % step_save == 0)
268 {
269 /* Write box shape, positions, velocities and forces */
270 if(tng_util_box_shape_write(traj, step, box_shape) != TNG_SUCCESS)
271 {
272 printf("Error writing box shape. %s: %d\n",
273 __FILE__, __LINE__);
274 exit(1);
275 }
276 if(tng_util_pos_write(traj, step, pos) != TNG_SUCCESS)
277 {
278 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
279 break;
280 }
281 if(tng_util_vel_write(traj, step, vel) != TNG_SUCCESS)
282 {
283 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
284 break;
285 }
286 if(tng_util_force_write(traj, step, force) != TNG_SUCCESS)
287 {
288 printf("Error adding data. %s: %d\n", __FILE__, __LINE__);
289 break;
290 }
291 }
292 update ( np, nd, pos, vel, force, acc, mass, dt );
293 }
294 wtime = omp_get_wtime ( ) - wtime;
295
296 printf ( "\n" );
297 printf ( " Elapsed time for main computation:\n" );
298 printf ( " %f seconds.\n", wtime );
299
300 free ( acc );
301 free ( box );
302 free ( box_shape );
303 free ( force );
304 free ( pos );
305 free ( vel );
306
307 /* Close the TNG output. */
308 tng_util_trajectory_close(&traj);
309
310 printf ( "\n" );
311 printf ( "MD_OPENMP\n" );
312 printf ( " Normal end of execution.\n" );
313
314 printf ( "\n" );
315 timestamp ( );
316
317 return 0;
318 }
319 /******************************************************************************/
320
compute(int np,int nd,float pos[],float vel[],float mass,float f[],float * pot,float * kin)321 void compute ( int np, int nd, float pos[], float vel[],
322 float mass, float f[], float *pot, float *kin )
323
324 /******************************************************************************/
325 /*
326 Purpose:
327
328 COMPUTE computes the forces and energies.
329
330 Discussion:
331
332 The computation of forces and energies is fully parallel.
333
334 The potential function V(X) is a harmonic well which smoothly
335 saturates to a maximum value at PI/2:
336
337 v(x) = ( sin ( min ( x, PI2 ) ) )**2
338
339 The derivative of the potential is:
340
341 dv(x) = 2.0 * sin ( min ( x, PI2 ) ) * cos ( min ( x, PI2 ) )
342 = sin ( 2.0 * min ( x, PI2 ) )
343
344 Licensing:
345
346 This code is distributed under the GNU LGPL license.
347
348 Modified:
349
350 21 November 2007
351
352 Author:
353
354 Original FORTRAN77 version by Bill Magro.
355 C version by John Burkardt.
356
357 Parameters:
358
359 Input, int NP, the number of particles.
360
361 Input, int ND, the number of spatial dimensions.
362
363 Input, float POS[ND*NP], the position of each particle.
364
365 Input, float VEL[ND*NP], the velocity of each particle.
366
367 Input, float MASS, the mass of each particle.
368
369 Output, float F[ND*NP], the forces.
370
371 Output, float *POT, the total potential energy.
372
373 Output, float *KIN, the total kinetic energy.
374 */
375 {
376 float d;
377 float d2;
378 int i;
379 int j;
380 int k;
381 float ke;
382 float pe;
383 float PI2 = 3.141592653589793 / 2.0;
384 float rij[3];
385
386 pe = 0.0;
387 ke = 0.0;
388
389 # pragma omp parallel \
390 shared ( f, nd, np, pos, vel ) \
391 private ( i, j, k, rij, d, d2 )
392
393
394 # pragma omp for reduction ( + : pe, ke )
395 for ( k = 0; k < np; k++ )
396 {
397 /*
398 Compute the potential energy and forces.
399 */
400 for ( i = 0; i < nd; i++ )
401 {
402 f[i+k*nd] = 0.0;
403 }
404
405 for ( j = 0; j < np; j++ )
406 {
407 if ( k != j )
408 {
409 d = dist ( nd, pos+k*nd, pos+j*nd, rij );
410 /*
411 Attribute half of the potential energy to particle J.
412 */
413 if ( d < PI2 )
414 {
415 d2 = d;
416 }
417 else
418 {
419 d2 = PI2;
420 }
421
422 pe = pe + 0.5 * pow ( sin ( d2 ), 2 );
423
424 for ( i = 0; i < nd; i++ )
425 {
426 f[i+k*nd] = f[i+k*nd] - rij[i] * sin ( 2.0 * d2 ) / d;
427 }
428 }
429 }
430 /*
431 Compute the kinetic energy.
432 */
433 for ( i = 0; i < nd; i++ )
434 {
435 ke = ke + vel[i+k*nd] * vel[i+k*nd];
436 }
437 }
438
439 ke = ke * 0.5 * mass;
440
441 *pot = pe;
442 *kin = ke;
443
444 return;
445 }
446 /******************************************************************************/
447
dist(int nd,float r1[],float r2[],float dr[])448 float dist ( int nd, float r1[], float r2[], float dr[] )
449
450 /******************************************************************************/
451 /*
452 Purpose:
453
454 DIST computes the displacement (and its norm) between two particles.
455
456 Licensing:
457
458 This code is distributed under the GNU LGPL license.
459
460 Modified:
461
462 21 November 2007
463
464 Author:
465
466 Original FORTRAN77 version by Bill Magro.
467 C version by John Burkardt.
468
469 Parameters:
470
471 Input, int ND, the number of spatial dimensions.
472
473 Input, float R1[ND], R2[ND], the positions of the particles.
474
475 Output, float DR[ND], the displacement vector.
476
477 Output, float D, the Euclidean norm of the displacement.
478 */
479 {
480 float d;
481 int i;
482
483 d = 0.0;
484 for ( i = 0; i < nd; i++ )
485 {
486 dr[i] = r1[i] - r2[i];
487 d = d + dr[i] * dr[i];
488 }
489 d = sqrt ( d );
490
491 return d;
492 }
493 /******************************************************************************/
494
initialize(int np,int nd,float box[],int * seed,float pos[],float vel[],float acc[])495 void initialize ( int np, int nd, float box[], int *seed, float pos[],
496 float vel[], float acc[] )
497
498 /******************************************************************************/
499 /*
500 Purpose:
501
502 INITIALIZE initializes the positions, velocities, and accelerations.
503
504 Licensing:
505
506 This code is distributed under the GNU LGPL license.
507
508 Modified:
509
510 21 November 2007
511
512 Author:
513
514 Original FORTRAN77 version by Bill Magro.
515 C version by John Burkardt.
516
517 Parameters:
518
519 Input, int NP, the number of particles.
520
521 Input, int ND, the number of spatial dimensions.
522
523 Input, float BOX[ND], specifies the maximum position
524 of particles in each dimension.
525
526 Input, int *SEED, a seed for the random number generator.
527
528 Output, float POS[ND*NP], the position of each particle.
529
530 Output, float VEL[ND*NP], the velocity of each particle.
531
532 Output, float ACC[ND*NP], the acceleration of each particle.
533 */
534 {
535 int i;
536 int j;
537 /*
538 Give the particles random positions within the box.
539 */
540 for ( i = 0; i < nd; i++ )
541 {
542 for ( j = 0; j < np; j++ )
543 {
544 pos[i+j*nd] = box[i] * r8_uniform_01 ( seed );
545 }
546 }
547
548 for ( j = 0; j < np; j++ )
549 {
550 for ( i = 0; i < nd; i++ )
551 {
552 vel[i+j*nd] = 0.0;
553 }
554 }
555 for ( j = 0; j < np; j++ )
556 {
557 for ( i = 0; i < nd; i++ )
558 {
559 acc[i+j*nd] = 0.0;
560 }
561 }
562 return;
563 }
564 /******************************************************************************/
565
r8_uniform_01(int * seed)566 float r8_uniform_01 ( int *seed )
567
568 /******************************************************************************/
569 /*
570 Purpose:
571
572 R8_UNIFORM_01 is a unit pseudorandom R8.
573
574 Discussion:
575
576 This routine implements the recursion
577
578 seed = 16807 * seed mod ( 2**31 - 1 )
579 unif = seed / ( 2**31 - 1 )
580
581 The integer arithmetic never requires more than 32 bits,
582 including a sign bit.
583
584 Licensing:
585
586 This code is distributed under the GNU LGPL license.
587
588 Modified:
589
590 11 August 2004
591
592 Author:
593
594 John Burkardt
595
596 Reference:
597
598 Paul Bratley, Bennett Fox, Linus Schrage,
599 A Guide to Simulation,
600 Springer Verlag, pages 201-202, 1983.
601
602 Bennett Fox,
603 Algorithm 647:
604 Implementation and Relative Efficiency of Quasirandom
605 Sequence Generators,
606 ACM Transactions on Mathematical Software,
607 Volume 12, Number 4, pages 362-376, 1986.
608
609 Parameters:
610
611 Input/output, int *SEED, a seed for the random number generator.
612
613 Output, float R8_UNIFORM_01, a new pseudorandom variate, strictly between
614 0 and 1.
615 */
616 {
617 int k;
618 float r;
619
620 k = *seed / 127773;
621
622 *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
623
624 if ( *seed < 0 )
625 {
626 *seed = *seed + 2147483647;
627 }
628
629 r = ( float ) ( *seed ) * 4.656612875E-10;
630
631 return r;
632 }
633 /******************************************************************************/
634
timestamp(void)635 void timestamp ( void )
636
637 /******************************************************************************/
638 /*
639 Purpose:
640
641 TIMESTAMP prints the current YMDHMS date as a time stamp.
642
643 Example:
644
645 31 May 2001 09:45:54 AM
646
647 Licensing:
648
649 This code is distributed under the GNU LGPL license.
650
651 Modified:
652
653 24 September 2003
654
655 Author:
656
657 John Burkardt
658
659 Parameters:
660
661 None
662 */
663 {
664 # define TIME_SIZE 40
665
666 static char time_buffer[TIME_SIZE];
667 const struct tm *tm;
668 time_t now;
669
670 now = time ( NULL );
671 tm = localtime ( &now );
672
673 strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
674
675 printf ( "%s\n", time_buffer );
676
677 return;
678 # undef TIME_SIZE
679 }
680 /******************************************************************************/
681
update(int np,int nd,float pos[],float vel[],float f[],float acc[],float mass,float dt)682 void update ( int np, int nd, float pos[], float vel[], float f[],
683 float acc[], float mass, float dt )
684
685 /******************************************************************************/
686 /*
687 Purpose:
688
689 UPDATE updates positions, velocities and accelerations.
690
691 Discussion:
692
693 The time integration is fully parallel.
694
695 A velocity Verlet algorithm is used for the updating.
696
697 x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt * dt
698 v(t+dt) = v(t) + 0.5 * ( a(t) + a(t+dt) ) * dt
699 a(t+dt) = f(t) / m
700
701 Licensing:
702
703 This code is distributed under the GNU LGPL license.
704
705 Modified:
706
707 17 April 2009
708
709 Author:
710
711 Original FORTRAN77 version by Bill Magro.
712 C version by John Burkardt.
713
714 Parameters:
715
716 Input, int NP, the number of particles.
717
718 Input, int ND, the number of spatial dimensions.
719
720 Input/output, float POS[ND*NP], the position of each particle.
721
722 Input/output, float VEL[ND*NP], the velocity of each particle.
723
724 Input, float F[ND*NP], the force on each particle.
725
726 Input/output, float ACC[ND*NP], the acceleration of each particle.
727
728 Input, float MASS, the mass of each particle.
729
730 Input, float DT, the time step.
731 */
732 {
733 int i;
734 int j;
735 float rmass;
736
737 rmass = 1.0 / mass;
738
739 # pragma omp parallel \
740 shared ( acc, dt, f, nd, np, pos, rmass, vel ) \
741 private ( i, j )
742
743 # pragma omp for
744 for ( j = 0; j < np; j++ )
745 {
746 for ( i = 0; i < nd; i++ )
747 {
748 pos[i+j*nd] = pos[i+j*nd] + vel[i+j*nd] * dt + 0.5 * acc[i+j*nd] * dt * dt;
749 vel[i+j*nd] = vel[i+j*nd] + 0.5 * dt * ( f[i+j*nd] * rmass + acc[i+j*nd] );
750 acc[i+j*nd] = f[i+j*nd] * rmass;
751 }
752 }
753
754
755 return;
756 }
757
758 #endif
759