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

subroutine interazione(pos,nbody,f,upot,box)
real(kind=kr), intent(in), dimension(:,:) :: pos
integer, intent(in) :: nbody
real(kind=kr), intent(in) :: box
real(kind=kr), intent(out) :: upot
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
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
end do
end do
upot = upot/2
end subroutine interazione
end module modulo1

program corpi3d
use modulo1, rk=>kr
implicit none
integer :: nstep,it,nbody,nsave,ios,i,j
real(kind=rk) :: dt,mepot,mekin,massa=1.,alfa=1.,vmax, box
real(kind=rk),dimension(:,:),allocatable :: pos
real(kind=rk),dimension(:),allocatable :: velcm
real(kind=rk),dimension(:,:),allocatable :: ekin,vel,f
logical :: dyn,anneal,cont

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(f(3,nbody))
allocate(velcm(3))

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=11) pos,vel
else
  do
    read(unit=10,fmt=*,iostat=ios) pos 
    if (ios<=0) exit
  end do
  if (dyn) then
!    do
!      read(unit=9,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

call interazione(pos,nbody,f,mepot,box)
do it = 1,nstep
  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)
    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)
    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
    end if
    if (anneal) vel=alfa*vel
  else
    do i = 1,nbody
      pos(:,i) = pos(:,i) + f(:,i)/massa * dt**2
    end do
    call interazione(pos,nbody,f,mepot,box)
    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
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

end program corpi3d