module accmodule
contains
subroutine acc(Elec,Proton,N,Ne,t,length,rx,ry,rz,ax,ay,az,LJ,a,b)   !calculate accelration & LJ @ t
	
	implicit none

	integer, intent(IN) :: Elec,Proton,N,Ne,t
	double precision, dimension(:), intent(inout) :: rx,ry,rz
	double precision, dimension(:), intent(inout) :: ax,ay,az
	double precision, dimension(:), intent(out) :: LJ
    double precision, intent(IN) :: length,a,b
	double precision rs,f,KJ,halfbox,dx,dy,dz,p1,p2,p3,p,mp,mne
	integer i,j
	


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


    halfbox=length/2.0

    KJ=0.0
	f=0.0
	ax=0.0
	ay=0.0
	az=0.0
   ! print*,"ax(50)=",ax(100)

    p=0

!!! e = 1,3 
!!! p = 2,4
!!! neutron = 5,6

!!! i should be proton, j should be electorn


do i=93,184		!i particle and calculating for this only
    do j=1,92		!j particle pairing with i for pairs (i,j)

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


            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


                    f=4.0*((1.0/(rs**(a)))-(8.0/(rs**(b))))  !!! two electrons velocity is high but resonably good)

                    !f=1.0*((1.0/(rs**(a)))-(8.0/(rs**(b))))

                    !ij interaction acc. calculation



                    p1=(dx*f*1.0)
                    p2=(dy*f*1.0)
                    p3=(dz*f*1.0)
                    !print*,"ji",p1,p2,p3


                    !i=1, proton....ax=p1/(mp*)

                    ax(i)=ax(i)+p1*1.0/(mp*1.0)  !acc. keeps on adding upon adding force for each pair (i,j)
                    ay(i)=ay(i)+p2*1.0/(mp*1.0)
                    az(i)=az(i)+p3*1.0/(mp*1.0)

                    !ji interaction acc. calculation

                    !the same force f also acts on j atom and so adding to is acc. since we dont repeat same pair for LJ cal. adding acc. to j atom now itself now "dx=-dx"


                    !j=2, electron.....ax=

                    ax(j)=ax(j)-p1*1.0  !acc. keeps on adding upon adding force for each pair (j,i)
                    ay(j)=ay(j)-p2*1.0
                    az(j)=az(j)-p3*1.0


!!! end of (2,1) pair !!!


         end do
     end do



    do i=93,183		!i particle and calculating for this only
        do j=i+1,184		!j particle pairing with i for pairs (i,j)



!!! pair (p-p) !!!
!i=proton
!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




f=4.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.0)))) !(works reasonably,now works well mainly on initial coordinates shooting of electrons leads to long rotation)

!f=2.0*((1.2/(rs**(1.4)))-(5.0/(rs**(1.0))))  (bigger aatom size)

!f=((1.9/(rs**(1.4)))-(5.0/(rs**(1.0))))  (works really well) (compared to first eq, only depth changes, one has to adjust it depending on helium nuclues potential energy)

!ij interaction acc. calculation



p1=(dx*f*1.0)
p2=(dy*f*1.0)
p3=(dz*f*1.0)
!print*,"ji",p1,p2,p3


!i=1, proton....ax=p1/(mp*)

ax(i)=ax(i)+p1*1.0/(mp*1.0)  !acc. keeps on adding upon adding force for each pair (i,j)
ay(i)=ay(i)+p2*1.0/(mp*1.0)
az(i)=az(i)+p3*1.0/(mp*1.0)

!ji interaction acc. calculation

!the same force f also acts on j atom and so adding to is acc. since we dont repeat same pair for LJ cal. adding acc. to j atom now itself now "dx=-dx"


!j=2, electron.....ax=

ax(j)=ax(j)-p1*1.0/(mp*1.0)  !acc. keeps on adding upon adding force for each pair (j,i)
ay(j)=ay(j)-p2*1.0/(mp*1.0)
az(j)=az(j)-p3*1.0/(mp*1.0)


!!! end of (2,4) pair !!!



        end do
    end do



    do i=93,184		!i particle and calculating for this only
        do j=185,330		!j particle pairing with i for pairs (i,j)



!!! pair (p-ne) !!!
!i=proton
!j=neutron


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


!only for oxygen atom
!f=4.0*((0.9/(rs**(3.0)))-(15.0/(rs**(2.0)))) !(works reasonably,now works well mainly on initial coordinates shooting of electrons leads to long rotation)


f=4.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.0)))) !(works reasonably,now works well mainly on initial coordinates shooting of electrons leads to long rotation)

!f=2.0*((1.2/(rs**(1.4)))-(5.0/(rs**(1.0))))  (bigger aatom size)

