! Input-output subroutines for the mdmorse code

subroutine ReadParams(initialT,m,size,tmax,periodic, &
     deltat,nneighbourlistupdate,rpotcut,rskincut,  &
     morseDe,morsealpha,morseRe,nmovieoutput, &
     btctau,desiredT,bpctau,bpcbeta,desiredP)
     
  USE GLOBAL
  implicit none

  real(double)  :: initialT,m,tmax
  type(vector)  :: size
  type(ivector) :: periodic
  real(double)  :: deltat,rpotcut,rskincut,morseDe,morsealpha,morseRe
  integer       :: nneighbourlistupdate,nmovieoutput
  real(double)  :: btctau,desiredT,bpctau,bpcbeta,desiredP
  
  integer :: i,ios
  character(80) :: line
  character(40) :: string
  real(double) :: x
  
  open(20,file='mdmorse.in',status='old',iostat=ios)
   
  if (ios /= 0) STOP 'Open error for file mdmorse.in'

  ! Initialize variables to sensible (?) default values
  initialT=0.0; m=63.546;
  size%x=0.0; size%y=0.0; size%z=0.0;

  deltat=1.0;
  nneighbourlistupdate=5;
  rpotcut=4.02; rskincut=5.0;
  ! Initialize Morse parameters for Cu
  ! Values from Girifalco and Weizer, Phys. Rev. 114 (1959) 687.
  morseDe=0.3429; morsealpha=1.3588; morseRe=2.866;

  btctau=0.0; desiredT=0.0;
  bpctau=0.0; bpcbeta=1e-3; desiredP=0.0;
  
  seed=23267761
  
  ! Loop through file and find all variables
  i=0
  do
     i=i+1

     ! First read in line in whole.
     read(20,fmt='(A80)',iostat=ios) line
     ! Check for probable end-of-file
     if (ios < 0) exit

     ! If not parameter line, cycle
     if (line(1:1) /= '-') cycle
     
     ! Attempt to parse line into command and value
     ! using internal formatted read
     read(unit=line,fmt=*,iostat=ios) string,x
     if (ios > 0) then
        print *,'ERROR: Data read in error on line',i
        stop 'Invalid parameter file'
     endif

     ! From here on line should be in variable format
     
     ! Look for variable-defining strings
     if (string=='-initialT') then
	initialT=x
     else if (string=='-desiredT') then
	desiredT=x
     else if (string=='-btctau') then
	btctau=x
     else if (string=='-mass') then
	m=x
     else if (string=='-xsize') then
	size%x=x
     else if (string=='-ysize') then
	size%y=x 
     else if (string=='-zsize') then
	size%z=x
     else if (string=='-periodicx') then
	periodic%x=int(x+0.5)
     else if (string=='-periodicy') then
	periodic%y=int(x+0.5)
     else if (string=='-periodicz') then
	 periodic%z=int(x+0.5)
     else if (string=='-seed') then
	seed=int(x+0.5)
     else if (string=='-tmax') then
	tmax=x
     else if (string=='-deltat') then
	deltat= x
     else if (string=='-nupdate') then
	nneighbourlistupdate=int(x+0.5)
     else if (string=='-rpotcut') then
	rpotcut=x 
     else if (string=='-rskincut') then
	rskincut=x     
     else if (string=='-morseDe') then
	morseDe=x
     else if (string=='-morsealpha') then
	morsealpha=x
     else if (string=='-morseRe') then
	 morseRe=x
     else if (string=='-bpctau') then
	 bpctau=x
     else if (string=='-bpcbeta') then
	 bpcbeta=x
     else if (string=='-desiredP') then
	 desiredP=x
     else if (string=='-nmovieoutput') then
	nmovieoutput=int(x+0.5)
     else
	print '(A,A)','Unknown parameter',string
	stop 'Parameter read in error'
     endif

     print '(A,A16,A,G13.6)','Read in parameter ',string,' value',x
    
  enddo

  close(20)

  ! Print out some general information about read results:

  print '(A,3I2)','Using periodics (1=on, 0=off)',     &
       periodic%x,periodic%y,periodic%z

  print '(A,3F12.6)','Morse potential parameters: De alpha Re', &
       morseDe,morsealpha,morseRe

  print '(A,I8,A)','Movie output selected every',nmovieoutput,' steps'

  if (btctau > 0.0) then 
     print '(/,A,2F10.3,/)','Doing Berendsen temperature control with tau T', &
	  btctau,desiredT
  endif
  
  if (bpctau > 0.0) then 
     print '(/,A,2F10.3,/)','Doing Berendsen pressure control with tau beta', &
	  bpctau,bpcbeta
  endif
  
  return
  
end subroutine ReadParams


subroutine ReadAtoms(x,atomname,N)
  USE GLOBAL
  implicit none

  integer :: N
  type(vector) :: x(*)
  character(5) atomname(*)
  
  character(80) :: comment
  integer :: i,ios,type

  open(30,file='atoms.in',status='old',iostat=ios)

  if (ios/=0) STOP 'File atoms.in open failed'

  read(30,fmt=*) N

  if (N > MAXAT) STOP 'Error: MAXAT overflow, increase it.'

  read(30,fmt='(A80)') comment
  print '(A,I7,A,A40)','Reading in ',N,' atoms described as ',comment
  
  do i=1,N

     read(30,fmt=*,iostat=ios) atomname(i),x(i)%x,x(i)%y,x(i)%z,type

     if (ios/=0) then
	print *,'Atom read in error for atom',i
	print *,atomname(i),x(i)%x,x(i)%y,x(i)%z,type
	STOP ' Atom read in error'
     endif

  enddo
  
  
end subroutine ReadAtoms



subroutine WriteAtoms(x,atomname,Ekin,Epot,N,time,size)
  USE GLOBAL
  implicit none

  type(vector) :: x(*)
  character(5) :: atomname(*)
  integer :: N
  real(double) :: Ekin(*),Epot(*),time
  type(vector) :: size
  
  character(80) :: comment
  integer :: i,ios,type

  logical, save :: firsttime=.true.

  if (firsttime) then
     open(31,file='atoms.out',status='replace',iostat=ios)
     if (ios/=0) STOP 'File atoms.out open failed'
     
     firsttime=.false.
  endif
     

  write(31,fmt='(I10)') N

  write(31,fmt='(A,F11.3,A,3F12.4)') 'mdmorse atom output at time ', &
       time,' fs boxsize ',size%x,size%y,size%z
  
  do i=1,N
     type=1
     write(31,fmt='(A2,1X,3F12.6,2F11.6)') atomname(i),x(i)%x,x(i)%y,x(i)%z, &
	  Ekin(i),Epot(i)
  enddo
  
end subroutine WriteAtoms