module modulo1
implicit none
private
public :: interazione,kr
integer,parameter::kr=selected_real_kind(12)
contains

subroutine interazione(pos,nbody,upot,box,press,super,radius,hard)
real(kind=kr), intent(in), dimension(:,:) :: pos
integer, intent(in) :: nbody
real(kind=kr), intent(in) :: box
real(kind=kr), intent(out) :: upot,press
real(kind=kr), dimension(size(pos,1)) :: posij, deltaij
real(kind=kr) :: rij, radius
integer :: i,j
logical, intent(out) :: super
logical, intent(in) :: hard

upot = 0
press = 0
super=.false.
do i=1,nbody
do j=1,nbody
  if( i==j ) cycle
  posij = pos(:,i)-pos(:,j)
  deltaij = mod(posij,box)
  deltaij = deltaij-box*int(2*deltaij/box)
  rij=sqrt(dot_product(deltaij,deltaij))
  if (hard) then
    if (rij.lt.(2.*radius)) then
      super=.true.
      exit
    end if 
  else
    upot = upot + 4*(rij**(-12)-rij**(-6)) 
    press=press+24*(2.*rij**(-12)-rij**(-6))  
  end if
end do
end do
upot = upot/2
press = press/6

end subroutine interazione

end module modulo1

program corpi3d
use modulo1, rk=>kr
implicit none
integer :: nstep,naccept=0,nchosen,it,nbody,nsave,ios,i,j,ig 
integer,parameter :: nh=100
real(kind=rk), parameter :: pi=4.*atan(1.) 
real(kind=rk) :: xrand, mcrand
real(kind=rk), dimension(3) :: srand, cdm0, cdm
real(kind=rk) :: dt,mepot,mepot0,alfa=1.,box,frac,dens, rsq, rad, del, part, radius
real(kind=rk) :: temp, press, press0, pressave, density
real(kind=rk),dimension(:,:),allocatable :: pos, pos0, trialpos
real(kind=rk),dimension(:,:),allocatable :: d
real(kind=rk),dimension(nh) :: g=0.
real(kind=rk), dimension(3) :: pij, dij
real(kind=rk) :: gij
logical :: dyn,anneal=.false.,cont,super,hard

write(unit=*,fmt="(a)",advance="no")"n of bodies?"
read*, nbody
write(unit=*,fmt="(a)",advance="no")"how many it in output?"
read*, nsave
write(unit=*,fmt="(a)",advance="no")"Hard?"
read*, hard
if (hard) then
  write(unit=*,fmt="(a)",advance="no")"radius?"
  read*, radius
  write(unit=*,fmt="(a)",advance="no")"Volume fraction?"
  read*, frac
  box=radius*(nbody*4.*pi/(3.*frac))**(1./3.)
else
  write(unit=*,fmt="(a)",advance="no")"temperature?"
  read*, temp 
  write(unit=*,fmt="(a)",advance="no")"annealing?"
  read*, anneal
  if (anneal) then
    write(unit=*,fmt="(a)",advance="no")"scaling parameter?"
    read*, alfa
  end if
  write(unit=*,fmt="(a)",advance="no") "Number density?"
  read*, dens
  box=(nbody/dens)**(1./3.)
end if

write (unit=*,fmt=*) "Box:", box
write (unit=*, fmt=*) "Volume:", box**3
density=nbody/box**3
write (unit=*, fmt=*) "Density:", density
if (hard) write (unit=*, fmt=*) "Volume fraction:", density*4.*pi*radius**3/3.

write(unit=*,fmt="(a)",advance="no")"step: "
read*,dt
write(unit=*,fmt="(a)",advance="no")"n.step: "
read*,nstep
write(unit=*,fmt="(a)",advance="no") "Continuation run?"
read*, cont

allocate(pos(3,nbody))
allocate(pos0(3,nbody))
allocate(trialpos(3,nbody))
allocate(d(3,nbody))
del=box/(2*nh)

