src/external/tng_io/src/tests/md_openmp.f

   1       program main
   2
   3 c*********************************************************************72
   4 c
   5 cc MAIN is the main program for MD_OPENMP.
   6 c
   7 c  Discussion:
   8 c
   9 c    The program implements a simple molecular dynamics simulation.
  10 c
  11 c    The program uses Open MP directives to allow parallel computation.
  12 c
  13 c    The velocity Verlet time integration scheme is used.
  14 c
  15 c    The particles interact with a central pair potential.
  16 c
  17 c    Output of the program is saved in the TNG format, which is why this
  18 c    code is included in the TNG API release.
  19 c
  20 c  Licensing:
  21 c
  22 c    This code is distributed under the GNU LGPL license.
  23 c
  24 c  Modified:
  25 c
  26 c    8 Jan 2013
  27 c
  28 c  Author:
  29 c
  30 c    Original FORTRAN90 version by Bill Magro.
  31 c    FORTRAN77 version by John Burkardt.
  32 c    TNG trajectory output by Magnus Lundborg.
  33 c
  34 c  Parameters:
  35 c
  36 c    None
  37 c
  38       implicit none
  39
  40       include 'omp_lib.h'
  41
  42       integer nd
  43       parameter ( nd = 3 )
  44       integer np
  45       parameter ( np = 250 )
  46       integer step_num
  47       parameter ( step_num = 1000 )
  48
  49       double precision acc(nd,np)
  50       double precision box(nd)
  51       double precision box_shape(9)
  52       double precision dt
  53       parameter ( dt = 0.0001D+00 )
  54       double precision e0
  55       double precision force(nd,np)
  56       integer i
  57       integer id
  58       double precision kinetic
  59       double precision mass
  60       parameter ( mass = 1.0D+00 )
  61       double precision pos(nd,np)
  62       double precision potential
  63       integer proc_num
  64       integer seed
  65       integer step
  66       integer step_print
  67       integer step_print_index
  68       integer step_print_num
  69       integer step_save
  70       integer*8 sparse_save
  71       integer thread_num
  72       double precision vel(nd,np)
  73       double precision wtime
  74
  75 c
  76 c  Cray pointers are not standard fortran 77, but must be used to allocate
  77 c  memory properly.
  78 c
  79       pointer (traj_p, traj)
  80       pointer (molecule_p, molecule)
  81       pointer (chain_p, chain)
  82       pointer (residue_p, residue)
  83       pointer (atom_p, atom)
  84       byte traj
  85       byte molecule
  86       byte chain
  87       byte residue
  88       byte atom
  89
  90 c
  91 c  The TNG functions expect 64 bit integers
  92 c
  93       integer*8 n_frames_per_frame_set
  94       integer*8 frames_saved_cnt
  95       integer*8 tng_n_particles
  96
  97 c
  98 c  These constants are also defined in tng_io.h, but need to
  99 c  set in fortran as well. This can be copied to any fortran
 100 c  source code using the tng_io library.
 101 c
 102       integer*8 TNG_UNCOMPRESSED
 103       parameter ( TNG_UNCOMPRESSED = 0)
 104       integer TNG_NON_TRAJECTORY_BLOCK
 105       parameter ( TNG_NON_TRAJECTORY_BLOCK = 0)
 106       integer TNG_TRAJECTORY_BLOCK
 107       parameter ( TNG_TRAJECTORY_BLOCK = 1)
 108       integer*8 TNG_GENERAL_INFO
 109       parameter ( TNG_GENERAL_INFO = 0 )
 110       integer*8 TNG_MOLECULES
 111       parameter ( TNG_MOLECULES = 1 )
 112       integer*8 TNG_TRAJECTORY_FRAME_SET
 113       parameter ( TNG_TRAJECTORY_FRAME_SET = 2 )
 114       integer*8 TNG_PARTICLE_MAPPING
 115       parameter ( TNG_PARTICLE_MAPPING = 3 )
 116       integer*8 TNG_TRAJ_BOX_SHAPE
 117       parameter ( TNG_TRAJ_BOX_SHAPE = 10000 )
 118       integer*8 TNG_TRAJ_POSITIONS
 119       parameter ( TNG_TRAJ_POSITIONS = 10001 )
 120       integer*8 TNG_TRAJ_VELOCITIES
 121       parameter ( TNG_TRAJ_VELOCITIES = 10002 )
 122       integer*8 TNG_TRAJ_FORCES
 123       parameter ( TNG_TRAJ_FORCES = 10003 )
 124       integer TNG_SKIP_HASH
 125       parameter ( TNG_SKIP_HASH = 0 )
 126       integer TNG_USE_HASH
 127       parameter ( TNG_USE_HASH = 1 )
 128       integer TNG_CHAR_DATA
 129       parameter ( TNG_CHAR_DATA = 0 )
 130       integer TNG_INT_DATA
 131       parameter ( TNG_INT_DATA = 1 )
 132       integer TNG_FLOAT_DATA
 133       parameter ( TNG_FLOAT_DATA = 2 )
 134       integer TNG_DOUBLE_DATA
 135       parameter ( TNG_DOUBLE_DATA = 3 )
 136
 137       call timestamp ( )
 138
 139       write ( *, '(a)' ) ' '
 140       write ( *, '(a)' ) 'MD_OPENMP'
 141       write ( *, '(a)' ) '  FORTRAN77/OpenMP version'
 142       write ( *, '(a)' ) ' '
 143       write ( *, '(a)' ) '  A molecular dynamics program.'
 144       write ( *, '(a)' ) ' '
 145       write ( *, '(a,i8)' )
 146      &  '  NP, the number of particles in the simulation is ', np
 147       write ( *, '(a,i8)' )
 148      &  '  STEP_NUM, the number of time steps, is ', step_num
 149       write ( *, '(a,g14.6)' )
 150      &  '  DT, the size of each time step, is ', dt
 151       write ( *, '(a)' ) ' '
 152       write ( *, '(a,i8)' )
 153      &  '  The number of processors = ', omp_get_num_procs ( )
 154       write ( *, '(a,i8)' )
 155      &  '  The number of threads    = ', omp_get_max_threads ( )
 156
 157       write ( *, '(a)' ) ' '
 158       write ( *, '(a)' ) '  Initializing trajectory storage.'
 159       call tng_trajectory_init(traj_p)
 160
 161 c
 162 c  N.B. The TNG output file should be modified according to needs
 163 c
 164       call tng_output_file_set(traj, "/tmp/tng_md_out_f77.tng")
 165
 166       write ( *, '(a)' ) '  Creating molecules in trajectory.'
 167       tng_n_particles = np
 168       call tng_molecule_add(traj, "water", molecule_p)
 169       call tng_molecule_chain_add(traj, molecule, "W", chain_p)
 170       call tng_chain_residue_add(traj, chain, "WAT", residue_p)
 171       call tng_residue_atom_add(traj, residue, "O", "O", atom_p)
 172       call tng_molecule_cnt_set(traj, molecule, tng_n_particles)
 173 c
 174 c  Set the dimensions of the box.
 175 c
 176       do i = 1, 9
 177         box_shape(i) = 0.0
 178       end do
 179       do i = 1, nd
 180         box(i) = 10.0D+00
 181         box_shape(i*nd + i) = box(i)
 182       end do
 183
 184 c
 185 c  Add the box shape data block
 186 c
 187       call tng_data_block_add(traj, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE",
 188      &  TNG_DOUBLE_DATA, TNG_NON_TRAJECTORY_BLOCK, int(1, 8),
 189      &  int(9, 8), int(1, 8), TNG_UNCOMPRESSED, box_shape)
 190
 191 c
 192 c  Write the file headers
 193 c
 194       call tng_file_headers_write(traj, TNG_USE_HASH)
 195
 196 c
 197 c  Set initial positions, velocities, and accelerations.
 198 c
 199       write ( *, '(a)' )
 200      &  '  Initializing positions, velocities, and accelerations.'
 201       seed = 123456789
 202       call initialize ( np, nd, box, seed, pos, vel, acc )
 203 c
 204 c  Compute the forces and energies.
 205 c
 206       write ( *, '(a)' ) ' '
 207       write ( *, '(a)' ) '  Computing initial forces and energies.'
 208
 209       call compute ( np, nd, pos, vel, mass, force, potential,
 210      &  kinetic )
 211
 212       e0 = potential + kinetic
 213
 214 c
 215 c  Saving frequency
 216 c
 217       step_save = 5
 218
 219       step_print = 0
 220       step_print_index = 0
 221       step_print_num = 10
 222       sparse_save = 10
 223
 224       frames_saved_cnt = 0
 225
 226 c
 227 c  This is the main time stepping loop.
 228 c
 229       write ( *, '(a,i4,a)' ) '  Every', step_save,
 230      &  ' steps particle positions, velocities and forces are'
 231       write ( *, '(a)' ) '  saved to a TNG trajectory file.'
 232       write ( *, '(a)' )
 233       write ( *, '(a)' )
 234      &  '  At each step, we report the potential and kinetic energies.'
 235       write ( *, '(a)' )
 236      &  '  The sum of these energies should be a constant.'
 237       write ( *, '(a)' )
 238      &  '  As an accuracy check, we also print the relative error'
 239       write ( *, '(a)' ) '  in the total energy.'
 240       write ( *, '(a)' ) ' '
 241       write ( *, '(a)' )
 242      &  '      Step      Potential       Kinetic        (P+K-E0)/E0'
 243       write ( *, '(a)' )
 244      &  '                Energy P        Energy K       ' //
 245      &  'Relative Energy Error'
 246       write ( *, '(a)' ) ' '
 247
 248       step = 0
 249       write ( *, '(2x,i8,2x,g14.6,2x,g14.6,2x,g14.6)' )
 250      &  step, potential, kinetic, ( potential + kinetic - e0 ) / e0
 251       step_print_index = step_print_index + 1
 252       step_print = ( step_print_index * step_num ) / step_print_num
 253
 254 c
 255 c  Create a frame set for writing data
 256 c
 257       call tng_num_frames_per_frame_set_get(traj,
 258      &  n_frames_per_frame_set)
 259       call tng_frame_set_new(traj, int(0, 8), n_frames_per_frame_set)
 260
 261 c
 262 c  Add empty data blocks.
 263 c
 264       call tng_particle_data_block_add(traj, TNG_TRAJ_POSITIONS,
 265      &  "POSITIONS", TNG_DOUBLE_DATA, TNG_TRAJECTORY_BLOCK,
 266      &  n_frames_per_frame_set, int(3, 8), int(1, 8), int(0, 8),
 267      &  tng_n_particles, TNG_UNCOMPRESSED, %VAL(int(0, 8)))
 268
 269       call tng_particle_data_block_add(traj, TNG_TRAJ_VELOCITIES,
 270      &  "VELOCITIES", TNG_DOUBLE_DATA, TNG_TRAJECTORY_BLOCK,
 271      &  n_frames_per_frame_set, int(3, 8), int(1, 8), int(0, 8),
 272      &  tng_n_particles, TNG_UNCOMPRESSED, %VAL(int(0, 8)))
 273
 274       call tng_particle_data_block_add(traj, TNG_TRAJ_FORCES,
 275      &  "FORCES", TNG_DOUBLE_DATA, TNG_TRAJECTORY_BLOCK,
 276      &  n_frames_per_frame_set, int(3, 8), int(1, 8), int(0, 8),
 277      &  tng_n_particles, TNG_UNCOMPRESSED, %VAL(int(0, 8)))
 278
 279 c
 280 c  The potential energy data block is saved sparsely.
 281 c
 282       call tng_data_block_add(traj, int(10101, 8),
 283      &  "POTENTIAL ENERGY", TNG_DOUBLE_DATA, TNG_TRAJECTORY_BLOCK,
 284      &  n_frames_per_frame_set, int(1, 8), sparse_save,
 285      &  TNG_UNCOMPRESSED, %VAL(int(0, 8)))
 286
 287
 288 c
 289 c  Write the frame set to disk
 290 c
 291       call tng_frame_set_write(traj, TNG_USE_HASH)
 292
 293       wtime = omp_get_wtime ( )
 294
 295       do step = 1, step_num
 296
 297         call compute ( np, nd, pos, vel, mass, force, potential,
 298      &    kinetic )
 299
 300         if ( step .eq. step_print ) then
 301
 302           write ( *, '(2x,i8,2x,g14.6,2x,g14.6,2x,g14.6)' )
 303      &      step, potential, kinetic, ( potential + kinetic - e0 ) / e0
 304
 305           step_print_index = step_print_index + 1
 306           step_print = ( step_print_index * step_num ) / step_print_num
 307
 308         end if
 309
 310 c
 311 c  Output to TNG file at regular intervals, specified by step_save
 312 c
 313         if ( step_save .EQ. 0 .OR. mod(step, step_save) .EQ. 0 ) then
 314           call tng_frame_particle_data_write(traj, frames_saved_cnt,
 315      &    TNG_TRAJ_POSITIONS, int(0, 8), tng_n_particles, pos,
 316      &    TNG_USE_HASH)
 317           call tng_frame_particle_data_write(traj, frames_saved_cnt,
 318      &    TNG_TRAJ_VELOCITIES, int(0, 8), tng_n_particles, vel,
 319      &    TNG_USE_HASH)
 320           call tng_frame_particle_data_write(traj, frames_saved_cnt,
 321      &    TNG_TRAJ_FORCES, int(0, 8), tng_n_particles, force,
 322      &    TNG_USE_HASH)
 323           frames_saved_cnt = frames_saved_cnt + 1
 324
 325           if (mod(step, step_save * sparse_save) .EQ. 0) then
 326             call tng_frame_data_write(traj, frames_saved_cnt,
 327      &      int(10101, 8), potential, TNG_USE_HASH)
 328           end if
 329
 330         end if
 331
 332         call update ( np, nd, pos, vel, force, acc, mass, dt )
 333
 334       end do
 335
 336       wtime = omp_get_wtime ( ) - wtime
 337
 338       write ( *, '(a)' ) ' '
 339       write ( *, '(a)' )
 340      &  '  Elapsed time for main computation:'
 341       write ( *, '(2x,g14.6,a)' ) wtime, ' seconds'
 342 c
 343 c  Terminate.
 344 c
 345       call tng_trajectory_destroy(traj_p)
 346
 347       write ( *, '(a)' ) ' '
 348       write ( *, '(a)' ) 'MD_OPENMP'
 349       write ( *, '(a)' ) '  Normal end of execution.'
 350
 351       write ( *, '(a)' ) ' '
 352       call timestamp ( )
 353
 354       stop
 355       end
 356       subroutine compute ( np, nd, pos, vel, mass, f, pot, kin )
 357
 358 c*********************************************************************72
 359 c
 360 cc COMPUTE computes the forces and energies.
 361 c
 362 c  Discussion:
 363 c
 364 c    The computation of forces and energies is fully parallel.
 365 c
 366 c  Licensing:
 367 c
 368 c    This code is distributed under the GNU LGPL license.
 369 c
 370 c  Modified:
 371 c
 372 c    31 July 2009
 373 c
 374 c  Author:
 375 c
 376 c    Original FORTRAN90 version by Bill Magro.
 377 c    FORTRAN77 version by John Burkardt.
 378 c
 379 c  Parameters:
 380 c
 381 c    Input, integer NP, the number of particles.
 382 c
 383 c    Input, integer ND, the number of spatial dimensions.
 384 c
 385 c    Input, double precision POS(ND,NP), the position of each particle.
 386 c
 387 c    Input, double precision VEL(ND,NP), the velocity of each particle.
 388 c
 389 c    Input, double precision MASS, the mass of each particle.
 390 c
 391       implicit none
 392
 393       integer np
 394       integer nd
 395
 396       double precision d
 397       double precision d2
 398       double precision dv
 399       double precision f(nd,np)
 400       integer i
 401       integer j
 402       integer k
 403       double precision kin
 404       double precision mass
 405       double precision PI2
 406       parameter ( PI2 = 3.141592653589793D+00 / 2.0D+00 )
 407       double precision pos(nd,np)
 408       double precision pot
 409       double precision rij(nd)
 410       double precision v
 411       double precision vel(nd,np)
 412
 413       pot = 0.0D+00
 414       kin = 0.0D+00
 415
 416 c$omp parallel
 417 c$omp& shared ( f, nd, np, pos, vel )
 418 c$omp& private ( d, d2, i, j, k, rij )
 419
 420 c$omp do reduction ( + : pot, kin )
 421       do i = 1, np
 422 c
 423 c  Compute the potential energy and forces.
 424 c
 425         do k = 1, nd
 426           f(k,i) = 0.0D+00
 427         end do
 428
 429         do j = 1, np
 430
 431           if ( i .ne. j ) then
 432
 433             call dist ( nd, pos(1,i), pos(1,j), rij, d )
 434 c
 435 c  Attribute half of the potential energy to particle J.
 436 c
 437             d2 = min ( d, pi2 )
 438
 439             pot = pot + 0.5D+00 * ( sin ( d2 ) )**2
 440
 441             do k = 1, nd
 442               f(k,i) = f(k,i) - rij(k) * sin ( 2.0D+00 * d2 ) / d
 443             end do
 444
 445           end if
 446
 447         end do
 448 c
 449 c  Compute the kinetic energy.
 450 c
 451         do k = 1, nd
 452           kin = kin + vel(k,i)**2
 453         end do
 454
 455       end do
 456 c$omp end do
 457
 458 c$omp end parallel
 459
 460       kin = kin * 0.5D+00 * mass
 461
 462       return
 463       end
 464       subroutine dist ( nd, r1, r2, dr, d )
 465
 466 c*********************************************************************72
 467 c
 468 cc DIST computes the displacement (and its norm) between two particles.
 469 c
 470 c  Licensing:
 471 c
 472 c    This code is distributed under the GNU LGPL license.
 473 c
 474 c  Modified:
 475 c
 476 c    13 November 2007
 477 c
 478 c  Author:
 479 c
 480 c    Original FORTRAN90 version by Bill Magro.
 481 c    FORTRAN77 version by John Burkardt.
 482 c
 483 c  Parameters:
 484 c
 485 c    Input, integer ND, the number of spatial dimensions.
 486 c
 487 c    Input, double precision R1(ND), R2(ND), the positions of the particles.
 488 c
 489 c    Output, double precision DR(ND), the displacement vector.
 490 c
 491 c    Output, double precision D, the Euclidean norm of the displacement.
 492 c
 493       implicit none
 494
 495       integer nd
 496
 497       double precision d
 498       double precision dr(nd)
 499       integer i
 500       double precision r1(nd)
 501       double precision r2(nd)
 502
 503       do i = 1, nd
 504         dr(i) = r1(i) - r2(i)
 505       end do
 506
 507       d = 0.0D+00
 508       do i = 1, nd
 509         d = d + dr(i)**2
 510       end do
 511       d = sqrt ( d )
 512
 513       return
 514       end
 515       subroutine initialize ( np, nd, box, seed, pos, vel, acc )
 516
 517 c*********************************************************************72
 518 c
 519 cc INITIALIZE initializes the positions, velocities, and accelerations.
 520 c
 521 c  Licensing:
 522 c
 523 c    This code is distributed under the GNU LGPL license.
 524 c
 525 c  Modified:
 526 c
 527 c    13 November 2007
 528 c
 529 c  Author:
 530 c
 531 c    Original FORTRAN90 version by Bill Magro.
 532 c    FORTRAN77 version by John Burkardt.
 533 c
 534 c  Parameters:
 535 c
 536 c    Input, integer NP, the number of particles.
 537 c
 538 c    Input, integer ND, the number of spatial dimensions.
 539 c
 540 c    Input, double precision BOX(ND), specifies the maximum position
 541 c    of particles in each dimension.
 542 c
 543 c    Input/output, integer SEED, a seed for the random number generator.
 544 c
 545 c    Output, double precision POS(ND,NP), the position of each particle.
 546 c
 547 c    Output, double precision VEL(ND,NP), the velocity of each particle.
 548 c
 549 c    Output, double precision ACC(ND,NP), the acceleration of each particle.
 550 c
 551       implicit none
 552
 553       integer np
 554       integer nd
 555
 556       double precision acc(nd,np)
 557       double precision box(nd)
 558       integer i
 559       integer j
 560       double precision pos(nd,np)
 561       double precision r8_uniform_01
 562       integer seed
 563       double precision vel(nd,np)
 564 c
 565 c  Give the particles random positions within the box.
 566 c
 567       do i = 1, nd
 568         do j = 1, np
 569           pos(i,j) = r8_uniform_01 ( seed )
 570         end do
 571       end do
 572
 573 c$omp parallel
 574 c$omp& shared ( acc, box, nd, np, pos, vel )
 575 c$omp& private ( i, j )
 576
 577 c$omp do
 578       do j = 1, np
 579         do i = 1, nd
 580           pos(i,j) = box(i) * pos(i,j)
 581           vel(i,j) = 0.0D+00
 582           acc(i,j) = 0.0D+00
 583         end do
 584       end do
 585 c$omp end do
 586
 587 c$omp end parallel
 588
 589       return
 590       end
 591       function r8_uniform_01 ( seed )
 592
 593 c*********************************************************************72
 594 c
 595 cc R8_UNIFORM_01 returns a unit pseudorandom R8.
 596 c
 597 c  Discussion:
 598 c
 599 c    This routine implements the recursion
 600 c
 601 c      seed = 16807 * seed mod ( 2**31 - 1 )
 602 c      r8_uniform_01 = seed / ( 2**31 - 1 )
 603 c
 604 c    The integer arithmetic never requires more than 32 bits,
 605 c    including a sign bit.
 606 c
 607 c    If the initial seed is 12345, then the first three computations are
 608 c
 609 c      Input     Output      R8_UNIFORM_01
 610 c      SEED      SEED
 611 c
 612 c         12345   207482415  0.096616
 613 c     207482415  1790989824  0.833995
 614 c    1790989824  2035175616  0.947702
 615 c
 616 c  Licensing:
 617 c
 618 c    This code is distributed under the GNU LGPL license.
 619 c
 620 c  Modified:
 621 c
 622 c    11 August 2004
 623 c
 624 c  Author:
 625 c
 626 c    John Burkardt
 627 c
 628 c  Reference:
 629 c
 630 c    Paul Bratley, Bennett Fox, Linus Schrage,
 631 c    A Guide to Simulation,
 632 c    Springer Verlag, pages 201-202, 1983.
 633 c
 634 c    Pierre L'Ecuyer,
 635 c    Random Number Generation,
 636 c    in Handbook of Simulation,
 637 c    edited by Jerry Banks,
 638 c    Wiley Interscience, page 95, 1998.
 639 c
 640 c    Bennett Fox,
 641 c    Algorithm 647:
 642 c    Implementation and Relative Efficiency of Quasirandom
 643 c    Sequence Generators,
 644 c    ACM Transactions on Mathematical Software,
 645 c    Volume 12, Number 4, pages 362-376, 1986.
 646 c
 647 c    Peter Lewis, Allen Goodman, James Miller,
 648 c    A Pseudo-Random Number Generator for the System/360,
 649 c    IBM Systems Journal,
 650 c    Volume 8, pages 136-143, 1969.
 651 c
 652 c  Parameters:
 653 c
 654 c    Input/output, integer SEED, the "seed" value, which should NOT be 0.
 655 c    On output, SEED has been updated.
 656 c
 657 c    Output, double precision R8_UNIFORM_01, a new pseudorandom variate,
 658 c    strictly between 0 and 1.
 659 c
 660       implicit none
 661
 662       double precision r8_uniform_01
 663       integer k
 664       integer seed
 665
 666       if ( seed .eq. 0 ) then
 667         write ( *, '(a)' ) ' '
 668         write ( *, '(a)' ) 'R8_UNIFORM_01 - Fatal error!'
 669         write ( *, '(a)' ) '  Input value of SEED = 0.'
 670         stop
 671       end if
 672
 673       k = seed / 127773
 674
 675       seed = 16807 * ( seed - k * 127773 ) - k * 2836
 676
 677       if ( seed .lt. 0 ) then
 678         seed = seed + 2147483647
 679       end if
 680 c
 681 c  Although SEED can be represented exactly as a 32 bit integer,
 682 c  it generally cannot be represented exactly as a 32 bit real number!
 683 c
 684       r8_uniform_01 = dble ( seed ) * 4.656612875D-10
 685
 686       return
 687       end
 688       subroutine timestamp ( )
 689
 690 c*********************************************************************72
 691 c
 692 cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
 693 c
 694 c  Licensing:
 695 c
 696 c    This code is distributed under the GNU LGPL license.
 697 c
 698 c  Modified:
 699 c
 700 c    12 January 2007
 701 c
 702 c  Author:
 703 c
 704 c    John Burkardt
 705 c
 706 c  Parameters:
 707 c
 708 c    None
 709 c
 710       implicit none
 711
 712       character * ( 8 ) ampm
 713       integer d
 714       character * ( 8 ) date
 715       integer h
 716       integer m
 717       integer mm
 718       character * ( 9 ) month(12)
 719       integer n
 720       integer s
 721       character * ( 10 ) time
 722       integer y
 723
 724       save month
 725
 726       data month /
 727      &  'January  ', 'February ', 'March    ', 'April    ',
 728      &  'May      ', 'June     ', 'July     ', 'August   ',
 729      &  'September', 'October  ', 'November ', 'December ' /
 730
 731       call date_and_time ( date, time )
 732
 733       read ( date, '(i4,i2,i2)' ) y, m, d
 734       read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm
 735
 736       if ( h .lt. 12 ) then
 737         ampm = 'AM'
 738       else if ( h .eq. 12 ) then
 739         if ( n .eq. 0 .and. s .eq. 0 ) then
 740           ampm = 'Noon'
 741         else
 742           ampm = 'PM'
 743         end if
 744       else
 745         h = h - 12
 746         if ( h .lt. 12 ) then
 747           ampm = 'PM'
 748         else if ( h .eq. 12 ) then
 749           if ( n .eq. 0 .and. s .eq. 0 ) then
 750             ampm = 'Midnight'
 751           else
 752             ampm = 'AM'
 753           end if
 754         end if
 755       end if
 756
 757       write ( *,
 758      &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' )
 759      &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm
 760
 761       return
 762       end
 763       subroutine update ( np, nd, pos, vel, f, acc, mass, dt )
 764
 765 c*********************************************************************72
 766 c
 767 cc UPDATE performs the time integration, using a velocity Verlet algorithm.
 768 c
 769 c  Discussion:
 770 c
 771 c    The time integration is fully parallel.
 772 c
 773 c  Licensing:
 774 c
 775 c    This code is distributed under the GNU LGPL license.
 776 c
 777 c  Modified:
 778 c
 779 c    13 November 2007
 780 c
 781 c  Author:
 782 c
 783 c    Original FORTRAN90 version by Bill Magro.
 784 c    FORTRAN77 version by John Burkardt.
 785 c
 786 c  Parameters:
 787 c
 788 c    Input, integer NP, the number of particles.
 789 c
 790 c    Input, integer ND, the number of spatial dimensions.
 791 c
 792 c    Input/output, double precision POS(ND,NP), the position of each particle.
 793 c
 794 c    Input/output, double precision VEL(ND,NP), the velocity of each particle.
 795 c
 796 c    Input, double precision MASS, the mass of each particle.
 797 c
 798 c    Input/output, double precision ACC(ND,NP), the acceleration of each
 799 c    particle.
 800 c
 801       implicit none
 802
 803       integer np
 804       integer nd
 805
 806       double precision acc(nd,np)
 807       double precision dt
 808       double precision f(nd,np)
 809       integer i
 810       integer j
 811       double precision mass
 812       double precision pos(nd,np)
 813       double precision rmass
 814       double precision vel(nd,np)
 815
 816       rmass = 1.0D+00 / mass
 817
 818 c$omp parallel
 819 c$omp& shared ( acc, dt, f, nd, np, pos, rmass, vel )
 820 c$omp& private ( i, j )
 821
 822 c$omp do
 823       do j = 1, np
 824         do i = 1, nd
 825
 826           pos(i,j) = pos(i,j)
 827      &      + vel(i,j) * dt + 0.5D+00 * acc(i,j) * dt * dt
 828
 829           vel(i,j) = vel(i,j)
 830      &      + 0.5D+00 * dt * ( f(i,j) * rmass + acc(i,j) )
 831
 832           acc(i,j) = f(i,j) * rmass
 833
 834         end do
 835       end do
 836 c$omp end do
 837
 838 c$omp end parallel
 839
 840       return
 841       end