! by Dongxiang Liao
! 
!  Program to calculate radial distribution function from *.crd file
! (particals are on a sphere)  
!
program gofrhd
  implicit none
  integer :: ni, nsets, iset, iatom, ibin,j
  integer :: natoms, nbins
  real :: sigma, dr, fine_r
  real, dimension(512,3) :: r
  real, dimension(10000) :: gofr, cgofr
  character (LEN=8) :: filename
  print*, 'calculate RDF, from *.crd file.'
  print*, 'how fine, i.e., dr/sigma = ?'
  read*, fine_r
!  fine_r = 0.1
  print*, 'how many points'
!  nbins = 100
  read*, nbins
  print*, 'data from ni set, total nsets'
  print*, 'ni = ?'
  read*, ni
!  ni = 2
  print*, 'nsets = ?'
  read*, nsets
!  nsets = 1
  print*, 'file name?'
  read*, filename
!  filename = "105.crd"
  open(12,file=filename,action="read")
! read and discard
  do iset = 1, ni-1 
     read(12,'(I4, 2X, f10.8)') natoms, sigma
     do iatom=1,natoms
        read(12,'(3(f10.8,1X))' ) ( r(iatom,j), j=1,3)
     end do
  end do
! useful data
  
  do iset = ni, ni+nsets-1
     read(12,'(I4, 2X, f10.8)') natoms, sigma
     dr = sigma * fine_r
     do iatom = 1, natoms
        read(12,'(3(f10.8,1X))' ) ( r(iatom,j), j=1,3 )
!        print*, r(iatom,:)
     end do
 !    print*, dr, " ", nbins
     call RDF(r, dr, nbins, gofr,natoms)
     cgofr = cgofr + gofr
  end do
  close(12)
  cgofr = cgofr / real(nsets)
  open(14,file="gofr.dat")
  do ibin = 1, nbins
     write(14,*) 1.0+ibin*fine_r, "  ", cgofr(ibin)
  end do
  close(14)


CONTAINS
subroutine RDF(r, dr, nbins, gofr, natoms) 
    ! calculate the Radial Distribution Function -- gofr
    real, parameter :: pi= 3.141592653
    real, intent (in) :: dr
    real :: dij,  N_ideal, area_ring, dotp, rho
    integer :: i, j, bin, natoms 
    integer, intent (in) :: nbins
    real, dimension(512,3),intent(in) :: r
    real, dimension(nbins), intent(out) :: gofr 
    do i = 1, nbins
       gofr(i) = 0.0
    end do

    do i = 1, natoms-1
       do j = i+1, natoms
          dotp = dot_product( r(i,:), r(j,:))
          if (dotp < - 1.0) dotp = -0.999999
          if (dotp > 1.0 ) dotp = 0.99999
          dij = acos( dotp)
          bin = int((dij-sigma)/dr) + 1 
          if (bin < nbins) then
             gofr(bin) = gofr(bin) + 2.0
          end if
       end do
    end do
    ! normalization against ideal gas
    do bin = 1, nbins       
       area_ring = 2*pi * 2*sin(dr/2.0)*sin((bin-0.5)*dr)
       ! = ( cos((i-1)*dr) - cos(i*dr) ), try to rid of round off error
       rho = real(natoms) / (4.0 * 3.141592653)
       N_ideal = area_ring * rho 
       gofr(bin) = gofr(bin)/N_ideal/real(natoms)
!       print*, N_ideal," ", gofr(bin) 
   end do
 
 
  end subroutine RDF



end program gofrhd