open (unit=1, file="posvel.d", status="replace")
if (.not.hard) open (unit=2, file="energies.d", status="replace")
open (unit=3, file="distances.d", status="replace")
open (unit=4, file="restart.d", status="unknown", form="unformatted")
open (unit=7, file="jmol.d", status="replace")
open (unit=8, file="squaredisp.d", status="replace")
open (unit=9, file="gofr.d", status="replace")
open (unit=10, file="inputcoord.d", status="old")
if (.not.hard) open (unit=22, file="pressure.d", status="replace")

if (cont) then
  read(unit=4) pos
else
  do
    read(unit=10,fmt=*,iostat=ios) pos 
    if (ios<=0) exit
  end do
  pos=pos*box
end if
do i=1,3
  cdm0(i) = sum(pos(i,:))/nbody
end do
pos0=pos
pressave=0.

call interazione(pos,nbody,mepot,box,press,super,radius,hard)
if (super) stop "Incorrect starting configuration"
time: do it = 1,nstep
  if (mod(it,10).eq.0) then
    write(unit=7,fmt=*) nbody
    write(unit=7,fmt=*)
  end if
  call random_number(xrand)
  nchosen=int(xrand*nbody)+1
  trialpos=pos
  call random_number(srand)
  trialpos(:,nchosen)=trialpos(:,nchosen)+dt*(srand-0.5)
  call interazione(trialpos,nbody,mepot0,box,press0,super,radius,hard)
  if (hard) then
    if (.not.super) then
      pos=trialpos
      naccept=naccept+1
    end if
  else
    call random_number(mcrand) 
    if ((mepot0.lt.mepot).or.(exp(-(mepot0-mepot)/temp).gt.mcrand)) then
      pos=trialpos
      mepot=mepot0
      press = (nbody*temp+press0)/box**3
      naccept=naccept+1
    end if
  end if
  do i=1,3
    cdm(i) = sum(pos(i,:))/nbody
  end do
  pos=pos+spread(cdm0-cdm,2,nbody)
  if (mod(it,nstep/nsave).eq.0) then
    write(unit=1,fmt=*)it,pos
    if (.not.hard) then 
      write(unit=2,fmt=*)it,mepot
      write(unit=22,fmt=*) it, press
    end if
  end if
  if (anneal) temp=alfa*temp
  rsq=sum((pos-pos0)**2)
  write(unit=8,fmt=*) it, rsq
  pressave=pressave+press

  if (mod(it,10).eq.0) then
    do i = 1,nbody
      d=mod(pos,box)
      where (d.lt.0.) d=d+box
      write (unit=7,fmt=*) "Ar", d(:,i)
    end do
  end if
  do i=1,nbody-1
    do j=i+1,nbody
      pij = pos(:,i)-pos(:,j)
      dij = mod(pij,box)
      dij = dij-box*int(2*dij/box)
      gij=sqrt(dot_product(dij,dij))
      if (gij.lt.box/2.) then
        ig=int(gij/del) 
        g(ig)=g(ig)+2
      end if
    end do
  end do
end do time

if (.not.hard) then
  write (unit=*,fmt=*) "Temperature:", temp
  write (unit=*,fmt=*) "Pressure:", pressave/nstep
end if
write (unit=*,fmt=*) "Acceptance ratio:", real(naccept)/nstep

do i=0,nh-1
  rad=del*(i+.5)
  part=4*pi*rad**2*del*nbody/box**3
  write(unit=9,fmt=*) rad, g(i+1)/(part*nbody*nstep)
end do
  
do i=1,nbody
do j=i+1,nbody
write(unit=3,fmt=*) i,j,sqrt(dot_product(pos(:,i)-pos(:,j),pos(:,i)-pos(:,j)))
end do
end do

write(unit=4) pos

close(1)
if (.not.hard) close(2)
close(3)
close(4)
close(7)
close(8)
close(9)
close(10)
close(11)
close(21)
if (.not.hard) close(22)

end program corpi3d