c    ===================================================================
        program md
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Molecular dynamics of Lennard-Jones atoms                      |
c    |  Atom number 1 has a hard-wired mass, different from the others |
c    ===================================================================

c    ** Subprograms referenced **
        external start, runner, finish

c    -------------------------------------------------------------------

        call start
	call runner
	call finish
	
	stop
	end
	
	
	
c    ===================================================================
        subroutine start
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Reads in run parameters and initial configuration              |
c    |  Calculates and prints out initial values                       |
c    ===================================================================

C    ** Local variables **
        integer         i, nn
	real            pot, vir, kin

c    ** Subprograms referenced **
        external getintg, getreal, force, kinetic
	
c    -------------------------------------------------------------------
	
c    ** Output program identification **	
        print *, '-----------------------------------------------------'	
	print *, 'Molecular dynamics simulation of Lennard-Jones atoms'
        print *, 'Atom number 1 has a different mass'
	print *, 'Cubic periodic boundary conditions'
	print *, 'Velocity Verlet algorithm'

c    ** Input run parameters **	
        print *, '-----------------------------------------------------'	
	call getintg ( 'number of blocks     ', nblock )
	call getintg ( 'number of steps/block', nstep  )
	call getreal ( 'time step            ', dt     )
	call getreal ( 'potential cutoff     ', rcut   )

c    ** Input starting configuration **	
        print *, '-----------------------------------------------------'	
        print *, 'Reading initial configuration from file md.old'
	
	open(unit=3,file='md.old',status='old')
	read(3,*) nn
	if ( nn .ne. n ) stop 'Error in number of particles'
	read(3,*) box
        do i = 1, n
	   read(3,*) rx(i), ry(i), rz(i), vx(i), vy(i), vz(i)
	enddo
	close(unit=3)

c    ** Output initial values **
	call force ( pot, vir )
	call kinetic ( kin )
	eng   = ( kin + pot ) / real(n)
	temp  = 2.0*kin/(3.0*real(n))
        dens  = real(n)/(box**3)
	press = dens*temp+vir/(box**3)
        print *, '-----------------------------------------------------'	
	print *, 'Initial values'
	write(*,'(1x,''Energy/atom = '',f10.4)') eng
	write(*,'(1x,''Density     = '',f10.4)') dens
	write(*,'(1x,''Pressure    = '',f10.4)') press
	write(*,'(1x,''Temperature = '',f10.4)') temp
		
	return
	end


	
c    ===================================================================
        subroutine finish
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Calculates and prints out final values                         |
c    |  Writes out final configuration                                 |
c    ===================================================================

C    ** Local variables **
	integer i
	real            pot, vir, kin

c    ** Subprograms referenced **
        external force, kinetic
	
c    -------------------------------------------------------------------

c    ** Output final values **
	call force ( pot, vir )
	call kinetic ( kin )
	eng   = ( kin + pot ) / real(n)
        temp  = 2.0*kin/( 3.0 * real(n) )
        dens  = real(n)/(box**3)
	press = dens*temp+vir/(box**3)
        print *, '-----------------------------------------------------'	
	print *, 'Final values'
	write(*,'(1x,''Energy/atom = '',f10.4)') eng
	write(*,'(1x,''Density     = '',f10.4)') dens
	write(*,'(1x,''Pressure    = '',f10.4)') press
	write(*,'(1x,''Temperature = '',f10.4)') temp

c    ** Output final configuration **
        print *, '-----------------------------------------------------'	
        print *, 'Writing final configuration to file md.new'	
	open(unit=3,file='md.new',status='unknown')
	write(3,'(1x,i5)') n
	write(3,'(1x,3f10.5)') box
        do i = 1, n
	   write(3,'(1x,6f10.5)') rx(i), ry(i), rz(i),
     :                            vx(i), vy(i), vz(i)
	enddo
	close(unit=3)
	
        print *, '-----------------------------------------------------'	
	print *, 'Program ends'
        print *, '-----------------------------------------------------'	

        return
	end
	
	
c    ===================================================================
        subroutine runner
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Conducts molecular dynamics run, nblock blocks of nstep steps  |
c    ===================================================================

