program kitty

    implicit none
    integer, parameter  :: n=120, natom = 40
    real*8  :: rx(n), ry(n), rz(n)
    real    :: k, l
    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)
    integer :: nn,mn,intgr
    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.02

    write(*,*)"which interaction?"
    read(*,*)nn, mn

    maxbin  = int(0.5*box/delr)

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

    do ii=1,ncountgr

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

        if(nn.eq.1.and.mn.eq.1) then
        do  i=1, n-1, 3
            do j = i+3, n, 3

!                if (ii.eq.1) then
!                write(*,*)i,j
!                end if
                rxij =  rx(i)-rx(j)
                ryij =  ry(i)-ry(j)
                rzij =  rz(i)-rz(j)

                if (rxij>box .or. rxij<0.0) then
                   rxij = rxij - (int(rxij/box))*box
                end if

                if (ryij>box .or. ryij<0.0) then
                    ryij = ryij - (int(ryij/box))*box
                end if

                if (rzij>box .or. rzij<0.0) then
                    rzij = rzij - (int(rzij/box))*box
                end if

                rijsq = rxij * rxij + ryij * ryij + rzij * rzij
                rij = dsqrt ( rijsq )
                bin = floor ( rij/delr )

                hist (bin) = hist (bin) + 1

            end do
        end do
        end if
    end do
! ==================================

    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) = hist(bin)
        gr2(bin) = hist(bin)/real(ncountgr)
        intgr = intgr + gr2(bin)
        gr3(bin) = hist(bin)/real(ncountgr)/real(natom)
        gr4(bin) = hist(bin)/real(ncountgr)/real(natom)/nideal
        write(120,*) rlower, gr1(bin), gr2(bin), gr3(bin), gr4(bin)
121     continue
    end do
    write(*,*)"integration of gr =", intgr

    end program