program md
!
!*******************************************************************************
!
!! MD implements a simple molecular dynamics simulation.
!
!
!  Discussion:
!
!    The program uses Open MP directives to allow parallel computation.
!
!    using the velocity Verlet time integration scheme. The particles
!    interact with a central pair potential.
!
!  Modified:
!
!    17 March 2002
!
!  Author:
!
!    Bill Magro, 
!    Kuck and Associates, Inc. (KAI), 1998
!
! THIS PROGRAM USES THE FORTRAN90 RANDOM_NUMBER FUNCTION AND ARRAY 
!   SYNTAX
!
!
!  Parameters:
!
  implicit none
!
  integer, parameter :: nd = 3
  integer, parameter :: np = 500
  integer, parameter :: nsteps = 1000
!
  double precision accel(nd,np)
  double precision box(nd)   ! dimensions of the simulation box
  double precision dt          ! time step
  double precision e0
  double precision force(nd,np)
  integer i
  double precision kinetic
  double precision mass        ! mass of the particles
  double precision pos(nd,np)
  double precision potential
  double precision velocity(nd,np)

  parameter(mass=1.0,dt=1.0e-4)
!
  box(1:nd) = 10.0D+00
!
!  Set initial positions, velocities, and accelerations.
!
  call initialize ( np, nd, box, pos, velocity, accel )
!
!  Compute the forces and energies.
!
  call compute ( np, nd, box, pos, velocity, mass, &
    force, potential, kinetic )

  e0 = potential + kinetic
!
!  This is the main time stepping loop.
!
  do i = 1, nsteps

    call compute ( np, nd, box, pos, velocity, mass, &
      force, potential, kinetic )

    write ( *, '(3g14.6)' ) potential, kinetic, 
      ( potential + kinetic - e0 ) / e0

    call update ( np, nd, pos, velocity, force, accel, mass, dt )

  end do

  stop
end
subroutine compute ( np, nd, box, pos, vel, mass, f, pot, kin )
!
!*******************************************************************************
!
!! COMPUTE computes the forces and energies, given positions, masses,
! and velocities
!
!
!  Discussion:
!
!    The computation of forces and energies is fully parallel.
!
!  Modified:
!
!    17 March 2002
!
!  Author:
!
!    Bill Magro, 
!    Kuck and Associates, Inc. (KAI), 1998
!
!  Parameters:
!
!    Input, integer NP, the number of particles.
!
!    Input, integer ND, the number of spatial dimensions.
!
!    Input, double precision BOX(ND), specifies the maximum position
!    of particles in each dimension.
!
!    Input, double precision POS(ND,NP), the position of each particle.
!
  implicit none

  integer np
  integer nd
!
  double precision box(nd)
  double precision d
  double precision dv
  double precision f(nd,np)
  integer i
  integer j
  integer k
  double precision kin
  double precision mass
  double precision, parameter :: PI2 = 3.14159265d0/2.0d0
  double precision pos(nd,np)
  double precision pot
  double precision rij(nd)
  double precision v
  double precision vel(nd,np)
  double precision x


  double precision dotr8
  external dotr8
!
! statement function for the pair potential and its derivative
! This potential is a harmonic well which smoothly saturates to a
! maximum value at PI/2.
!
  v(x) = sin ( min ( x, PI2 ) )**2.
  dv(x) = 2.0 * sin ( min (x,PI2) ) * cos(min(x,PI2) )

  pot = 0.0
  kin = 0.0

!$omp  parallel do
!$omp& default(shared)
!$omp& private(i,j,k,rij,d)
!$omp& reduction(+ : pot, kin)
  do i = 1, np
!
!  Compute the potential energy and forces.
!
    f(1:nd,i) = 0.0
    do j = 1, np
         if ( i /= j ) then
           call dist ( nd, pos(1,i), pos(1,j), rij, d )
!
!  Attribute half of the potential energy to particle J.
!
           pot = pot + 0.5 * v(d)

           do k = 1, nd
             f(k,i) = f(k,i) - rij(k) * dv(d) / d
           end do

         end if
    end do