C    ** Local variables **
        integer         block, step
	real            pot, vir, kin

c    ** Subprograms referenced **
	external force, movea, moveb
	external blkzer, blkacc, blkout
	external runzer, runacc, runout
	
c    -------------------------------------------------------------------
	
        write(*,'(/1x,60(''=''))')
        write(*,'(1x,''          Avg.Period'',
     :               ''    Energy   Density  Pressure      Temp'')')
        write(*,'(1x,60(''=''))')
	
c    ** Zero run accumulators **
        call runzer

c    ** Perform simulation block-by-block **
        do block = 1, nblock
	   
c       ** Zero block accumulators **
           call blkzer

c       ** Perform required steps in block **
           do step = 1, nstep
	   
c          ** First part of velocity Verlet algorithm **
              call movea

c          ** Calculate forces and potential terms **	      
	      call force ( pot, vir )

c          ** Second part of velocity Verlet algorithm **
              call moveb
	      	      
c          ** Values for this step **
              call kinetic ( kin )
	      eng   = ( kin + pot ) / real(n)
	      temp  = 2.0 * kin / ( 3.0 * real(n) )
              dens  = real(n)/(box**3)
	      press = dens*temp+vir/(box**3)
	      
c          ** Accumulate block averages **
              call blkacc

           enddo
c       ** End of loop over steps in block

c       ** Normalize and print block averages **
           call blkout ( block )

c       ** Accumulate run averages **
           call runacc

        enddo
c    ** End of loop over blocks

c    ** Normalize and print run averages **
        write(*,'(1x,60(''=''))')
        call runout
        write(*,'(1x,60(''='')/)')
	
	return
	end


	
c    ===================================================================
        subroutine force ( pot, vir )
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Calculates forces, potential energy and virial                 |
c    |  Returned values:    pot           .... total potential energy  |
c    |                      vir           .... total virial            |
c    ===================================================================

c    ** Arguments **
        real            pot, vir

C    ** Local variables **
	integer i, j
	real            rxi, ryi, rzi, fxi, fyi, fzi
	real            rxij, ryij, rzij, rijsq, rcutsq, vcut
	real            r2ij, r6ij, r12ij, fij, vij
	real            fxij, fyij, fzij
	
c    -------------------------------------------------------------------

        rcutsq = rcut**2
	r2ij = 1.0/rcutsq
	r6ij = r2ij*r2ij*r2ij
	r12ij = r6ij**2
	vcut = r12ij - r6ij
	
C    ** Zero energies and forces **
        pot = 0.0
	vir = 0.0
	do i = 1, n
	   fx(i) = 0.0
	   fy(i) = 0.0
	   fz(i) = 0.0
	enddo

C    ** Outer loop over i **
        do i = 1, n-1
	   rxi = rx(i)
	   ryi = ry(i)
	   rzi = rz(i)
	   fxi = fx(i)
	   fyi = fy(i)
	   fzi = fz(i)
	   
C       ** Inner loop over j **
           do j = i+1, n
	   
c          ** Calculate interatomic vector **
	      rxij = rxi - rx(j)
	      ryij = ryi - ry(j)
	      rzij = rzi - rz(j)
	      rxij = rxij - anint(rxij/box)*box
	      ryij = ryij - anint(ryij/box)*box
	      rzij = rzij - anint(rzij/box)*box
              rijsq = rxij**2 + ryij**2 + rzij**2

c          ** Apply potential cutoff **		 
              if ( rijsq .lt. rcutsq ) then
		 r2ij = 1.0/rijsq
		 r6ij = r2ij*r2ij*r2ij
		 r12ij = r6ij*r6ij
		 vij = r12ij - r6ij
		 fij = (vij + r12ij)*r2ij
		 vij = vij - vcut
		 pot = pot + vij
		 fxij = fij*rxij
		 fyij = fij*ryij
		 fzij = fij*rzij
		 vir = vir + rxij*fxij + ryij*fyij + rzij*fzij
		 fxi = fxi + fxij
		 fyi = fyi + fyij
		 fzi = fzi + fzij
		 fx(j) = fx(j) - fxij
		 fy(j) = fy(j) - fyij
		 fz(j) = fz(j) - fzij
	      endif
	      
	   enddo
