module initialconfignewmodule
contains
subroutine initialconfignew(N,length,seed,rx,ry,rz,vx,vy,vz)

    !use velrescalemodule
	implicit none
	integer, intent(IN) :: N,seed 	!global parameters
	double precision, intent(IN) :: length
	double precision, dimension(:), intent(OUT) :: rx,ry,rz,vx,vy,vz   !global parameters
	double precision q,halfbox,v,r,s                  !,Vmax=0.2	!local parameters
	integer :: i,j,k,g,p,b,m,ab 
    double precision, dimension(4) :: x1,y1,z1

    halfbox=length/2.0             !for shifiting axis to (h,h,h), @ centre of cube
	p=(N**(1./3.))		!fixing number of particles to place in each direction
	!q=L*1.0/p			!cal. equal spacing for p particles in L space
	g=1



    !simple placing of atoms in a cube edges with a space, starts from 0.5 units
    !m=5
    !x1=(/0.75,0.25,0.25,0.75/)
    !y1=(/0.75,0.25,0.75,0.25/)
    !z1=(/0.75,0.75,0.25,0.25/)
    !q=length/(m*1.0)
!print*,"q=",q

	!do i=1,m
	!	do j=1,m
	!		do k=1,m
     !           do ab=1,4
!
 !                   if (g<=N) then
  !                      rx(g)=(i-x1(ab))*q       !placing starts from 0.5
   !                     ry(g)=(j-y1(ab))*q       !for y direction
    !                    rz(g)=(k-z1(ab))*q       !for z direction
     !                   g=g+1
      !              end if

       !         end do  !counting no of particles placed so far
		!	end do
		!end do
!	end do

!shifting the coordinates and making center of cube(h,h,h) as origin

  !  do i=1,N
   !     rx(i)=rx(i)-halfbox
    !    ry(i)=ry(i)-halfbox
     !   rz(i)=rz(i)-halfbox
        !write(8,*)"Ar",xx(i),xy(i),xz(i)
   ! end do




!!assigning initial coordinates to particles

!!! electron coordinates
s=2.5
r=3.0
v=4.5

!electrons for first nuclues
    do i=1,8


        !rx(i)=10.5*length*RAND()      !placing at centre of box
        !ry(i)=20.5*length*RAND()
        !rz(i)=5.5*length*RAND()

rx(i)=5.5+15.0*RAND()      !placing at centre of box
ry(i)=5.5+15.0*RAND()
rz(i)=5.5+15.0*RAND()

        s=s+0.4
        r=r+0.3
        v=v-0.4

    end do

!electrons for second nuclues
    do i=9,16

        !rx(i)=10.5*length*RAND()      !placing at centre of box
        !ry(i)=20.5*length*RAND()
        !rz(i)=5.5*length*RAND()

rx(i)=75.5+10.0*RAND()      !placing at centre of box
ry(i)=75.5+10.0*RAND()
rz(i)=75.5+10.0*RAND()

        s=s+0.4
        r=r+0.3
        v=v-0.4

    end do


!!! proton and neutron coordinates for oxygen atom 1

    do i=17,24

        rx(i)=5.0+0.1*RAND()      !placing at centre of box
        ry(i)=5.0+0.1*RAND()
        rz(i)=5.0+0.01*RAND()

    end do

    do i=33,40

        rx(i)=5.0+0.1*RAND()      !placing at centre of box
        ry(i)=5.0+0.1*RAND()
        rz(i)=5.0+0.01*RAND()

    end do

!!! proton and neutron coordinates for oxygen atom 2

    do i=25,32

        rx(i)=75.0+0.1*RAND()      !placing at centre of box
        ry(i)=75.0+0.1*RAND()
        rz(i)=75.0+0.01*RAND()

    end do

    do i=41,48

        rx(i)=75.0+0.1*RAND()      !placing at centre of box
        ry(i)=75.0+0.1*RAND()
        rz(i)=75.0+0.01*RAND()

    end do




	!giving initial velocities acc. to boltzmann normal distribution with scaling temp T
	b=seed
	do i=1,N
        b=seed
		call boltzmann(i,b,v)   !i,b serves as a seed
        vx(i)=v
        b=b+8223
		call boltzmann(i,b,v)
        vy(i)=v
        b=b+9334321
        call boltzmann(i,b,v)
        vz(i)=v

	end do


!!! velocity of proton n neutron to zero
    do i=17,48
        vx(i)=0.0
        vy(i)=0.0
        vz(i)=0.0
    end do

!!! end of velocity to zero



 !call velrescale(N,vx,vy,vz,Temp)



end subroutine initialconfignew




    subroutine boltzmann(i,b,v)
        implicit none
        integer, intent(IN) :: i
        integer, intent(inout) :: b
        double precision, intent(out) :: v
        double precision a,x,y
        integer u,w

        u=i+b
        w=i+b
		do

			x=rand((i+207*i)*u)
			y=rand((i+223*i)*w)
			x=2.0*x-1.0                 !interval of -1 to 1
			y=2.0*y-1.0                 !interval of -1 to 1
			a=(x*x)+(y*y)       !a should be in between 0 to 1
			if(a>0.0.AND.a<1.0)exit
            u=int(x*10000)
            w=int(y*10000)

		end do
		v=x*sqrt(-2.0*log(a)/a)
	end subroutine boltzmann
		

end module initialconfignewmodule
	  
