#!bin/bash
for seed in 10987 5454 53345 909 435532 6757575 3232 456 1189 7331
do
echo "Computing Rh for traj" $seed
cat > compute_hydrodynamic.f90 << END
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_${seed}.xyz")
OPEN(UNIT=11, FILE="hydrodynamic_${seed}.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
END
gfortran compute_hydrodynamic.f90
./a.out
cat hydrodynamic_${seed}.txt >> HYDRODYNAMIC_ALL
echo "Done for seed" $seed
done
awk '{print $2}' HYDRODYNAMIC_ALL > junk
gfortran average.f95
./a.out junk
./histogram_manual junk dist.rh 0.01


