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

c    ===================================================================
c    |  Molecular dynamics of Lennard-Jones atoms                      |
c    ===================================================================

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

c    **declarations **
c    -------------------------------------------------------------------

        call start
	call runner
	call finish
	
	stop
	end
	
	
	
c    ===================================================================
        subroutine start
	implicit none
	include 'md2.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*8            pot, vir, kin, ptot, pimp,timp

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

	open(unit=15,file='mts.dat',status='unknown')
	open(unit=16,file='mtsn.dat',status='unknown')
	
c    ** Output program identification **	
        print *, '-----------------------------------------------------'	
	print *, 'Molecular dynamics simulation of Lennard-Jones atoms'
	print *, 'Cubic periodic boundary conditions'
	print *, 'Canonical ensemble 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   )
	call getreal ( 'temp. dec. rate',nut )
	call getreal ( 'required temp.',extem)


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
	call kinetic(kin, ptot)
c	timp=ptot
	close(unit=3)

c    ** Set initial Nose'variables **

	eta=1.0   
	beta=0.0
	
c    ** Set initial impulse equal to zero **

	pimp=ptotx**2+ptoty**2+ptotz**2
	if(pimp.ne.0.0)then
c   	timp=pimp 
	 do i=1,n
        vx(i)=vx(i)-ptotx
        vy(i)=vy(i)-ptoty 
        vz(i)=vz(i)-ptotz
        enddo
	call kinetic(kin, ptot)
	endif


c    ** Output initial values **
10	call force ( pot, vir )
	call kinetic ( kin,ptot )
	eng   = ( kin + pot ) / real(n)
	temp  = 2.0*kin/(3.0*(real(n)-1))
        dens  = real(n)/(box**3)
	press = dens*temp+vir/(box**3)
	co =(real(n)-1)*1.5*extem*(beta/nut)**2+3.0*(real(n)-1)*extem*eta
	co =co/real(n)
	nco=eng+co
	
	 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
	write(*,'(1x,''Const.of motion = '',f10.4)') nco	
c 	write(*,'(1x,''Scal imp(^2)= '',e10.4)') ptot	
	write(*,'(1x,''px= '',e10.4)') ptotx	
	write(*,'(1x,''py= '',e10.4)') ptoty
	write(*,'(1x,''pz= '',e10.4)') ptotz
	end


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

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

C    ** Local variables **
	integer i
	real*8           pot, vir, kin, ptot, pimp

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

c    ** Output final values **
	call force ( pot, vir )
	call kinetic ( kin, ptot )
	eng   = ( kin + pot ) / real(n)
        temp  = 2.0*kin/( 3.0 * (real(n)-1) )
        dens  = real(n)/(box**3)
	press = dens*temp+vir/(box**3)
	co =(real(n)-1)*1.5*extem*(beta/nut)**2+3.0*(real(n)-1)*extem*eta  
        co=co/real(n)
	nco=eng+co  
 
        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
	write(*,'(1x,''Cons.of motion= '',f10.4)') nco
	write(*,'(1x,''px= '',e10.4)') ptotx
        write(*,'(1x,''py= '',e10.4)') ptoty
        write(*,'(1x,''pz= '',e10.4)') ptotz
	
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)
	close(15)

	
        print *, '-----------------------------------------------------'	
	print *, 'Program ends'
        print *, '-----------------------------------------------------'	

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

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

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

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



	write(*,'(1x,70(''=''))')
	
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, ptot )
	      eng   = ( kin + pot ) / real(n)
	      temp  = 2.0 * kin / ( 3.0 * (real(n)-1) )
              dens  = real(n)/(box**3)
	      press = dens*temp+vir/(box**3)
          co=(real(n)-1)*1.5*extem*(beta/nut)**2 &+3.0*(real(n)-1)*extem*eta 
	      co=co/real(n) 
	      nco= eng +co
	  write(15,*)kin/real(n),pot/real(n),eng
    	  write (16,*)eta,feta,nco 

 
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,70(''=''))')
        call runout
        write(*,'(1x,70(''='')/)')
	
	return
	end


	
c    ===================================================================
        subroutine force ( pot, vir )
	implicit none
	include 'md2.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*8           pot, vir

C    ** Local variables **
	integer i, j
	real*8           rxi, ryi, rzi, fxi, fyi, fzi
	real*8           rxij, ryij, rzij, rijsq, rcutsq, vcut
	real*8           r2ij, r6ij, r12ij, fij, vij
	real*8           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, ptot )
	implicit none
	include 'md2.inc'
c    ===================================================================

c    ===================================================================
c    |  Calculates kinetic energy and sq of total impulse              |
c    |  Returned value:     kin             .... total kinetic energy  |
c    |                      ptot            ....square of tot impulse  |
c    ===================================================================

