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     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
 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
 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
 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
 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
 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
 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