module velverletmodule
contains
subroutine velverlet(Elec,Proton,N,Ne,dt,t,length,rx,ry,rz,vx,vy,vz,ax,ay,az,LJ,a,b)
	
	use accmodule
	implicit none	
	integer, intent(in) :: Elec,Proton,N,Ne,t
	double precision, intent(in) :: dt,length,a,b
    double precision, dimension(:), intent(inout) :: ax,ay,az
    double precision, dimension(:), intent(inout) :: rx,ry,rz
    double precision, dimension(:), intent(out) :: vx,vy,vz
    double precision, dimension(:), intent(inout) :: LJ
    double precision, dimension(N) :: axold,ayold,azold,rxold,ryold,rzold,vxold,vyold,vzold
    integer i,j
    double precision halfbox


    halfbox=length/2.0

    !a=-1.0*L  !for periodic BC

	rxold=rx	!preserving xx(t-1) value in xxold to use in VV formula
	ryold=ry
	rzold=rz
!print*,"xxold(100)",xxold(100)
    !open(unit=2,file="check.txt")

    do i=1,N
        !do
                rx(i)=rxold(i)+(dt*vx(i))+(0.5*dt*dt*ax(i))  
                ry(i)=ryold(i)+(dt*vy(i))+(0.5*dt*dt*ay(i))
                rz(i)=rzold(i)+(dt*vz(i))+(0.5*dt*dt*az(i))

    end do

        !updated new x values at t and now cal. velocity at t

	   vxold=vx    !preserving v(t-1) value to vxold to use in VV formula
	   vyold=vy
	   vzold=vz

	   axold=ax    !preserving a(t-1) value to axold to use in VV formula
	   ayold=ay
	   azold=az
		
	   call acc(Elec,Proton,N,Ne,t,length,rx,ry,rz,ax,ay,az,LJ,a,b)  !new a(t) will be stored in ax,ay,az

	do i=1,N
	      
	   vx(i)=(vxold(i)+(0.5*dt*(axold(i)+ax(i))))		
	   vy(i)=(vyold(i)+(0.5*dt*(ayold(i)+ay(i))))		
	   vz(i)=(vzold(i)+(0.5*dt*(azold(i)+az(i))))

        !vx(1)=0.0
        !vy(1)=0.0
        !vz(1)=0.0
		!print*,vx(2),vy(2),vy(2)

	end do

end subroutine velverlet
end module velverletmodule	      	
	

