**************************************************************************** C LONG BUI, DARYL EBANKS, JONATHON WILLIAMS C HIGH PERFORMANCE COMPUTING - CSC 435 C DR POUNDS C APRIL 2007 C C MOLECULAR DYNAMICS C C THIS PROGRAM DISPLAYS THE X,Y,Z COORDINATES FOR A TRAJECTORY OF 5 ATOMS C BEING MOVED BY THE FORCES AMONGST 2000 ATOMS IN A BOX ***************************************************************************** program proj4 include '/usr/local/include/mpif.h' parameter (MASTER = 0) parameter (FROM_MASTER = 1) parameter (FROM_WORKER = 2) parameter (a=20) !a==atom_count parameter (total=2d0)!timetotal parameter (h=0.01d0)!timestep parameter (dim=100d0)!size of box integer status(4), dest, source, offsetavg, term, atoms integer numworkers, offset, taskid, numtasks, ierr, mtype double precision x(a),y(a),z(a),vx(a),vy(a),vz(a), TEMP(3), & fx(a,a),fy(a,a),fz(a,a),fSUMx(a),fSUMy(a),fSUMz(a), & prevFxSUM(a),prevFySUM(a),prevFzSUM(a) double precision atomX(total,5),atomY(total,5),atomZ(total,5) DOUBLE PRECISION FORCECALC, DISTANCE EXTERNAL FORCECALC, DISTANCE integer time, n, ISEED(4) n=a C INITIALIZE SEED VALUES FOR RANDOM GENERATOR ISEED(1) = 133 ISEED(2) = 31 ISEED(3) = 1333 ISEED(4) = 131 C INITIALIZE MPI CODE SECTION AND DECLARE RANK OF PROCESSOR AND NUMBER OF PROCESSORS call MPI_INIT( ierr ) call MPI_COMM_RANK( MPI_COMM_WORLD, taskid, ierr ) call MPI_COMM_SIZE( MPI_COMM_WORLD, numtasks, ierr ) numworkers = numtasks-1 offsetavg = n / numworkers offset = offsetavg+1 extra = mod(n,numworkers) C *************************** master task ************************************* if (taskid .eq. MASTER) then C INITIALIZE POSITIONS UNIFORMLY RANDOM WITHIN THE CORRECT RANGE (0,DIM) do i=1,n call dlaruv(ISEED,3,TEMP) x(i) = TEMP(1) * DIM y(i) = TEMP(2) * DIM z(i) = TEMP(3) * DIM enddo do i=1,5 print 1001, x(i)," ",y(i)," ",z(i) enddo do time = 1,total/h C CALCULATE FORCES! !longs calcs (small number of iterations, better done serially) do j=1,n do k=j+1,n C This calculates all the X forces FX(J,K) = FORCECALC(1,J,K,X,Y,Z,N) FX(K,J) = 1.0D0 * FX(J,K) C This calculates all the Y forces FY(J,K) = FORCECALC(2,J,K,X,Y,Z,N) FY(K,J) = 1.0D0 * FY(J,K) C This calculates all the Z forces FZ(J,K) = FORCECALC(3,J,K,X,Y,Z,N) FZ(K,J) = 1.0D0 * FZ(J,K) enddo enddo! END CALCULATING FORCES ! CALCULATE VELOCITY (AND NEW FSUM) ! daryls calcs !get an array position and size of data so that it distributes amongst the processors aveatoms = N/numworkers extra = mod(N, numworkers) offset = 1 mtype = FROM_MASTER do dest=1, numworkers atoms = aveatoms if (dest .le. extra) then atoms = atoms + 1 endif ! send data to the workers call MPI_SEND( offset, 1, MPI_INTEGER, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( atoms, 1, MPI_INTEGER, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( FX, N*N, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( FY, N*N, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( FZ, N*N, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr ) offset = offset + atoms enddo!endsend !receive data calculated from workers mtype = FROM_WORKER do source = 1,numworkers call MPI_RECV( offset, 1, MPI_INTEGER, source, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV( atoms, 1, MPI_INTEGER, source, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV(vx(offset), atoms, MPI_double_precision, source, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV(vy(offset), atoms, MPI_double_precision, source, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV(vz(offset), atoms, MPI_double_precision, source, & mtype, MPI_COMM_WORLD, status, ierr ) enddo!endrecv ! CALCULATE NEW POSITION !jons calcs !get an array position and size of data so that it distributes amongst the processors term=1 do dest=1,numworkers if (dest.le.extra) then offset = offsetavg+1 else offset = offsetavg endif C Send data to the worker tasks mtype = FROM_MASTER call MPI_SEND( term, 1, MPI_INTEGER, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( offset, 1, MPI_INTEGER, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( x(term), offset, MPI_DOUBLE_PRECISION, & dest, mtype, MPI_COMM_WORLD, ierr ) call MPI_SEND( y(term), offset, MPI_DOUBLE_PRECISION, & dest, mtype, MPI_COMM_WORLD, ierr ) call MPI_SEND( z(term), offset, MPI_DOUBLE_PRECISION, & dest, mtype, MPI_COMM_WORLD, ierr ) call MPI_SEND( FSUMX, n, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( FSUMY, n, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( FSUMZ, n, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( vx, n, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( vy, n, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( vz, n, MPI_DOUBLE_PRECISION, dest, mtype, & MPI_COMM_WORLD, ierr )!i need a lot of data sent term = term + offset enddo!endsend C Receive data from the worker tasks mtype = FROM_WORKER do source=1,numworkers call MPI_RECV( term, 1, MPI_INTEGER, source, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( offset, 1, MPI_INTEGER, source, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( x(term), offset, MPI_DOUBLE_PRECISION, & source, mtype, MPI_COMM_WORLD, status, ierr ) !recv the portion of x newly calculated call MPI_RECV( y(term), offset, MPI_DOUBLE_PRECISION, & source, mtype, MPI_COMM_WORLD, status, ierr ) !recv the portion of y newly calculated call MPI_RECV( z(term), offset, MPI_DOUBLE_PRECISION, & source, mtype, MPI_COMM_WORLD, status, ierr ) !recv the portion of z newly calculated enddo!endrecv ! END OUTTER POSITION SUM do i = 1,5 atomx(time,i) = x(i) !KEEP TRACK OF 5 PARTICLES atomy(time,i) = y(i) atomz(time,i) = z(i) print 1001, atomX(time,i)," ",atomY(time,i)," ",atomZ(time,i) enddo !end store enddo!endtimestep C do time = 1,total/h C do i = 1,5 C print 1001, atomX(time,i)," ",atomY(time,i)," ",atomZ(time,i) C enddo C enddo endif C *************************** worker task ************************************* if (taskid > MASTER) then !velocity calculations worker (daryls calcs continued) do time = 1,total/h,numworkers mtype = FROM_MASTER call MPI_RECV( offset, 1, MPI_INTEGER, MASTER, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV( atoms, 1, MPI_INTEGER, MASTER, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV( FX, N*N, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV( FY, N*N, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV( FZ, N*N, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, status, ierr ) do i=1,N VX(i) = 0D0 VY(i) = 0D0 VZ(i) = 0D0 enddo c calculate velocity matrices CALL CALCVELOCITY(OFFSET, ATOMS, N, Fx, prevfxsum, Vx) CALL CALCVELOCITY(OFFSET, ATOMS, N, Fy, prevfysum, Vy) CALL CALCVELOCITY(OFFSET, ATOMS, N, Fz, prevfzsum, Vz) do i=OFFSET, OFFSET + ATOMS-1 Vx(i) = h*Vx(i) Vy(i) = h*Vy(i) Vz(i) = h*Vz(i) !print *, i, Vx(i), Vy(i), Vz(i) enddo mtype = FROM_WORKER call MPI_SEND( offset, 1, MPI_INTEGER, MASTER, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( atoms, 1, MPI_INTEGER, MASTER, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( VX(offset), atoms, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, ierr ) call MPI_SEND( VY(offset), atoms, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, ierr ) call MPI_SEND( VZ(offset), atoms, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, ierr ) !position calculations worker (jons work continued) C Receive matrix data from master task mtype = FROM_MASTER call MPI_RECV( term, 1, MPI_INTEGER, MASTER, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( offset, 1, MPI_INTEGER, MASTER, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( x(term), offset, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV( y(term), offset, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV( z(term), offset, MPI_DOUBLE_PRECISION, MASTER, & mtype, MPI_COMM_WORLD, status, ierr ) call MPI_RECV( fSUMx, n, MPI_DOUBLE_PRECISION, MASTER, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( fSUMy, n, MPI_DOUBLE_PRECISION, MASTER, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( fSUMz, n, MPI_DOUBLE_PRECISION, MASTER, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( vx, n, MPI_DOUBLE_PRECISION, MASTER, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( vy, n, MPI_DOUBLE_PRECISION, MASTER, mtype, & MPI_COMM_WORLD, status, ierr ) call MPI_RECV( vz, n, MPI_DOUBLE_PRECISION, MASTER, mtype, & MPI_COMM_WORLD, status, ierr ) C OUTER POSITION SUM (CONT.) do worker's portion do i = term,offset C INNER SUM done previously as fsum(1:n) x(i) = x(i) + h*vx(i) + h**2*fSUMx(i) y(i) = y(i) + h*vy(i) + h**2*fSUMy(i) z(i) = z(i) + h*vz(i) + h**2*fSUMz(i) C BOUNDS CHECK x(i) = mod(x(i),dim) x(i) = mod(x(i) + dim,dim) y(i) = mod(y(i),dim) y(i) = mod(y(i) + dim,dim) z(i) = mod(z(i),dim) z(i) = mod(z(i) + dim,dim) enddo !END OUTER POSITION SUM (CONT.) mtype = FROM_WORKER !send data back to master call MPI_SEND( term, 1, MPI_INTEGER, MASTER, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( offset, 1, MPI_INTEGER, MASTER, mtype, & MPI_COMM_WORLD, ierr ) call MPI_SEND( x(term), offset, MPI_DOUBLE_PRECISION, & MASTER, mtype, MPI_COMM_WORLD, status, ierr ) call MPI_SEND( y(term), offset, MPI_DOUBLE_PRECISION, & MASTER, mtype, MPI_COMM_WORLD, status, ierr ) call MPI_SEND( z(term), offset, MPI_DOUBLE_PRECISION, & MASTER, mtype, MPI_COMM_WORLD, status, ierr ) enddo!end timestep within worker endif !end worker task 1001 format (3(f20.17,a1)) call MPI_FINALIZE(ierr) !end MPI CODE SECTION end !end program DOUBLE PRECISION FUNCTION FORCECALC(D, I, J, X, Y, Z, N) INTEGER N, I, J, D DOUBLE PRECISION X(N), Y(N), Z(N) DOUBLE PRECISION DIST DOUBLE PRECISION DISTANCE EXTERNAL DISTANCE DIST = DISTANCE(X(I), X(J), Y(I), Y(J), Z(I), Z(J)) ! If we're calculating X component of force, then D = 1 IF (D .EQ. 1) THEN FORCECALC = (X(I) - X(J)) * (DIST ** -14 - DIST ** -8) ! If we're calculating Y component of force, then D = 2 ELSE IF (D .EQ. 2) THEN FORCECALC = (Y(I) - Y(J)) * (DIST ** -14 - DIST ** -8) ! If we're calculating Z component of force, then D != 1 && D != 2 ELSE FORCECALC = (Z(I) - Z(J)) * (DIST ** -14 - DIST ** -8) ENDIF RETURN END ! end forcecalc subroutine DOUBLE PRECISION FUNCTION DISTANCE(X1, X2, Y1, Y2, Z1, Z2) DOUBLE PRECISION X1, X2, Y1, Y2, Z1, Z2 DISTANCE = DSQRT((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1)+(Z2-Z1)*(Z2-Z1)) RETURN END !end distance subroutine SUBROUTINE CALCVELOCITY(OFFSET, ATOMS, & N, F, prevfsum, V) INTEGER I, J, N, OFFSET, ATOMS DOUBLE PRECISION F(N,N), prevfsum(N), fsum(N), V(N) do j=OFFSET,OFFSET+ATOMS-1 !sum forces !to save time i save the prev F-sums fsum(j) = 0 do i = 1,N !sum of the forces on particle j fsum(j) = fsum(j) + F(i,j) enddo V(j) =+ .05D0*(fsum(j) + prevfsum(j)) !prev sum is stored prevfsum(j) = fsum(j) enddo RETURN END !end calcvelocity subroutine