c    ** Argument **
        real*8           kin, ptot

C    ** Local variables **
	integer i
c	real*8 ptotx,ptoty,ptotz
	
c    -------------------------------------------------------------------

	kin = 0.0
	ptot=0.0
	ptotx=0.0
	ptoty=0.0
	ptotz=0.0	

	do i = 1, n
	   kin = kin + vx(i)**2 + vy(i)**2 + vz(i)**2
	   ptotx=ptotx+vx(i)/real(n)
	   ptoty=ptoty+vy(i)/real(n)
	   ptotz=ptotz+vz(i)/real(n)		

       
	 enddo
	kin = kin / 2.0
	ptot=ptotx**2+ptoty**2+ptotz**2
	
	return
	end


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

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

C    ** Local variables **
	integer i
	real*8   dt2,dtsq2 , nutsq 

c    -------------------------------------------------------------------
        dt2 = dt / 2.0
	dtsq2 = dt*dt2
	nutsq=nut**2
	feta=(temp/extem-1)    

	etadot= beta+dt2*nutsq*feta
	es=etadot*dt2
	do i = 1, n
	rx(i) = rx(i) + dt*(vx(i)+dt*fx(i)*0.25)*exp(-es)+ dtsq2*0.5*fx(i)
	ry(i) = ry(i) + dt*(vy(i)+dt*fy(i)*0.25)*exp(-es)+ dtsq2*0.5*fy(i)
	rz(i) = rz(i) + dt*(vz(i)+dt*fz(i)*0.25)*exp(-es)+ dtsq2*0.5*fz(i)
	vx(i) = (vx(i)+dt*0.25*fx(i))*exp(-2.0*es)+ 0.5*dt2*exp(-es)*fx(i)
	vy(i) = (vy(i)+dt*0.25*fy(i))*exp(-2.0*es)+0.5 *dt2*exp(-es)*fy(i)
	vz(i) = (vz(i)+dt*0.25*fz(i))*exp(-2.0*es)+0.5* dt2*exp(-es)*fz(i)
	enddo


	eta=eta+dt*beta+dtsq2*nutsq*feta
	beta=beta+dt2*nutsq*feta	
	return
	end


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

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

C    ** Local variables **
	integer i
	
	real*8   dt2, nutsq
	real*8   kin, ptot
c    **Subprograms referenced** 
	external kinetic	
c    -------------------------------------------------------------------
	nutsq=nut**2
        dt2 = dt / 2.0
	do i = 1, n
	   vx(i) = vx(i) +0.5* dt2*exp(-es)*fx(i)+0.25*dt*fx(i)
	   vy(i) = vy(i) +0.5* dt2*exp(-es)*fy(i)+0.25*dt*fy(i)
	   vz(i) = vz(i) +0.5* dt2*exp(-es)*fz(i)+0.25*dt*fz(i)
	enddo
	call kinetic(kin, ptot)
	temp=2.0*kin/(3.0*(real(n)-1))
	feta=(temp/extem-1)
	beta=beta+dt2*nutsq*feta
	return
	end
	
	

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

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

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

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

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

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

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

  	runcons= 0.0d0 
        runeng = 0.0d0
        runden = 0.0d0
        runprs = 0.0d0
        runtmp = 0.0d0

	errcons= 0.0d0	
        erreng = 0.0d0
        errden = 0.0d0
        errprs = 0.0d0
        errtmp = 0.0d0
	
	return
	end

	

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

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

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

        blkeng = blkeng + eng
        blkden = blkden + dens
        blkprs = blkprs + press
        blktmp = blktmp + temp
	blkcons= blkcons+nco

	
	return
	end

	

c    ===================================================================
        subroutine runacc
	implicit none
	include 'md2.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
	runcons= runcons+ blkcons

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


	

c    ===================================================================
        subroutine blkout ( block )
	implicit none
	include 'md2.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)
	blkcons= blkcons/dble(nstep)

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

	

c    ===================================================================
        subroutine runout
	implicit none
	include 'md2.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)
	runcons= runcons/ dble(nblock)

        erreng = erreng / dble(nblock)
        errden = errden / dble(nblock)
        errprs = errprs / dble(nblock)
        errtmp = errtmp / dble(nblock)
      	errcons = errcons/ 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)
	errcons= ( errcons - runcons**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)
	if (errcons.gt.0.0d0) errcons=sqrt(errcons)	
	write(*,'(1x,''Run avg.  '',i10,5f10.4)') 
     :     nblock*nstep, 
     :     runeng, runden, runprs, runtmp, runcons
	write(*,'(1x,''Error est.'',i10,5f10.4)') 
     :     nblock*nstep,
     :     erreng, errden, errprs, errtmp, errcons

	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*8           value
	
c    -------------------------------------------------------------------

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