!
! Simple MD program made for demonstration purposes
! 
!
! Internal units:
!
! All lengths in Å = 1e-10 m
! All times in fs = 1e-15 s
! All velocities and accelerations in Å/fs and Å/fs^2
! All forces in eV/Å
! All energies in eV
!
! Simulation cell centered at the origin.
!

program mdmorse
  USE GLOBAL
  implicit none
  
  ! Atom data containing arrays
  !
  ! x, v and a are the atom positions, velocities and accelerations
  ! Epot and Ekin are the potential and kinetic energy
  ! atomname is simply the atom name
  !
  type(vector)  :: x(MAXAT),v(MAXAT),a(MAXAT)
  real(double)  :: Epot(MAXAT),Ekin(MAXAT)
  character(5)  :: atomname(MAXAT)
  integer       :: neighbourlist(MAXAT*MAXNEIGHBOURS)
  
  ! Other parameters and variables
  !
  ! time is the current time, tmax the maximum time, deltat the time step
  !   and itime the integer number of the current time step
  ! T is the Temperature (in K) and m is the atom mass
  ! Ekinsum, Epotsum and Etot are the total energies for all atoms
  ! N is the number of atoms
  ! nneighbourlistupdate determines the interval between neighbour list
  !    updates, and nmovieoutput the interval between atom position output.
  ! rpotcut and rskincut are the potential and neighbour list cutoff radii
  ! size is the simulation cell size
  ! periodic determines whether periodic boundaries will be used (1=yes, 0=no)
  !
  real(double)  :: time,tmax,deltat,T,Ekinsum,Epotsum,Etot,m
  !
  integer       :: N,nneighbourlistupdate,itime,nmovieoutput
  real(double)  :: rpotcut,rskincut
  type(vector)  :: size
  type(ivector) :: periodic

  ! Morse potential parameters
  real(double) :: morseDe,morsealpha,morseRe

  ! Parameters for Berendsen temperature control [J. Chem. Phys.
  !  81 (1984) 3684]
  ! btctau is the Berendsen temperature control tau, initialT the initial
  ! temperature and desiredT the desired temperature.
  real(double) :: btctau,initialT,desiredT

  ! Parameters for Berendsen pressure control [J. Chem. Phys.
  !  81 (1984) 3684]
  ! bpctau is the Berendsen pressure control tau, bpcbeta 1/B,
  ! P the instantaneous pressure, and desiredP the desired pressure
  real(double) :: bpctau,bpcbeta,P,desiredP,virial,mu
  
  integer :: i

  print *,''
  print *,'--------------- mdmorse03 V1.2 --------------------'
  print *,''
  
  !
  ! Initialize the simulation
  !
  call ReadParams(initialT,m,size,tmax,periodic,  &
       deltat,nneighbourlistupdate,rpotcut,rskincut,  &
       morseDe,morsealpha,morseRe,nmovieoutput,    &
       btctau,T,bpctau,bpcbeta,desiredP)

  ! btc: Store desired T for temperature control
  desiredT=T;
  
  call ReadAtoms(x,atomname,N)
  
  ! Factor 2 because of energy equipartition theorem
  ! btc: changed T to initialT in call
  call SetTemperature(v,N,m,2.0*initialT)
  ! Check that we got it right
  call GetTemperature(v,N,m,T)
  print *,'Initial atom temperature is ',T

  ! Initialize accelerations to 0
  do i=1,N
     a(i)%x=0.0; a(i)%y=0.0; a(i)%z=0.0;
  enddo
  
  ! Start main loop over times
  time=0.0
  itime=0
  P=0.0
  do 
     if (mod(itime,nneighbourlistupdate)==0) then
	call UpdateNeighbourlist(x,N,size,periodic,neighbourlist,rskincut)
     endif
     
     call Solve1(x,v,a,N,size,periodic,deltat)
     
     call GetForces(x,a,N,size,periodic,m,neighbourlist,rpotcut, &
	  morseDe,morsealpha,morseRe,Epot,virial)

     call GetTemperature(v,N,m,T)
     
     call Solve2(v,a,N,deltat,btctau,T,desiredT)
     
     ! Update time
     time=time+deltat
     
     ! Actual solution done, now compute result quantities

     call GetEnergies(v,N,m,Ekinsum,Epotsum,Etot,Ekin,Epot)

     ! Compute pressure in units of kbar
     P=(N*kB*T+1.0/3.0*virial*e)/(size%x*size%y*size%z*1e-30)/1e8;

     ! Output pressure in units of kbar
     print '(A,F10.3,1X,F10.3,4F12.5)','ec ',time,T,Ekinsum/N,  &
	  Epotsum/N,Etot/N,P

     if (mod(itime,nmovieoutput)==0) then
	print '(A,F10.3)','Outputting atom movie at t = ',time
	call WriteAtoms(x,atomname,Ekin,Epot,N,time,size)
     endif

     ! Do Berendsen pressure control
     if (bpctau > 0.0) then
	mu=(1.0-bpcbeta*deltat/bpctau*(desiredP-P))**(1.0/3.0)

	size%x=size%x*mu
	size%y=size%y*mu
	size%z=size%z*mu
	do i=1,N
	   x(i)%x=x(i)%x*mu
	   x(i)%y=x(i)%y*mu
	   x(i)%z=x(i)%z*mu
	enddo
	if (mod(itime,10)==1) then
	   print '(A,3F10.6,3F12.5)','bpc ',size%x,size%y,size%z, &
		size%x*size%y*size%z,P,mu
	endif
     endif
	   
     ! Check whether we are done
     if (time >= tmax) exit     

     itime=itime+1

  end do

  ! Always output final energy and atom coordinates
  ! using high precision 
  print '(A,F12.5,1X,F12.5,3F14.7)','ec ',time,T,Ekinsum/N,Epotsum/N,Etot/N

  call WriteAtoms(x,atomname,Ekin,Epot,N,time,size)

  ! Close movie file 
  close(31)
     
end program mdmorse