PROGRAM hydrodynamic
IMPLICIT NONE
INTEGER, PARAMETER :: natoms = 80
CHARACTER(len=20) :: junk
INTEGER, PARAMETER :: snapshots = 19502
INTEGER :: natom
CHARACTER(len=2) :: atomlabel
REAL :: x,y,z
REAL :: bondl
INTEGER :: no_of_pairs
REAL, DIMENSION(snapshots) :: sum_distance
REAL :: local_distance
INTEGER :: i,j,k
TYPE coords
REAL, DIMENSION(3) :: pos
END TYPE coords

REAL :: ensemble_average
REAL :: hydro

TYPE(coords), DIMENSION(natoms,snapshots) :: coordinates

OPEN(UNIT=10, FILE="sop_prod_7331.xyz")
OPEN(UNIT=11, FILE="hydrodynamic_7331.txt")

!Read coordinates
DO j = 1, snapshots
READ(10,*) natom
READ(10,*) junk
DO i = 1, natoms
READ(10,*) atomlabel, x,y,z
coordinates(i,j)%pos(1) = x
coordinates(i,j)%pos(2) = y
coordinates(i,j)%pos(3) = z
ENDDO
ENDDO

no_of_pairs = 0
sum_distance(:) = 0.
ensemble_average = 0.
!Calculate pairwise distances
DO k = 1, snapshots
no_of_pairs = 0
DO i = 1, natoms-1
DO j = i+1, natoms
IF(i.eq.j)CYCLE
local_distance = bondl(coordinates(i,k)%pos(:),coordinates(j,k)%pos(:))
!write(*,*) local_distance, 1/local_distance
sum_distance(k) = sum_distance(k) + (1/local_distance)
no_of_pairs = no_of_pairs + 1
ENDDO
ENDDO
!hydro = sum_distance(k)/(no_of_pairs)
!write(*,*) "Pairs", no_of_pairs
hydro = (sum_distance(k)/no_of_pairs)
WRITE(11,*) k, 1/hydro
!no_of_pairs = 0
ENDDO

!DO k = 1, 19502
!ensemble_average = ensemble_average + sum_distance(k)
!ENDDO





!ensemble_average = ensemble_average/19502
!hydro = ensemble_average/real(natoms*natoms)

!WRITE(*,*) 1/hydro


CLOSE(UNIT=10)
CLOSE(UNIT=11)




END PROGRAM



REAL FUNCTION bondl(r1,r2)
IMPLICIT NONE
REAL, DIMENSION(3), INTENT(IN) :: r1
REAL, DIMENSION(3), INTENT(IN) :: r2
REAL :: distance_x, distance_y, distance_z


!WRITE(*,*) r2(1), r2(2), r2(3)
distance_x = r1(1)-r2(1)
distance_y = r1(2)-r2(2)
distance_z = r1(3)-r2(3)
bondl = sqrt(distance_x**2 + distance_y**2 + distance_z**2)/10.0
END FUNCTION bondl
