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=17,32		!i particle and calculating for this only
    do j=1,16		!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=6.0*((1.0/(rs**(a)))-(8.0/(rs**(b))))  !original O2 interaction

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

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


!f=4.0*((0.7/(rs**(3.5)))-(8.0/(rs**(2.5))))  !!! two electrons velocity is high but resonably good)

!f=8.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.0))))   !!! works well
!f=20.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.40))))
!f=8.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.6))))
!f=8.0*((1.65/(rs**(1.6)))-(5.0/(rs**(1.6))))

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

                    !ij interaction acc. calculation

!!!! making short-range interaction zero
if (rs<=50) then
f=0
end if



                    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


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!adding extra P_E short range strong interaction for localising electrons to each nuclues
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

do i=17,32        !i particle and calculating for this only
    do j=1,16        !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


!original O2 interaction
f=40.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.0))))   !!! works well for 8 *(), 16*()

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




!f=20.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.40))))
!f=4.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.6))))
!f=8.0*((1.65/(rs**(1.6)))-(5.0/(rs**(1.6))))
!f=4.0*((0.7/(rs**(3.5)))-(8.0/(rs**(2.5))))  !!! two electrons velocity is high but resonably good)

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

!ij interaction acc. calculation


!!! making long-range interaction zero beyond 2 !!!
if (rs<=30) then
f=0
end if
if (rs>=50) then
f=0
end if



                    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



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!adding extra P_E short range strong interaction for localising electrons to each nuclues
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

do i=17,32        !i particle and calculating for this only
    do j=1,16        !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


!original O2 interaction
f=80.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.0))))   !!! works well for 8 *(), 16*()

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




!f=20.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.40))))
!f=4.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.6))))
!f=8.0*((1.65/(rs**(1.6)))-(5.0/(rs**(1.6))))
!f=4.0*((0.7/(rs**(3.5)))-(8.0/(rs**(2.5))))  !!! two electrons velocity is high but resonably good)

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

!ij interaction acc. calculation


!!! making long-range interaction zero beyond 2 !!!
if (rs>=30) then
f=0
end if



                    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






!!! P-P for nuclues 1

    do i=17,23		!i particle and calculating for this only
        do j=i+1,24		!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=12.0*((1.65/(rs**(1.4)))-(35.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


!!! making long-range interaction zero beyond 2 !!!
if (rs>=2) then
f=0
end if




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

!!! P-P for nuclues 2

    do i=25,31        !i particle and calculating for this only
        do j=i+1,32        !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=12.0*((1.65/(rs**(1.4)))-(35.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


!!! making long-range interaction zero beyond 2 !!!
if (rs>=2) then
f=0
end if




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




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!! p-p same as like hydrogen molecule  !!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   do i=17,24        !i particle and calculating for this only
        do j=25,32        !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



!!! works well, old equation !!!
!f=4.0*((2.5**38/(rs**(13)))-(8.0**2.9/(rs**(1.8))))

f=4.0*((2.87**57.0/(rs**(13.0)))-1.0*(8.0**4.8/(rs**(1.8))))

!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

!!! making short-range interaction zero beyond 2 !!!
if (rs<=100) then
f=30            !!! can adjust based on O-O bond length
end if

if (rs>=130) then
f=0
end if


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


!!!!!!!end of p-p hydrogen molecule interaction!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!! N-N same as p-p same as like hydrogen molecule!!!!
!!!!!!!!!! N-N for nuclues 1 to nuclues 2 !!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   do i=33,40        !i particle and calculating for this only
        do j=41,48        !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



!!! works well, old equation !!!
!f=4.0*((2.5**38/(rs**(13)))-(8.0**2.9/(rs**(1.8))))

f=4.0*((2.87**57.0/(rs**(13.0)))-1.0*(8.0**4.8/(rs**(1.8))))



!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

!!! making short-range interaction zero beyond 2 !!!
if (rs<=100) then
f=30
end if

if (rs>=130) then
f=0
end if


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


!!!!!!!end of p-p hydrogen molecule interaction!!!!!!





!!!!! P-N for nucleus 1

    do i=17,24		!i particle and calculating for this only
        do j=33,40		!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)))-(35.0/(rs**(1.0))))
f=12.0*((1.65/(rs**(1.4)))-(35.0/(rs**(1.0))))
!f=4.0*((1.65/(rs**(1.4)))-(15.0/(rs**(1.0)))) !(works reasonably,now works well mainly on initial coordinates shooting of electrons leads to long rotation)
!f=0
!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



!!!!! P-N for nucleus 2

    do i=25,32        !i particle and calculating for this only
        do j=41,48        !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)))-(35.0/(rs**(1.0))))
f=12.0*((1.65/(rs**(1.4)))-(35.0/(rs**(1.0))))
!f=4.0*((1.65/(rs**(1.4)))-(15.0/(rs**(1.0)))) !(works reasonably,now works well mainly on initial coordinates shooting of electrons leads to long rotation)
!f=0
!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



!!! Ne-Ne for nucleus 1

    do i=33,39		!i particle and calculating for this only
        do j=i+1,40		!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=12.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.0))))
!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=0
!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



!!! Ne-Ne for nucleus 2

    do i=41,47        !i particle and calculating for this only
        do j=i+1,48        !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=12.0*((1.65/(rs**(1.4)))-(5.0/(rs**(1.0))))
!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=0
!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,15		!i particle and calculating for this only
        do j=i+1,16		!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=0.5*((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