!f=((1.9/(rs**(1.4)))-(5.0/(rs**(1.0))))  (works really well) (compared to first eq, only depth changes, one has to adjust it depending on helium nuclues potential energy)

!ij interaction acc. calculation



p1=(dx*f*1.0)
p2=(dy*f*1.0)
p3=(dz*f*1.0)
!print*,"ji",p1,p2,p3


!i=1, proton....ax=p1/(mp*)

ax(i)=ax(i)+p1*1.0/(mp*1.0)  !acc. keeps on adding upon adding force for each pair (i,j)
ay(i)=ay(i)+p2*1.0/(mp*1.0)
az(i)=az(i)+p3*1.0/(mp*1.0)

!ji interaction acc. calculation

!the same force f also acts on j atom and so adding to is acc. since we dont repeat same pair for LJ cal. adding acc. to j atom now itself now "dx=-dx"


!j=2, electron.....ax=

ax(j)=ax(j)-p1*1.0/(mne*1.0)  !acc. keeps on adding upon adding force for each pair (j,i)
ay(j)=ay(j)-p2*1.0/(mne*1.0)
az(j)=az(j)-p3*1.0/(mne*1.0)


!!! end of (2,5) pair !!!


        end do
    end do



    do i=185,329		!i particle and calculating for this only
        do j=i+1,330		!j particle pairing with i for pairs (i,j)



!!! pair (ne-ne) !!!
!i=neutron
!j=neutron


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



!only for oxygen atom
!f=4.0*((0.8/(rs**(6.0)))-(25.0/(rs**(0.5)))) !(works reasonably,now works well mainly on initial coordinates shooting of electrons leads to long rotation)

f=4.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.0)))) !(works reasonably,now works well mainly on initial coordinates shooting of electrons leads to long rotation)

!f=2.0*((1.2/(rs**(1.4)))-(5.0/(rs**(1.0))))  !(bigger aatom size)

!f=((1.9/(rs**(1.4)))-(5.0/(rs**(1.0))))  (works really well) (compared to first eq, only depth changes, one has to adjust it depending on helium nuclues potential energy)

!ij interaction acc. calculation



p1=(dx*f*1.0)
p2=(dy*f*1.0)
p3=(dz*f*1.0)
!print*,"ji",p1,p2,p3


!i=1, proton....ax=p1/(mp*)

ax(i)=ax(i)+p1*1.0/(mne*1.0)  !acc. keeps on adding upon adding force for each pair (i,j)
ay(i)=ay(i)+p2*1.0/(mne*1.0)
az(i)=az(i)+p3*1.0/(mne*1.0)

!ji interaction acc. calculation

!the same force f also acts on j atom and so adding to is acc. since we dont repeat same pair for LJ cal. adding acc. to j atom now itself now "dx=-dx"


!j=2, electron.....ax=

ax(j)=ax(j)-p1*1.0/(mne*1.0)  !acc. keeps on adding upon adding force for each pair (j,i)
ay(j)=ay(j)-p2*1.0/(mne*1.0)
az(j)=az(j)-p3*1.0/(mne*1.0)


!!! end of (5,6) pair !!!


        end do
    end do











    do i=1,91		!i particle and calculating for this only
        do j=i+1,92		!j particle pairing with i for pairs (i,j)



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


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


!f=4.0*((1.0/(rs**(3.5)))-(8.0/(rs**(1.0))))

!f=4.0*((2.10/(rs**(6)))-(8.0/(rs**(2.0))))

f=4.0*((1.0/(rs**(1.2)))-(1.5/(rs**(1.0))))

!f=4.0*((1.0/(rs**(1.2)))-(1.2/(rs**(1.0))))  !!!original one used in all cases  !!! for carbon case, two electrons move little bit away, but reasonably good !!!

!ij interaction acc. calculation



p1=(dx*f*1.0)
p2=(dy*f*1.0)
p3=(dz*f*1.0)
!print*,"ji",p1,p2,p3


!i=1, proton....ax=p1/(mp*)

ax(i)=ax(i)+p1*1.0  !acc. keeps on adding upon adding force for each pair (i,j)
ay(i)=ay(i)+p2*1.0
az(i)=az(i)+p3*1.0

!ji interaction acc. calculation

!the same force f also acts on j atom and so adding to is acc. since we dont repeat same pair for LJ cal. adding acc. to j atom now itself now "dx=-dx"


!j=2, electron.....ax=

ax(j)=ax(j)-p1*1.0  !acc. keeps on adding upon adding force for each pair (j,i)
ay(j)=ay(j)-p2*1.0
az(j)=az(j)-p3*1.0


!!! end of (1,3) pair !!!





         end do
     end do


	LJ(t)=KJ	!for calculating LJ potential used for PE at time t




end subroutine acc

end module accmodule
