module atomtotalenergymodule
contains
subroutine atomtotalenergy(t,rx,ry,rz,vx,vy,vz,atomEtotal,atomPotential,atomKE)

	implicit none
    integer, intent(in) :: t
    double precision, dimension(:), intent(in) :: rx,ry,rz,vx,vy,vz
    double precision, intent(inout) :: atomEtotal,atomPotential,atomKE
    integer p1,p2,e1,e2,n1,n2,i,j
    double precision mp,mne,V,KE,dx,dy,dz,rs 


mp=1836.1526    !proton/electron mass ratio
mne=1838.684    !neutron/electron mass ratio

p1=10
p2=18
e1=1
e2=9
n1=19
n2=28


atomEtotal=0.0


!!!!! Calculating Potential Energy  !!!!!!!!

atomPotential=0.0

!Electron-Proton interaction

do i=e1,e2        !i particle and calculating for this only
    
    do j=p1,p2        !j particle pairing with i for pairs (i,j)

!!! pair (p-e) !!!
!i=electron
!j=proton

            dx=rx(i)-rx(j)
            dy=ry(i)-ry(j)
            dz=rz(i)-rz(j)



            rs=(dx*dx)+(dy*dy)+(dz*dz)     !new rs will consider MIC rs's also that are qualifed for cutoff filtrations
            rs=sqrt(rs)
            !p=rs

            V=4.0*((1.0/(1.5*(rs**1.5)))+(8.0*rs))
      
            atomPotential=atomPotential+V

    end do
end do


!Electron-Electron interaction

do i=e1,e2-1        !i particle and calculating for this only
    do j=i+1,e2        !j particle pairing with i for pairs (i,j)

!!! pair (e-e) !!!
!i=electron
!j=electron

        !if (i/=j) then

            dx=rx(i)-rx(j)
            dy=ry(i)-ry(j)
            dz=rz(i)-rz(j)



            rs=(dx*dx)+(dy*dy)+(dz*dz)     !new rs will consider MIC rs's also that are qualifed for cutoff filtrations
            rs=sqrt(rs)
            !p=rs

            V=4.0*((-1*(rs**0.8)/0.8)+(1.5*rs))
      
            atomPotential=atomPotential+V

        !end if

    end do
end do



!proton-proton, neutron-neutron, proton-neutron interactions

do i=p1,n2-1        !i particle and calculating for this only
    do j=i+1,n2        !j particle pairing with i for pairs (i,j)

!!! pair (e-e) !!!
!i=electron
!j=electron

        !if (i/=j) then

            dx=rx(i)-rx(j)
            dy=ry(i)-ry(j)
            dz=rz(i)-rz(j)



            rs=(dx*dx)+(dy*dy)+(dz*dz)     !new rs will consider MIC rs's also that are qualifed for cutoff filtrations
            rs=sqrt(rs)
            !p=rs

            V=4.0*((-2.75*(rs**0.6))+(5.0*rs))
      
            atomPotential=atomPotential+V

        !end if

    end do
end do






!!!!! Calculating Kinetic Energy  !!!!!!!!

atomKE=0.0

!electrons

do i=e1,e2        !i particle and calculating for this only


    KE=0.5*(vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i))

    atomKE=atomKE+KE
        

end do


!protons

do i=p1,p2        !i particle and calculating for this only


    KE=0.5*mp*(vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i))

    atomKE=atomKE+KE
        

end do



!neutrons

do i=n1,n2        !i particle and calculating for this only


    KE=0.5*mne*(vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i))

    atomKE=atomKE+KE
        

end do


atomEtotal=atomPotential+atomKE







!files for Atom's Energy Quantised observations (TE,PE,KE)
!files for Atom's Energy Quantised observations (TE,PE,KE)
!files for Atom's Energy Quantised observations (TE,PE,KE)

i=t  !!

if (i>10000) then
    if (i<=500000) then

        write(28,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if



if (i>500000) then
    if (i<=1500000) then

        write(17,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>1500000) then
    if (i<=2500000) then

        write(18,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>2500000) then
    if (i<=3500000) then

        write(19,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>3500000) then
    if (i<=4500000) then
 
        write(20,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>4500000) then
    if (i<=5500000) then

        write(21,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>5500000) then
    if (i<=6500000) then

        write(22,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>6500000) then
    if (i<=7500000) then
 
        write(23,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>7500000) then
    if (i<=8500000) then
 
        write(24,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>8500000) then
    if (i<=9500000) then

        write(25,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

if (i>9500000) then
    if (i<=10000000) then

        write(26,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE
    end if
end if

!end of files for Atom's Energy Quantised observations (TE,PE,KE)
!end of files for Atom's Energy Quantised observations (TE,PE,KE)
!end of files for Atom's Energy Quantised observations (TE,PE,KE)





end subroutine atomtotalenergy
end module atomtotalenergymodule
	  
