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

subroutine interazione(pos,nbody,f,upot,box,press)
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), intent(out), dimension(:,:) :: f
real(kind=kr), dimension(size(pos,1)) :: posij, deltaij
real(kind=kr) :: rij
integer :: i,j

upot = 0
f = 0
press = 0
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))
  upot = upot + 4*(rij**(-12)-rij**(-6)) 
  f(:,i) = f(:,i) + 24*(2.*rij**(-14)-rij**(-8))*deltaij
  press=press+24*(2.*rij**(-12)-rij**(-6))  
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,it,nbody,nsave,ios,i,j,ig 
integer,parameter :: nh=100
real(kind=rk), parameter :: pi=4.*atan(1.) 
real(kind=rk) :: dt,mepot,mekin,massa=1.,alfa=1.,vmax, box, rsq, rad, del, part
real(kind=rk) :: temp, tempave, press, pressave
real(kind=rk),dimension(:,:),allocatable :: pos, pos0
real(kind=rk),dimension(:),allocatable :: velcm, vbox
real(kind=rk),dimension(:,:),allocatable :: d
real(kind=rk),dimension(:,:),allocatable :: ekin,vel,f
real(kind=rk),dimension(nh) :: g=0.
real(kind=rk), dimension(3) :: pij, dij
real(kind=rk) :: gij
logical :: dyn,anneal,cont

open (unit=1, file="posvel.d", status="replace")
open (unit=2, file="energies.d", status="replace")
open (unit=3, file="distances.d", status="replace")
open (unit=4, file="restart.d", status="replace", 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")
!open (unit=11, file="inputvel.d", status="old")
open (unit=21, file="temperature.d", status="replace")
open (unit=22, file="pressure.d", status="replace")


write(unit=*,fmt="(a)",advance="no")"box?"
read*, box
write(unit=*,fmt="(a)",advance="no")"dynamics?"
read*, dyn
if (dyn) then
  write(unit=*,fmt="(a)",advance="no")"annealing?"
  read*, anneal
  if (anneal) then
    write(unit=*,fmt="(a)",advance="no")"scaling parameter?"
    read*, alfa
  end if
end if
write(unit=*,fmt="(a)",advance="no")"n of bodies?"
read*, nbody
write(unit=*,fmt="(a)",advance="no")"how many it in output?"
read*, nsave

allocate(vel(3,nbody))
allocate(ekin(3,nbody))
allocate(pos(3,nbody))
allocate(pos0(3,nbody))
allocate(f(3,nbody))
allocate(velcm(3))
allocate(d(3,nbody))
vbox=spread(box,1,size(pos,1))
del=box/(2*nh)

write(unit=*,fmt="(a)",advance="no")"time step: "
read*,dt
write(unit=*,fmt="(a)",advance="no")"n.step: "
read*,nstep
!print*,"mass of the bodies: "
!read*,massa
write(unit=*,fmt="(a)",advance="no") "Continuation run?"
read*, cont

vel=0.
if (cont) then
  read(unit=4) pos,vel
else
  do
    read(unit=10,fmt=*,iostat=ios) pos 
    if (ios<=0) exit
  end do
  pos=pos*box
  if (dyn) then
!    do
!      read(unit=11,fmt=*,iostat=ios) vel 
!      if (ios<=0) exit
!    end do
    call random_number(vel)
    write(unit=*,fmt="(a)",advance="no")"vmax? "
    read*, vmax 
    vel=2*vmax*(vel-0.5)
    do i=1,3
      velcm(i) = sum(vel(i,:))/nbody
      vel(i,:) = vel(i,:) - velcm(i)
    end do
  end if
end if
write (unit=*, fmt=*) "Volume:", box**3
write (unit=*, fmt=*) "Density:", nbody/box**3
pos0=pos
tempave=0.
pressave=0.

call interazione(pos,nbody,f,mepot,box,press)
time: do it = 1,nstep
  write(unit=7,fmt=*) nbody
  write(unit=7,fmt=*)
  if (dyn) then
    do i = 1,nbody
      pos(:,i) = pos(:,i) + vel(:,i) * dt + 0.5* f(:,i)/massa * dt**2
      vel(:,i) = vel(:,i) + 0.5 * dt * f(:,i)/massa
    end do
    call interazione(pos,nbody,f,mepot,box,press)
    do i=1,nbody
      vel(:,i) = vel(:,i) + 0.5 * dt * f(:,i)/massa
      ekin(:,i) = 0.5 * massa * (vel(:,i))**2
    end do
    mekin = sum(ekin)
    temp = 2*mekin/(3*nbody)
    press = (nbody*temp+press)/box**3
    if (mod(it,nstep/nsave).eq.0) then
      write(unit=1,fmt=*)it,it*dt,pos,vel
      write(unit=2,fmt=*)it,mekin,mepot,mekin+mepot
      write(unit=21,fmt=*) it, temp
      write(unit=22,fmt=*) it, press
    end if
    if (anneal) vel=alfa*vel
    rsq=sum((pos-pos0)**2)
    write(unit=8,fmt=*) it, rsq
    tempave=tempave+temp
    pressave=pressave+press
  else
    do i = 1,nbody
      pos(:,i) = pos(:,i) + f(:,i)/massa * dt**2
    end do
    call interazione(pos,nbody,f,mepot,box,press)
    if (mod(it,nstep/nsave).eq.0) write(unit=1,fmt=*)it,pos
    if (mod(it,nstep/nsave).eq.0) write(unit=2,fmt=*)it,mepot
  end if
  do i = 1,nbody
    d=mod(pos,box)
    where (d.lt.0.) d=d+box
    write (unit=7,fmt=*) "Ar", d(:,i)
  end do
  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

write (unit=*,fmt=*) "Temperature:", tempave/nstep
write (unit=*,fmt=*) "Pressure:", pressave/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,vel

close(1)
close(2)
close(3)
close(4)
close(7)
close(8)
close(9)
close(10)
close(11)
close(21)
close(22)

end program corpi3d