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