
program mc
	!use randommodule
    use inputmodule
    use initialconfigmodule
    use energypotmodule
    use trailchangemodule
    use radiuswritemodule


	Implicit none
	integer seed,N,cycles
    double precision, allocatable, dimension(:) :: rx,ry,rz,energy,vi
    integer :: v(6)  !v(1) is valance electron, v(6) is core electron
    integer a,t,i,j,p,q
    double precision cumr,radius,c,avgE,length !variables for random number for choosing random atom


    call input(N,length,seed,cycles,rx,ry,rz)



    allocate (rx(N),ry(N),rz(N),energy(cycles),vi(6))

    !call random(seed,N,rdmx,rdmy,rdmz)  !generating rdm numbers for x,y,z direction
    call initialconfig(N,length,rx,ry,rz) !giving inital random cooridinates
    call energypot(N,rx,ry,rz,1,energy)   !for 1st cycle cal LJ pot.

    !print*,rx(301)*3.4,ry(301)*3.4,rz(301)*3.4
    !print*,energy(1)
    !stop
    
    open(unit=1,file="MC.xyz")
    write(1,*) "18"
   ! do j=1,N
   !     write(1,*) "Ar",rx(j),ry(j),rz(j)
   ! end do
    

    open(unit=3,file="graph.txt")
    open(unit=4,file="radius.txt")
    open(unit=5,file="all_radius.txt")




    write(3,*) "1",energy(1)
    !print*,energy(1)
    !stop
    !open(unit=8,file="MC1.xyz")     !writing xyz file
   

    c=0.0
    radius=0
    cumr=0 !cummulative of radius value

    q=10000  !results collection data, for the last q steps
    p=cycles-q   !data/properties collection variable


    do t=2,cycles


        call trailchange(N,length,t,seed,rx,ry,rz,energy) !addjusting position and E @ same time

        !print*,"t=",t,"E=",energy(t)
     
        if (t>p) then
            c=c+energy(t)
            write(3,*) t,energy(t)
        
            !write(8,*) "512"
            !write(8,*) "Argon 512 atoms","step:", t
       
            do j=1,N
                write(1,*) "Ar",rx(j),ry(j),rz(j)
            end do


       
        
            !!!! For calculating Radius  !!!!!
            call radiuswrite(rx,ry,rz,radius,v) !v(1) is valance electron, v(6) is core electron
            write(4,*) radius
            cumr=cumr+radius  !cummulative of radius value

        end if



    end do

    avgE=c/(t*1.0)
   

    radius=cumr/(q*1.0)  !average of 10000 values

    write(3,*) "AverageE", avgE
    write(4,*) "Average radius   ",radius


    close(unit=1)
    close(unit=3)
    close(unit=4)
    close(unit=5)

end program mc

				
		
