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