program mc  ! for Argon and Krypton
  ! An MC simulation for a cube of side L of liquid or gas.
  ! Units are eV and Angstroms. 
  use defs; use faces; implicit none
  character(len=20)::element
  character(len=20)::infile,outfile,start,plotfile,form
  character(len=20)::newdatafile,newinfile,newoutfile
  integer(i4b)::i,j,k,Natoms,Nsweeps,isweep,atom,Ndone
  integer(i4b)::nAccepts,run
  real(dp)::dx,duh1,Eold,Enew,myHeatVap
  real(dp)::rho,rhog,kT,beta,mol,heatVap,L
  real(dp),parameter::kB=8.617342D-05 ! eV/K
  real(dp),parameter::Joule = 1.0D19/1.602176462 ! eV
  real(dp),parameter::T0 = 273.15 ! Zero Celsius 
  real(dp),parameter::TbAr = -185.85 ! Celsius for Argon
  real(dp),parameter::TbKr = -153.22 ! Celsius for Kr (-153.4 Air Liquide)
  real(dp),parameter::TArBK = T0 + TbAr ! Boiling point of Argon (K)
  real(dp),parameter::TKrBK = T0 + TbKr ! Boiling point of Krypton (K)
  real(dp),parameter::kTAr = kB*TArBK
  real(dp),parameter::kTKr = kB*TKrBK
  real(dp),parameter::rhoAr = 1392.8D03 ! g/m^3 Density, Liquid @BP @1atm
  real(dp),parameter::rhoKr = 2413.0D03 ! g/m^3 Density, Liquid @BP @1atm
  real(dp),parameter::rhogAr = 5853.0 ! g/m3 density, gas @BP @1atm
  real(dp),parameter::rhogKr = 8520.0 ! g/m3 density, gas @BP @1atm
  real(dp),parameter::molAr = 39.948 ! grams
  real(dp),parameter::molKr = 83.80 ! grams
  real(dp),parameter::NA = 6.02214199D23 ! atoms per mol
  real(dp),parameter::MheatVapAr = 160.81 ! J/g = kJ/kg
  real(dp),parameter::MheatVapKr = 107.81 ! J/g
  real(dp),parameter::heatVapAr = MheatVapAr*Joule*molAr/NA
  real(dp),parameter::heatVapKr = MheatVapKr*Joule*molKr/NA
  real(dp),dimension(3)::xold,xnew
  real(dp)::density,T ! degrees Kelvin
  real(dp)::potEnergy,enthalpy,sigma
  real(dp),dimension(:),allocatable::sweepEnergy
  real(dp),dimension(:,:),allocatable::x,duh
  ! Read in or initialize data 
  read(5,*);read(5,*)element,form,run
  read(5,*);read(5,*)start, Natoms, Nsweeps, Ndone
  read(5,*);read(5,*)dx, infile, outfile, plotfile
  read(5,*);read(5,*)newinfile, newoutfile, newdatafile
  allocate(x(3,Natoms),duh(3,Natoms),sweepEnergy(Nsweeps))
  ! Program starts with mol wt, density, beta, & radius:
  if ( element == 'Ar' ) then
     mol = molAr; heatVap = heatVapAr
     rho = rhoAr; rhog = rhogAr
     kT = kTAr; beta = 1./kT
     density = (rho*NA)/(mol*1.0d30) ! .02099639 atoms/A^3 for Ar
  else if ( element == 'ArGas' ) then
     mol = molAr; heatVap = heatVapAr
     rho = rhoAr; rhog = rhogAr; rho = rhog
     kT = kTAr; beta = 1./kT
     density = (rho*NA)/(mol*1.0d30) ! 0.000088234 atoms/A^3 for Ar gas
  else if ( element == 'Kr' ) then
     mol = molKr; heatVap = heatVapKr
     rho = rhoKr; rhog = rhogKr
     kT = kTKr; beta = 1./kT
     density = (rho*NA)/(mol*1.0d30) ! 0.01734 atoms/A^3 for Kr
  else if ( element == 'KrGas' ) then
     mol = molKr; heatVap = heatVapKr
     rho = rhoKr; rhog = rhogKr; rho = rhog
     kT = kTKr; beta = 1./kT
     density = (rho*NA)/(mol*1.0d30) ! 0.000061228 atoms/A^3 for Kr gas
  else
     open(99,file='badElement')
     write(99,*)'Element',element,' is inappropriate.'
     close(99); stop
  end if
  L = (Natoms/density)**(1.0/3.0) ! length of box
  call random_seed()
  if ( start == 'cold' ) then
     call random_number(duh)
     x = L*duh
  else
     open(4,file=infile)
     read(4,*) x
  end if
  potEnergy = 0.; sweepEnergy = 0.; nAccepts = 0
  main: do isweep = 1, Nsweeps
     call random_number(duh)
     drop: do atom = 1, Natoms
        xold(:) = x(:,atom)
        xnew(:) = mod(xold(:) + dx*(duh(:,atom) - 0.5)+L,L)
        ! Do MC step
        call energy(element,form,atom,Natoms,xold,x,L,Eold)
        call energy(element,form,atom,Natoms,xnew,x,L,Enew)
        if ( Enew < Eold ) then  ! accept
           x(:,atom) = xnew; nAccepts = nAccepts + 1
           potEnergy = potEnergy + Enew
           sweepEnergy(isweep) = sweepEnergy(isweep) + Enew
        else ! accept conditionally
           call random_number(duh1)
           if ( duh1 < exp(-beta*(Enew - Eold)) ) then
              x(:,atom) = xnew; nAccepts = nAccepts + 1 ! accept
              potEnergy = potEnergy + Enew
              sweepEnergy(isweep) = sweepEnergy(isweep) + Enew
           else
              x(:,atom) = xold ! reject
              potEnergy = potEnergy + Eold
              sweepEnergy(isweep) = sweepEnergy(isweep) + Eold
           end if
        end if
     end do drop
     if ( mod(isweep,100) == 0 ) then
        open(7,file=trim(outfile))
        do i = 1, Natoms
           write(7,*)x(:,i)
        end do
        close(7)
     end if
  end do main
  potEnergy = 0.5*potEnergy/real(Natoms*Nsweeps,dp)
  sweepEnergy = kT*(1.0 - rhog/rho) - 0.5*sweepEnergy/real(Natoms,dp)
  open(7,file=trim(outfile))
  do i = 1, Natoms
     write(7,*)x(:,i)
  end do
  close(7)
  Ndone = Nsweeps + Ndone
  myHeatVap = -potEnergy + kT*(1.0 - rhog/rho)
  call stat(Nsweeps,sweepEnergy,enthalpy,sigma)
  open(8,file=trim(newdatafile))
  write(8,'(a)')'element     form   run'
  write(8,'(a,5x,a,i8)')trim(element), trim(form), run+1
  write(8,'(a)')'start  Natoms  Nsweeps  Ndone'
  write(8,'(a,3(i10))')'hot', Natoms, Nsweeps, Ndone
  write(8,'(a)')'dx   infile   outfile   plotfile'
  write(8,'(f7.3,3(2x,a,i4))')dx,'Xfile',run+1,'Xfile',run+2,'plot',run+1
  write(8,'(a)') 'newinfile   newoutfile   newdatafile'
  write(8,'(3(2x,a,i4))')'Xfile',run+2,'Xfile',run+3,'data',run+2
  write(8,'(a,f10.7)')'In the last run, PE per atom was:',potEnergy
  write(8,'(a,f10.7,a,f10.6)')'In the last run, totalHeatVap per atom was:',&
       enthalpy,' +-',sigma
  write(8,'(a,f9.3,a,2x,a,f10.7)')'% error = ',100.*(enthalpy-heatVap)/heatVap,'%',&
       'acceptance ratio =',real(nAccepts,dp)/real(Natoms*Nsweeps,dp)
  close(8)
  open(9,file=trim(plotfile))
  do i = 1, Nsweeps
     write(9,*)i, sweepEnergy(i)
  end do
  close(9)
  open(10,file=trim(plotfile)//'stats')
  write(10,'(a,f12.5,a,f12.9)')'Length of box =', L,' density =',density
  write(10,'(a,f10.7)')'experimental latent heat of vap =',heatVap
  write(10,'(a,f10.7,a,f10.7)')'enthalpy =', enthalpy,' +-',sigma
  write(10,'(a,f10.7)')'acceptance ratio =', &
       real(nAccepts,dp)/real(Natoms*Nsweeps,dp)
end program mc