c       ** End of inner loop over j **

	   fx(i) = fxi
	   fy(i) = fyi
	   fz(i) = fzi
	   
	enddo
c    ** End of outer loop over i

c    ** Incorporate correct numerical factors **
        pot = pot * 4.0
	vir = vir * 24.0/3.0
        do i = 1, n
	   fx(i) = fx(i)*24.0
	   fy(i) = fy(i)*24.0
	   fz(i) = fz(i)*24.0
	enddo
	
	return
	end


	
c    ===================================================================
        subroutine kinetic ( kin )
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Calculates kinetic energy                                      |
c    |  Returned value:     kin             .... total kinetic energy  |
c    ===================================================================

c    ** Argument **
        real            kin

C    ** Local variables **
	integer i
        real    mass1
        parameter ( mass1 = 0.5 )
	
c    -------------------------------------------------------------------

	kin = 0.0
	
c    ** Atom 1 is special **
        kin = kin + mass1*(vx(1)**2 + vy(1)**2 + vz(1)**2)
	do i = 2, n
	   kin = kin + vx(i)**2 + vy(i)**2 + vz(i)**2
	enddo
	kin = kin / 2.0
	
	return
	end


	
c    ===================================================================
        subroutine movea
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  First part of velocity Verlet algorithm                        |
c    ===================================================================

C    ** Local variables **
	integer i
	real    dt2, dtsq2
        real    mass1
        parameter ( mass1 = 0.5 )
	
c    -------------------------------------------------------------------

        dt2 = dt / 2.0
	dtsq2 = dt*dt2
	
c    ** Atom 1 is special **
	rx(1) = rx(1) + dt*vx(1) + dtsq2*fx(1)/mass1
	ry(1) = ry(1) + dt*vy(1) + dtsq2*fy(1)/mass1
	rz(1) = rz(1) + dt*vz(1) + dtsq2*fz(1)/mass1
	vx(1) = vx(1) + dt2*fx(1)/mass1
	vy(1) = vy(1) + dt2*fy(1)/mass1
	vz(1) = vz(1) + dt2*fz(1)/mass1
	do i = 2, n
	   rx(i) = rx(i) + dt*vx(i) + dtsq2*fx(i)
	   ry(i) = ry(i) + dt*vy(i) + dtsq2*fy(i)
	   rz(i) = rz(i) + dt*vz(i) + dtsq2*fz(i)
	   vx(i) = vx(i) + dt2*fx(i)
	   vy(i) = vy(i) + dt2*fy(i)
	   vz(i) = vz(i) + dt2*fz(i)
	enddo
	
	return
	end


	
c    ===================================================================
        subroutine moveb
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Second part of velocity Verlet algorithm                       |
c    ===================================================================

C    ** Local variables **
	integer i
	
	real    dt2
        real    mass1
        parameter ( mass1 = 0.5 )
	
c    -------------------------------------------------------------------

        dt2 = dt / 2.0
	
c    ** Atom 1 is special **
	vx(1) = vx(1) + dt2*fx(1)/mass1
	vy(1) = vy(1) + dt2*fy(1)/mass1
	vz(1) = vz(1) + dt2*fz(1)/mass1
	do i = 2, n
	   vx(i) = vx(i) + dt2*fx(i)
	   vy(i) = vy(i) + dt2*fy(i)
	   vz(i) = vz(i) + dt2*fz(i)
	enddo
	
	return
	end
	
	

c    ===================================================================
        subroutine blkzer
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Sets to zero all block accumulators                            |
c    ===================================================================

c    -------------------------------------------------------------------

        blkeng = 0.0d0
        blkden = 0.0d0
        blkprs = 0.0d0
        blktmp = 0.0d0
	
	return
	end
	

c    ===================================================================
        subroutine runzer
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Sets to zero all run accumulators                              |
c    ===================================================================

c    -------------------------------------------------------------------

        runeng = 0.0d0
        runden = 0.0d0
        runprs = 0.0d0
        runtmp = 0.0d0
	
        erreng = 0.0d0
        errden = 0.0d0
        errprs = 0.0d0
        errtmp = 0.0d0
	
	return
	end

	