!
!  Compute the kinetic energy.
!
    kin = kin + dotr8 ( nd, vel(1,i), vel(1,i) )

  end do
!$omp  end parallel do

  kin = kin * 0.5 * mass
  
  return
end
subroutine initialize ( np, nd, box, pos, vel, acc )
!
!*******************************************************************************
!
!! INITIALIZE initializes the positions, velocities, and accelerations.
!
!
!  Modified:
!
!    18 March 2002
!
!  Author:
!
!    Bill Magro, 
!    Kuck and Associates, Inc. (KAI), 1998
!
!  Parameters:
!
!    Input, integer NP, the number of particles.
!
!    Input, integer ND, the number of spatial dimensions.
!
!    Input, double precision BOX(ND), specifies the maximum position
!    of particles in each dimension.
!
!    Output, double precision POS(ND,NP), the position of each particle.
!
  implicit none
!
  integer np
  integer nd
!
  double precision acc(nd,np)
  double precision box(nd)
  integer i
  double precision pos(nd,np)
  double precision vel(nd,np)
  double precision x
!
!  Give the particles random positions within the box.
!
  call random_number ( harvest = pos(1:nd,1:np) )

  do i = 1, nd
    pos(i,1:np) = pos(i,1:np) * box(i)
  end do

  vel(1:nd,1:np) = 0.0D+00
  acc(1:nd,1:np) = 0.0D+00

  return
end
subroutine dist ( nd, r1, r2, dr, d )
!
!*******************************************************************************
!
!! DIST computes the displacement (and its norm) between two particles.
!
!
!  Modified:
!
!    17 March 2002
!
!  Author:
!
!    Bill Magro, 
!    Kuck and Associates, Inc. (KAI), 1998
!
!  Parameters:
!
!    Input, integer ND, the number of spatial dimensions.
!
  implicit none
!
  integer nd
!
  double precision box(nd)
  double precision d
  double precision dr(nd)
  double precision r1(nd)
  double precision r2(nd)
!
  dr(1:nd) = r1(1:nd) - r2(1:nd)

  d = sqrt ( sum ( dr(1:nd)**2 ) )

  return
end
function dotr8 ( n, x, y )
!
!*******************************************************************************
!
!! DOTR8 returns the dot product of two double precision vectors.
!
!
!  Modified:
!
!    17 March 2002
!
!  Author:
!
!    Bill Magro, 
!    Kuck and Associates, Inc. (KAI), 1998
!
!  Parameters:
!
  implicit none
!
  integer n
!
  double precision dotr8
  integer i
  double precision x(n)
  double precision y(n)
!
  dotr8 = 0.0
  do i = 1, n
    dotr8 = dotr8 + x(i)*y(i)
  end do

  return
end
subroutine update ( np, nd, pos, vel, f, a, mass, dt )
!
!*******************************************************************************
!
!! UPDATE performs the time integration, using a velocity Verlet algorithm.
!
!
!  Discussion:
!
!    The time integration is fully parallel.
!
!  Modified:
!
!    17 March 2002
!
!  Author:
!
!    Bill Magro, 
!    Kuck and Associates, Inc. (KAI), 1998
!
!  Parameters:
!
!    Input, integer NP, the number of particles.
!
!    Input, integer ND, the number of spatial dimensions.
!
!    Input/output, double precision POS(ND,NP), the position of each particle.
!
  implicit none

  integer np
  integer nd
!
  double precision a(nd,np)
  double precision dt
  double precision f(nd,np)
  integer i
  integer j
  double precision mass
  double precision pos(nd,np)
  double precision rmass
  double precision vel(nd,np)
!
  rmass = 1.0 / mass
!
!$omp  parallel do
!$omp& default(shared)
!$omp& private(i,j)
  do j = 1, np
    do i = 1, nd
      pos(i,j) = pos(i,j) + vel(i,j) * dt + 0.5 * a(i,j) * dt**2
      vel(i,j) = vel(i,j) + 0.5 * dt * ( f(i,j) * rmass + a(i,j))
      a(i,j) = f(i,j) * rmass
    end do
  end do
!$omp  end parallel do

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end