program CalculateInteractionEnergy
implicit none
real :: sigma_CO, epsilon_CO
real :: distance, lj_energy, total_interaction_energy,Fx,Fy,Fz,Dx,Dy,Dz,F_z,F_y,F_x
integer :: snapshot, num_snapshots, num_atoms,N_o, i, j,atom_id(5000),ii,k, time
character(4) :: line
character(100), parameter :: file_path = 'data1.xyz'
integer, parameter :: max_atoms = 1000
real*8 x(5000),y(5000),z(5000)
sigma_CO = 3.289 ! Angstrom
epsilon_CO = 0.09242*0.0433634 ! KCal/mol**0.0433634=eV
num_snapshots = 50
open(unit=10, file=file_path, status='old')
do k = 1, num_snapshots
read(10, *) num_atoms
read(10,*) line
N_o=0
do j = 1, num_atoms
read(10, *) ii,atom_id(j) , x(j), y(j), z(j)
if(atom_id(j).eq. 2) N_o=N_o+1
end do
do i=1,num_atoms
do j=1,num_atoms
if(atom_id(j).eq. 1 .and. atom_id(i).eq. 2) then
distance=sqrt((x(i)-x(j))**2+(y(i)-y(j))**2+(z(i)-z(j))**2)
lj_energy = lj_energy+4.0 * epsilon_CO * ((sigma_CO / distance)**12 - (sigma_CO / distance)**6)
Dx=x(i)-x(j); Dy=y(i)-y(j); Dz=z(i)-z(j)
Fx=Fx+(-6*(Dx*(sigma_CO**6 / distance**8))+12*(Dx*(sigma_CO**12 / distance**14)))*4.0 * epsilon_CO
Fy=Fy+(-6*(Dy*(sigma_CO**6 / distance**8))+12*(Dy*(sigma_CO**12 / distance**14 )))*4.0 * epsilon_CO
Fz=Fz+(-6*(Dz*(sigma_CO**6 / distance**8))+12*(Dz*(sigma_CO**12 / distance**14 )))*4.0 * epsilon_CO
endif
enddo
enddo
write(170,*)lj_energy/N_o,N_o,Fz
write(6,*)lj_energy/N_o,N_o,Fz,k
total_interaction_energy=total_interaction_energy+lj_energy/N_o
F_x=F_x+Fx !/N_o
F_y=F_y+Fy !/N_o
F_z=F_z+Fz !/N_o
lj_energy=0.0d0
Fx=0.0d0
Fy=0.0d0
Fz=0.0d0
write(120, *) N_o
write(120,*) 'atoms'
do j = 1, num_atoms
if(atom_id(j) .eq. 2) write(120, 10) atom_id(j) , x(j), y(j), z(j)
end do
end do
write(6,*) 'average='
write(6,*) total_interaction_energy/50.0, N_o !eV
write(6,*) F_x/50,F_y/50,F_z/50 !eV
close(unit=10)
10 format(I4,3E25.16)
stop
end program CalculateInteractionEnergy