        program md
	implicit none
        include 'mpif.h'
        integer NCELLS,N
        parameter(NCELLS=50)
        parameter(N=NCELLS*NCELLS*NCELLS)
        double precision x(N),y(N),z(N)
        double precision fx(N),fy(N),fz(N)
        double precision epot
        double precision dx,dy,dz,dr2,dr
        double precision myrandom
        double precision sigma, epsilon, BoxL,BoxL2
        double precision tzero, tnow, mytimer
        double precision atdensity
        double precision cutoff
        double precision mw, density, molvolume, boxvolume
        double precision sumf2
        integer i,j
c. additional vars : e4=4*e, e24=24*e
        double precision e4,e24
c.vars to keep sigma/r^6 and sigma/r ^12
c. cutsq square of cutoff
        double precision sigr6,sigr12, cutsq, sigr126
        double precision e24sigr126,dxdr,dydr,dzdr
c
        integer me,nprocs,ierror
        integer Status(MPI_STATUS_SIZE)
        integer msgtag,src
        integer k,len,ierr
        character*(MPI_MAX_PROCESSOR_NAME) hn
        double precision tmp1, tmpA(N)
        integer ii

        me = -1

c       
        CALL MPI_Init(ierror)

        if ( ierror.ne.MPI_SUCCESS ) then
            write(*,*)'MPI_Init() failed. Error_Code : ',ierror
            stop
        endif

        CALL MPI_Comm_Rank(MPI_COMM_WORLD,me,ierror)

        if ( ierror .ne. MPI_SUCCESS) Then
            write(*,*)'MPI_Comm_Rank() failed. Error_Code : ',ierror
            CALL MPI_Finalize(ierror)
        endif
            
        CALL MPI_Comm_Size(MPI_COMM_WORLD,nprocs,ierror)
        if ( ierror .ne. MPI_SUCCESS) Then
            write(*,*)'MPI_Comm_Size() failed. Error_Code : ',ierror
            CALL MPI_Finalize(ierror)
        endif
    
        call MPI_GET_PROCESSOR_NAME(hn, len, ierr) 

        write(*,'(A4,1X,I4,1X,A4,1X,I4,1X,A4,1X,A10)')
     $ 'proc ',me,' of ',nprocs,' on ',hn


        
        mw = 131.290d0
	density = 1.779800d0
        molvolume =  mw / density
        molvolume = molvolume * dble(N)*1.0d24 / 6.0225d23
c molvolume in Angstrom^3       
        boxL = molvolume**(1.0d0/3.0d0)
        epsilon = 221.0d0
        sigma   = 3.9d0
        cutoff = 6.0 * sigma
        BoxL2=boxL/2.0d0
        cutsq = cutoff * cutoff
       

        e4 = 4.0d0 * epsilon
        e24 = 24.0d0 * epsilon

c Initialize Cubic Cell coordinates, zero energy, zero forces
        epot=0.0d0
        if (me.eq.0) then
        call makeInitialConf(N,boxL,x,y,z)

        do i = 1, N
            fx(i)=0.0000000d0
            fy(i)=0.0000000d0
            fz(i)=0.0000000d0
        enddo

       endif

c        Broadcast theese data to all slaves

       call MPI_Bcast(x,N,MPI_DOUBLE_PRECISION,0,
     $                MPI_COMM_WORLD,ierror)
       if ( ierror.ne.MPI_SUCCESS ) then
         write(*,*)'MPI_Bcast() failed. Error_Code : ',ierror
         CALL MPI_Finalize(ierror)
       endif

       call MPI_Bcast(y,N,MPI_DOUBLE_PRECISION,0,
     $                MPI_COMM_WORLD,ierror)
       if ( ierror.ne.MPI_SUCCESS ) then
         write(*,*)'MPI_Bcast() failed. Error_Code : ',ierror
         CALL MPI_Finalize(ierror)
       endif

       call MPI_Bcast(z,N,MPI_DOUBLE_PRECISION,0,
     $                MPI_COMM_WORLD,ierror)
       if ( ierror.ne.MPI_SUCCESS ) then
         write(*,*)'MPI_Bcast() failed. Error_Code : ',ierror
         CALL MPI_Finalize(ierror)
       endif


       call MPI_Bcast(fx,N,MPI_DOUBLE_PRECISION,0,
     $                MPI_COMM_WORLD,ierror)
       if ( ierror.ne.MPI_SUCCESS ) then
         write(*,*)'MPI_Bcast() failed. Error_Code : ',ierror
         CALL MPI_Finalize(ierror)
       endif

       call MPI_Bcast(fy,N,MPI_DOUBLE_PRECISION,0,
     $                MPI_COMM_WORLD,ierror)
       if ( ierror.ne.MPI_SUCCESS ) then
         write(*,*)'MPI_Bcast() failed. Error_Code : ',ierror
         CALL MPI_Finalize(ierror)
       endif

       call MPI_Bcast(fz,N,MPI_DOUBLE_PRECISION,0,
     $                MPI_COMM_WORLD,ierror)
       if ( ierror.ne.MPI_SUCCESS ) then
         write(*,*)'MPI_Bcast() failed. Error_Code : ',ierror
         CALL MPI_Finalize(ierror)
       endif

