
program md
    use inputmodule
	use initialconfignewmodule
	use accmodule
	!use averagepropmodule
    !use velrescalemodule
	use velverletmodule
    use radiuswritemodule
    use valenceenergymodule
    use nucleusradiusmodule
    use transformmodule
    use volumeslicemodule
    use electrondistancemodule
    use velocitymodule
    use atomtotalenergymodule



	Implicit none
	integer :: Elec,Proton,N,Ne,t,itr,seed,i,j,h,p,q,rescale,total!,v1,v2
    !integer :: v3,v4,v5, v6
    integer :: v(12)  !v(1) is valance electron, v(6) is core electron
    double precision dt,length,Temp,c,a,b,radius,cumr,cumnucr,nucradius,atomEtotal,atomPotential,atomKE
        !,val1Etotal,cumval1Etotal
    !double precision val2Etotal,cumval2Etotal,nucradius
    !double precision val3Etotal,cumval3Etotal
    !double precision val4Etotal,cumval4Etotal
    !double precision val5Etotal,cumval5Etotal
    !double precision val6Etotal,cumval6Etotal
    double precision :: valEtotal(12),cumvalEtotal(12),valPot(12),cumvalPot(12),valKE(12),cumvalKE(12)
	double precision, allocatable, dimension(:) :: rx,ry,rz,vx,vy,vz,ax,ay,az
	double precision, allocatable, dimension(:) :: Potential,Etotal,caltemp,KE,LJ   !average temperature at each instant time




    call input(Elec,Proton,N,Ne,length,dt,seed,Temp,itr,a,b,q)



	allocate (rx(N),ry(N),rz(N),vx(N),vy(N),vz(N),ax(N),ay(N),az(N))

	allocate (Caltemp(itr),LJ(itr),Potential(12),KE(12),Etotal(12))



	call initialconfignew(N,length,seed,rx,ry,rz,vx,vy,vz)  !that will give initial value of x,v at time t = 1


    open(unit=6,file="newmd.xyz")


    write(6,*) "36"
    write(6,*) "particle dynamics p & e"

    !!!!!    Writing coordinates !!!!!!!

 !   do j=1,N
 !      write(6,*) "Ar",rx(j),ry(j),rz(j)
 !   end do


    !!!!!    End of Writing coordinates !!!!!!!

    call acc(Elec,Proton,N,Ne,1,length,rx,ry,rz,ax,ay,az,LJ,a,b)  !to call acc. at initial value t=1

    call atomtotalenergy(1,rx,ry,rz,vx,vy,vz,atomEtotal,atomPotential,atomKE) !calculating atom's energies at t=1


    i=2

    open(unit=7,file="radius.txt")
    open(unit=8,file="Etotal.txt")
    open(unit=9,file="Potential.txt")
    open(unit=10,file="KE.txt")
    open(unit=11,file="nucradius.txt")
    !open(unit=10,file="sorting.txt")
    open(unit=12,file="volumeslice.txt")
    open(unit=13,file="electrondistances.txt")
    open(unit=14,file="velocity.txt")
    open(unit=15,file="atomenergy_all.txt")
    open(unit=16,file="all_radius.txt")

!files for Phase shift observations (TE,PE,KE)
open(unit=28,file="atomenergy_0.txt")
open(unit=17,file="atomenergy_1.txt")
open(unit=18,file="atomenergy_2.txt")
open(unit=19,file="atomenergy_3.txt")
open(unit=20,file="atomenergy_4.txt")
open(unit=21,file="atomenergy_5.txt")
open(unit=22,file="atomenergy_6.txt")
open(unit=23,file="atomenergy_7.txt")
open(unit=24,file="atomenergy_8.txt")
open(unit=25,file="atomenergy_9.txt")
open(unit=26,file="atomenergy_10.txt")
!end of files for Phase shift observations (TE,PE,KE)
    
    open(unit=27,file="first_PE.txt")


    write(12,*) "36"
    write(12,*) "particle dynamics p & e"

    write(15,*) "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE


    radius=0
    !r=0
    cumr=0 !cummulative of radius value
    cumnucr=0  !cummulative nucleus radius

    do j=1,9
        valEtotal(j)=0
        cumvalEtotal(j)=0
        valPot(j)=0
        cumvalPot(j)=0
        valKE(j)=0
        cumvalKE(j)=0
    end do


    p=itr-q   !data/properties collection
    h=itr-10000  !coordinates file data collection for last 10,000 time steps

    do i=2,itr

           call velverlet(Elec,Proton,N,Ne,dt,i,length,rx,ry,rz,vx,vy,vz,ax,ay,az,LJ,a,b)   !will give new values for x,v at t

            !first PE data collection
        
            if (i<10000) then

                call atomtotalenergy(i,rx,ry,rz,vx,vy,vz,atomEtotal,atomPotential,atomKE) !calculating atom's energies at i timestep
        
    
                write(27,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE

            end if
            !end of first PE data collection


            !files for Phase shift observations (TE,PE,KE)

            !if (i>500000) then

                call atomtotalenergy(i,rx,ry,rz,vx,vy,vz,atomEtotal,atomPotential,atomKE) !calculating atom's energies at i timestep
                
                !all timestep files
                !write(15,*) "time step: ",i, "AtomEtotal     ",atomEtotal,"   AtomPotential   ",atomPotential,"    AtomKE    ",atomKE

            !end if

            !end of files for Phase shift observations (TE,PE,KE)


            !call transform(rx,ry,rz)  !!!transform nucleus COM to origin


            if (i>h) then

                !!!!!    Writing coordinates !!!!!!!

                !write(6,*) "C",rx(1),ry(1),rz(1)
                !write(6,*) "H",rx(2),ry(2),rz(2)
                !write(6,*) "N",rx(3),ry(3),rz(3)
                !write(6,*) "O",rx(4),ry(4),rz(4)
                !write(6,*) "P",rx(5),ry(5),rz(5)
                !write(6,*) "He",rx(6),ry(6),rz(6)
               
                do j=1,N
                    write(6,*) "Ar",rx(j),ry(j),rz(j)
                end do

                !!!!!    End of Writing coordinates !!!!!!!
        
 
    
            end if

            !!! calculating cummulative radius !!!


            !!! end of calculating cummulative radius !!!


	end do


    write(7,*) "Average radius   ",radius
    write(8,*) "AvgEtotal from val to core  ",valEtotal
    write(9,*) "AvgPot from val to core  ",valPot
    write(10,*) "AvgKE from val to core  ",valKE
    write(11,*) "Average nucleus radius   ",nucradius

    close(unit=6)
    close(unit=7)
    close(unit=8)
    close(unit=9)
    close(unit=10)
    close(unit=11)
    close(unit=12)
    close(unit=13)
    close(unit=14)
    close(unit=15)
    close(unit=16)
    close(unit=17)
    close(unit=18)
    close(unit=19)
    close(unit=20)
    close(unit=21)
    close(unit=22)
    close(unit=23)
    close(unit=24)
    close(unit=25)
    close(unit=26)
    close(unit=27)
    close(unit=28)


end program md


		