c    ===================================================================
        subroutine blkacc
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Increments all block accumulators                              |
c    ===================================================================

c    -------------------------------------------------------------------

        blkeng = blkeng + eng
        blkden = blkden + dens
        blkprs = blkprs + press
        blktmp = blktmp + temp
	
	return
	end

	

c    ===================================================================
        subroutine runacc
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Increments all run accumulators (including error estimators)   |
c    ===================================================================

c    -------------------------------------------------------------------

        runeng = runeng + blkeng
        runden = runden + blkden
        runprs = runprs + blkprs
        runtmp = runtmp + blktmp

        erreng = erreng + blkeng**2
        errden = errden + blkden**2
        errprs = errprs + blkprs**2
        errtmp = errtmp + blktmp**2
	
	return
	end


	

c    ===================================================================
        subroutine blkout ( block )
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Normalizes and prints out block averages                       |
c    |  Supplied argument: block  .............. current block number  |
c    ===================================================================

c    ** Arguments **
        integer block
	
c    -------------------------------------------------------------------

        blkeng = blkeng / dble(nstep)
        blkden = blkden / dble(nstep)
        blkprs = blkprs / dble(nstep)
        blktmp = blktmp / dble(nstep)
	
	write(*,'(1x,''Block  '',i3,i10,4f10.4)') 
     :     block, nstep, blkeng, blkden, blkprs, blktmp
     
	return
	end

	

c    ===================================================================
        subroutine runout
	implicit none
	include 'md.inc'
c    ===================================================================

c    ===================================================================
c    |  Normalizes and prints out run averages with error estimates    |
c    ===================================================================

c    -------------------------------------------------------------------

        runeng = runeng / dble(nblock)
        runden = runden / dble(nblock)
        runprs = runprs / dble(nblock)
        runtmp = runtmp / dble(nblock)

        erreng = erreng / dble(nblock)
        errden = errden / dble(nblock)
        errprs = errprs / dble(nblock)
        errtmp = errtmp / dble(nblock)
        erreng = ( erreng - runeng**2 ) / dble(nblock)
        errden = ( errden - runden**2 ) / dble(nblock)
        errprs = ( errprs - runprs**2 ) / dble(nblock)
        errtmp = ( errtmp - runtmp**2 ) / dble(nblock)
	if (erreng.gt.0.0d0) erreng=sqrt(erreng)
	if (errden.gt.0.0d0) errden=sqrt(errden)
	if (errprs.gt.0.0d0) errprs=sqrt(errprs)
	if (errtmp.gt.0.0d0) errtmp=sqrt(errtmp)
	
	write(*,'(1x,''Run avg.  '',i10,4f10.4)') 
     :     nblock*nstep, 
     :     runeng, runden, runprs, runtmp
	write(*,'(1x,''Error est.'',i10,4f10.4)') 
     :     nblock*nstep,
     :     erreng, errden, errprs, errtmp

	return
	end



c    ===================================================================
        subroutine getintg ( string, value )
	implicit none
c    ===================================================================

c    ===================================================================
c    |  Prompts for an integer value with a specified string           |
c    |  Echoes the supplied value in the output file                   |
c    |  Supplied arguments: string        .... string used in prompt   |
c    |  Returned values:    value         .... the input value         |
c    ===================================================================

c    ** Arguments **
        character*(*) string
	integer       value
	
c    -------------------------------------------------------------------

        print *, 'Enter ',string
	read  *, value
	print *, 'Value of ',string,' = ', value
	return
	end
	


c    ===================================================================
        subroutine getreal ( string, value )
	implicit none
c    ===================================================================

c    ===================================================================
c    |  Prompts for a  real value with a specified string              |
c    |  Echoes the supplied value in the output file                   |
c    |  Supplied arguments: string        .... string used in prompt   |
c    |  Returned values:    value         .... the input value         |
c    ===================================================================

c    ** Arguments **
        character*(*)   string
	real            value
	
c    -------------------------------------------------------------------

        print *, 'Enter ',string
	read  *, value
	print *, 'Value of ',string,' = ', value
	return
	end