src/external/tng_io/src/tests/md_openmp_util.c

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