c  Timer
	tzero = MPI_WTime()

        do ii = 1 + me, N - 1, nprocs
          if ( mod(ii,2).eq.0 ) then
                i=N -1 -ii
          else
               i = ii
         endif
            do j = i + 1, N
               dx = x(i) - x(j)
               dy = y(i) - y(j)
               dz = z(i) - z(j)
c  Boundary conditions
               if ( dx .lt. -BoxL2 ) then
                    dx = dx + BoxL
               endif
               if ( dy .lt. -BoxL2 ) then
                    dy = dy + BoxL
               endif
               if ( dz .lt. -BoxL2 ) then
                    dz = dz + BoxL
               endif

               if ( dx .gt. BoxL2 ) then
                    dx = dx - BoxL
               endif
               if ( dy .gt. BoxL2 ) then
                    dy = dy - BoxL
               endif
               if ( dz .gt. BoxL2 ) then
                    dz = dz - BoxL
               endif

               dr2 = dx*dx + dy*dy + dz*dz

c... Cutoff
               if ( dr2.le.cutsq ) then
                 dr  = dsqrt(dr2)
                 sigr6 = (sigma/dr)**6
                 sigr12 = sigr6 * sigr6
                 sigr126=2.0*sigr12 - sigr6
                 dxdr = dx/dr
                 dydr = dy/dr
                 dzdr = dz/dr
                 epot = epot + e4* (sigr12 - sigr6 )
                 e24sigr126 = e24 * sigr126
                 fx(i) = fx(i) - e24sigr126*dxdr
                 fx(j) = fx(j) + e24sigr126*dxdr
                 fy(i) = fy(i) - e24sigr126*dydr
                 fy(j) = fy(j) + e24sigr126*dydr
                 fz(i) = fz(i) - e24sigr126*dzdr
                 fz(j) = fz(j) + e24sigr126*dzdr
               endif
            enddo
        enddo

cccc Now collect the results : epot and fx, fy, fz
        
        if (me.eq.0 ) then
           msgtag = MPI_ANY_TAG 
           do i = 1, nprocs - 1
             CALL MPI_Recv(tmp1,1,MPI_DOUBLE_PRECISION,i,msgtag,
     %    MPI_COMM_WORLD,Status,ierror)
          epot =epot + tmp1

             CALL MPI_Recv(tmpA,N,MPI_DOUBLE_PRECISION,i,msgtag,
     %    MPI_COMM_WORLD,Status,ierror)
              do j = 1, N
               fx(i) = fx(i) + tmpA(i)
              enddo

             CALL MPI_Recv(tmpA,N,MPI_DOUBLE_PRECISION,i,msgtag,
     %    MPI_COMM_WORLD,Status,ierror)
              do j = 1, N
               fy(i) = fy(i) + tmpA(i)
              enddo

             CALL MPI_Recv(tmpA,N,MPI_DOUBLE_PRECISION,i,msgtag,
     %    MPI_COMM_WORLD,Status,ierror)
              do j = 1, N
               fz(i) = fz(i) + tmpA(i)
              enddo

          enddo
        else
          msgtag = 1234
              CALL MPI_Send(epot,1,MPI_DOUBLE_PRECISION,0,msgtag,
     %    MPI_COMM_WORLD,ierror)
              CALL MPI_Send(fx,N,MPI_DOUBLE_PRECISION,0,msgtag,
     %    MPI_COMM_WORLD,ierror)
              CALL MPI_Send(fy,N,MPI_DOUBLE_PRECISION,0,msgtag,
     %    MPI_COMM_WORLD,ierror)
              CALL MPI_Send(fz,N,MPI_DOUBLE_PRECISION,0,msgtag,
     %    MPI_COMM_WORLD,ierror)
        endif

        if ( me .eq. 0 ) then
        tnow = MPI_WTime() - tzero


        epot = epot / dble(N)
        write(*,1)'Epot     : ',epot
        write(*,2)'WallTime : ',tnow
        endif
        
        CALL MPI_Finalize(ierror)

1       format(A10,1X,F15.4)
2       format(A10,1X,F15.4,' sec')
        end

 
c............................................
	Subroutine makeInitialConf(N,boxL,x,y,z)
	implicit none
        integer N
	double precision x(N),y(N),z(N),boxL,un
        integer ndivs,tmp_index
        integer i,j,k,jj
        double precision xa,ya,za
	
        ndivs=int( (N)**(1.0/3.0) + 0.01)
        tmp_index = ndivs*ndivs*ndivs
        if ( tmp_index .ne. N) then
            write(*,*)'N does not corrspond to Cubic Cell'
            stop
        endif
        un = 1.0 / dble(ndivs)

         jj = 0
         xa = -0.5 - 0.5*un
         do i = 1, ndivs
              xa = xa + un
              ya = -0.5 - 0.5*un
           do j = 1,ndivs
              ya = ya + un
              za = -0.5 - 0.5*un
              do k= 1,ndivs
                za = za + un
                jj = jj + 1
                x(jj) = xa
                y(jj) = ya
                z(jj) = za
              enddo
           enddo
         enddo


        do i = 1, N
            x(i) = x(i) * boxL
            y(i) = y(i) * boxL
            z(i) = z(i) * boxL
        enddo

        return
	end

