Module Random !======================================================================= ! Portable pseudo-random integer generator, especially for ! microcomputers with a limitation of 16 bit integers. Translated from ! original Pascal version(1) to Fortran 77 by H. D. Knoble, PSU. ! ! Corrected typo (change 30308.0 to 30307.0) courtesy of Ed Wynn ! and Dave Dodson, 3 May 2002. ! ! The supporting paper is: ! (1) B. Wichmann & D. Hill, "Building a Random-Number Generator", ! BYTE, March, 1987, 12(3):127-128. ! ! Also see the following related works: ! (2) Wichmann, B.A. and Hill, I.D., "An Efficient and Portable", ! Pseudo-random Number Generator", Applied Statistics, ! Volume 31, 1982, pages 188-190. ! (3) Haas, Alexander, "The Multiple Prime Random Number Generator", ! ACM Transactions on Mathematical Software; December, ! 1987, 13(4):368-381. ! (4) L'Ecuyer, Pierre, "Efficient and Portable Combined Random Number ! Generators", Communications of the ACM; June, 1988, ! 31(6):742-749,774. ! ! Use... ! CALL RAND(U,N) ! To generate a sequence, U, of N Uniform(0,1) numbers. ! Cycle length is ((30269-1)*(30307-1)*(30323-1))/4 or ! 6953607871644 > 6.95E+12. ! ! To access the SEED vector in the calling program use statements: ! use Ramdom ! ! The common variable SEED is the array of three current seeds. INTEGER :: X,Y,Z contains !======================================================================= SUBROUTINE RAND(U,N) REAL :: U(:),V IF(N <= 0) RETURN DO I=1,N X=171*MOD(X,177)-2*(X/177) IF(X.LT.0) X=X+30269 Y=172*MOD(Y,176)-35*(Y/176) IF(Y.LT.0) Y=Y+30307 Z=170*MOD(Z,178)-63*(Z/178) IF(Z.LT.0) Z=Z+30323 V=X/30269.0 + Y/30307.0 + Z/30323.0 U(I)=V-INT(V) END DO RETURN end subroutine RAND end module random program md !===================================================================== ! This program implements a simple molecular dynamics simulation, ! using the velocity Verlet time integration scheme. The particles ! interact with a central pair potential. ! ! Author: Bill Magro, Kuck and Associates, Inc. (KAI), 1998 ! ! Parallelism is implemented via OpenMP directives. ! THIS PROGRAM USES THE FORTRAN90 RANDOM_NUMBER FUNCTION AND ARRAY ! SYNTAX !===================================================================== ! Removed REAL*8 and added real time clock timings. ! Added own random number generator to facilitate identical results ! accross compilers and platforms. ! Original program may be found at: ! http://www.openmp.org/drupal/samples/md.html ! hdk !===================================================================== implicit none ! simulation parameters: ! ndim = dimensionality of the physical space ! nparts = number of particles ! nsteps = number of time steps in the simulation ! mass = mass of the particles ! dt = time step integer, parameter :: ndim=3, nparts=500, nsteps=1000 integer, parameter :: DP = selected_real_kind(15) real(kind=DP) :: mass=1.0_DP real(kind=DP) :: dt=1.0e-4_DP ! simulation variables real(kind=DP) :: box(ndim) ! dimensions of the simulation box real(kind=DP) :: position(ndim,nparts), velocity(ndim,nparts) real(kind=DP) :: force(ndim,nparts), accel(ndim,nparts) real(kind=DP) :: potential, kinetic, E0 real :: T1,T2, Seconds integer :: i, pulses, PPS box(1:ndim) = 10.0 ! set initial positions, velocities, and accelerations call initialize(nparts,ndim,box,position,velocity,accel) CALL SYSTEM_CLOCK(COUNT=Pulses,COUNT_RATE=PPS) T1 = REAL(Pulses,DP)/PPS ! compute the forces and energies call compute(nparts,ndim,position,velocity,mass, & force,potential,kinetic) E0 = potential + kinetic ! This is the main time stepping loop do i=1,nsteps call compute(nparts,ndim,position,velocity,mass, & force,potential,kinetic) print *, potential, kinetic,(potential + kinetic - E0)/E0 call update(nparts,ndim,position,velocity,force,accel,mass,dt) enddo CALL SYSTEM_CLOCK(COUNT=Pulses,COUNT_RATE=PPS) T2 = REAL(Pulses,DP)/PPS Seconds=T2-T1 print *," Number of particles: ", nparts print *," Number of time steps:", nsteps print *," Time in seconds for the simulation: ",Seconds end program md subroutine compute(np,nd,pos,vel,mass,f,pot,kin) !--------------------------------------------------------------------- ! Compute the forces and energies, given positions, masses, ! and velocities !--------------------------------------------------------------------- implicit none integer :: np, nd, i, j, k integer, parameter :: DP = selected_real_kind(15) real(kind=DP) :: pos(nd,np), vel(nd,np), f(nd,np) real(kind=DP) :: mass, pot, kin, dotr8 real(kind=DP) :: rij(nd), d, v, dv real(kind=DP) :: PI2=3.141592653589793238462643_DP/2.0_DP pot = 0.0 kin = 0.0 ! The computation of forces and energies is fully parallel. !$OMP PARALLEL DO & !$OMP DEFAULT(SHARED) & !$OMP PRIVATE(I,J,K,RIJ,D) & !$OMP REDUCTION(+ : POT, KIN) do i=1,np ! compute 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_DP*v(d) do k=1,nd f(k,i) = f(k,i) - rij(k)*dv(d)/d enddo endif enddo ! compute kinetic energy kin = kin + dotr8(nd,vel(1,i),vel(1,i)) enddo !$OMP END PARALLEL DO kin = kin*0.5_DP*mass end subroutine compute subroutine initialize(np,nd,box,pos,vel,acc) !--------------------------------------------------------------------- ! Initialize the positions, velocities, and accelerations. ! The Fortran90 random_number function is used to choose positions. !--------------------------------------------------------------------- use random implicit none integer :: np, nd, i, j integer, parameter :: DP = selected_real_kind(15) real(kind=DP) :: box(nd), pos(nd,np), vel(nd,np), acc(nd,np) real :: xx(1) ! Initialize Random Number Generator X = 1 Y = 10000 Z = 3000 do i=1,np do j=1,nd call rand(xx,1) pos(j,i) = box(j)*xx(1) vel(j,i) = 0.0 acc(j,i) = 0.0 enddo enddo end subroutine initialize ! Compute the displacement vector (and its norm) ! between two particles. subroutine dist(nd,r1,r2,dr,d) implicit none integer :: nd, i integer, parameter :: DP = selected_real_kind(15) real(kind=DP) :: r1(nd), r2(nd), dr(nd), d d = 0.0 do i=1,nd dr(i) = r1(i) - r2(i) d = d + dr(i)**2 enddo d = sqrt(d) end subroutine dist ! Return the dot product between two vectors of type ! Double Precision and length n function dotr8(n,x,y) result(DotResult) implicit none integer, parameter :: DP = selected_real_kind(15) integer :: n real(kind=DP) :: x(n), y(n), DotResult DotResult = Sum(x*y) end function dotr8 subroutine update(np,nd,pos,vel,f,a,mass,dt) !--------------------------------------------------------------------- ! Perform the time integration, using a velocity Verlet algorithm !--------------------------------------------------------------------- implicit none integer :: np, nd, i, j integer, parameter :: DP = selected_real_kind(15) real(kind=DP) :: pos(nd,np), vel(nd,np), f(nd,np), a(nd,np) real(kind=DP) :: mass, dt, rmass rmass = 1.0_DP/mass ! The time integration is fully parallel !$OMP PARALLEL DO & !$OMP DEFAULT(SHARED) & !$OMP PRIVATE(I,J) do i = 1,np do j = 1,nd pos(j,i) = pos(j,i) + vel(j,i)*dt + 0.5_DP*dt*dt*a(j,i) vel(j,i) = vel(j,i) + 0.5_DP*dt*(f(j,i)*rmass + a(j,i)) a(j,i) = f(j,i)*rmass enddo enddo !$OMP END PARALLEL DO return end function v(x) result(vr) ! function v for the pair potential and its derivative, dv. ! This potential is a harmonic well which smoothly saturates to a ! maximum value at PI/2. integer, parameter :: DP = selected_real_kind(15) real(kind=DP) vr, x real(kind=DP) :: PI2=3.141592653589793238462643_DP/2.0_DP vr = sin(min(x,PI2))**2 end function v function dv(x) result(dvr) integer, parameter :: DP = selected_real_kind(15) real(kind=DP) dvr, x real(kind=DP) :: PI2=3.141592653589793238462643_DP/2.0_DP dvr = 2.0_DP*sin(min(x,PI2))*cos(min(x,PI2)) end function dv