program kitty

    implicit none
    integer, parameter  :: n=120
    real*8  :: rx(n), ry(n), rz(n)
    real    :: k, l
    integer :: flag, intgr
    real*8  :: rxij, ryij, rzij, rijsq, rij
    real*8  :: delr, rupper, rlower, nideal, box, box3, const, dens, pi
    integer :: bin, maxbin, ncountgr, i, j, ii, hist(100000)
    real*8  :: gr1(100000),gr2(100000),gr3(100000),gr4(100000)

    open ( unit = 7, file = 'out.out', status = 'unknown' )

    pi  = 4.0d0*tan(1.0d0)
    box = 6.09d0
    write(*,*)"box in nm= ", box
    box3 = box*box*box
    write(*,*) "box3= ", box3
    dens = n/box3
    write(*,*)"density= ", dens
    ncountgr = 25001
    delr  = 0.04

!    maxbin  = int(0.5*box/delr)
    maxbin = 100005

    do i=1,maxbin
        hist(i)=0
    end do
    flag = 0

    do ii=1,ncountgr
        do i=1, n
                    flag = flag+1
            read(7,*) k,l,rx(i),ry(i),rz(i)
        end do

        do  i=1, n-1

            do j = i+1, n

                rxij =  rx(i)-rx(j)
                ryij =  ry(i)-ry(j)
                rzij =  rz(i)-rz(j)
                if (rxij>box .or. rxij<0.0) then
                    rxij = modulo(rxij, box)
                end if
                if (ryij>box .or. ryij<0.0) then
                    ryij = modulo(ryij, box)
                end if
                if (rzij>box .or. rzij<0.0) then
                    rzij = modulo(rzij, box)
                end if

                rijsq = rxij * rxij + ryij * ryij + rzij * rzij
                rij = dsqrt ( rijsq )
                bin = int ( rij/delr )
                hist (bin) = hist (bin) + 1
            end do
        end do
    end do
   write(*,*)"flag=",flag/ncountgr

! Calculation of solute solvent g(r)
! ==================================

    const=4.0*pi*dens/3.0   ! constant for calculating the volume
    intgr = 0
    do  bin = 1, maxbin
        rlower = real(bin)*delr
        rupper = rlower + delr
        nideal = const*( rupper**3 - rlower**3 )
        gr1(bin) = real(hist(bin))
        gr2(bin) = real(hist(bin))/real(ncountgr)
        intgr    = intgr  + hist(bin)
        gr3(bin) = real(hist(bin))/real(ncountgr)/real(n)
        gr4(bin) = real(hist(bin))/real(ncountgr)/real(n)/nideal
        write(120,*) rlower, gr1(bin), gr2(bin), gr3(bin), gr4(bin)
    end do
         write(*,*)"intgr=",intgr/ncountgr

    end program