

module trailchangemodule
contains
subroutine trailchange(N,length,t,seed,rx,ry,rz,energy)

use energypotmodule
implicit none

integer, intent(IN) :: N,seed,t   	!global parameters
double precision, intent(IN) :: length
double precision, dimension(:), intent(INOUT) :: rx,ry,rz   !global parameters
double precision, dimension(:), intent(INOUT) :: energy
integer i,j,q,b,a,count1s,count2s,count2p
double precision a1,m,m1,r,r1,t1,dx,rxnew,rynew,rznew,rxold,ryold,rzold,E1,E2,temp,dE,p,halfbox
double precision, dimension(6) :: rdm
double precision, dimension(3) :: u
double precision comx,comy,comz,mp,mne,rnew,rold,z,ao,p1snew,p1sold,p2snew,p2sold,p2pnew,p2pold

u=1.01
z=6
ao=68.02109

halfbox=length/2.0
        !generating N random number for selecting particles and changes with each cycle, t

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


r1=seed*t
a1=15765
m1=265487
i=1


count1s=0
count2s=0
count2p=0

do b=1,N        !for one cycle to complete N molecules
    a=b
    !a=z(b)      !a = random particle choosen (that covers N particle in N runs
    !print*,a
    !for generating 4 random nubers for moving xx,xy,xz

    a1=83645.0*b*167*(t+23)
    m=638243.0*b*17*(t+13)
    r=seed*1.0*t/5.0*b
        do i=1,6

            t1=mod((a1*r),m)
            rdm(i)=t1/m
            r=t1
            !print*,rdm(i)

        end do



    !end of generating 5 random number in rdm(1), rdm(2), rdm(3), rdm(4), rdm(5)

    !a=INT(rdm(5)*N+1)  !choosing particle at random with random number rdmx
    !print*,a
    !dx=0.1   !how much to move to new trail position

    dx=0.05
    if (a<=6) then
        dx=50.5      !*RAND()    !how much to move to new trail position for electrons
    end if
   ! if (a==18) then
    !    dx=0.000      !how much to move to new trail position for nuclues
    !end if

if (b==1) then
energy(t)=energy(t-1)
end if


    rxold=rx(a)    !storing old values of 'a' molecule
    ryold=ry(a)
    rzold=rz(a)
    E1=energy(t)  !storing old E value for comparision

    !print*,"xxold=",xxold,"xyold=",xyold,"xzold=",xzold


        rxnew=rx(a)+(2.0*rdm(1)-1.0)*dx*RAND()
        rynew=ry(a)+(2.0*rdm(2)-1.0)*dx*RAND()
        rznew=rz(a)+(2.0*rdm(3)-1.0)*dx*RAND()
        !print*,"xxnew=",xxnew,"xynew=",xynew,"xznew=",xznew
 

    rx(a)=rxnew  !assigning new position into xx,xy,xz array for cal energypot
    ry(a)=rynew
    rz(a)=rznew

!if (t>=12) then
!    print*,"before energypot",t, energy(t)
!end if


    call energypot(N,rx,ry,rz,t,energy)
    E2=energy(t)  !storing new energy pot to E2 for comparing

!if (t>=12) then
!    print*,"after energypot",t, energy(t)
!end if

    !dE=E2-E1 !change in energy due to trail position
    !print*,"dE=",dE,"E1=",E1,"E2=",E2

        
        !!! its by default, if E2 < E1 accept the move
        
      !  if(E2<E1) then
      !      energy(t)=E2
      !      rx(a)=rxnew  !assigning new position into xx,xy,xz array for cal energypot
      !      ry(a)=rynew
      !      rz(a)=rznew
      !  end if


!if (t>=12) then
!    print*,"after accepting move",t, energy(t)
!end if


        !!! for nucleus trailchange based on energy minimization !!!
    
        if(E2>E1) then
   

!if (t>=12) then
!    print*,"after rejecting move",t, energy(t)
!end if
   
            !for nuclues (a=7,18), if E2>E1 we reject without probability
            if (a>6) then

                energy(t)=E1
                rx(a)=rxold  !assigning old position into xx,xy,xz array for cal energypot
                ry(a)=ryold
                rz(a)=rzold

            end if
    
        end if



!!! only for electrons trail change !!!
        

            !accepting rejected move based on probability is only for electrons (which have r)
            if (a<=6) then
                
!Center of Mass
comx=(rx(7)*mp+rx(8)*mp+rx(9)*mp+rx(10)*mp+rx(11)*mp+rx(12)*mp+rx(13)*mne+rx(14)*mne+rx(15)*mne+rx(16)*mne+rx(17)*mne+rx(18)*mne)/&
        (6.0*mp+6.0*mne)
comy=(ry(7)*mp+ry(8)*mp+ry(9)*mp+ry(10)*mp+ry(11)*mp+ry(12)*mp+ry(13)*mne+ry(14)*mne+ry(15)*mne+ry(16)*mne+ry(17)*mne+ry(18)*mne)/&
        (6.0*mp+6.0*mne)
comz=(rz(7)*mp+rz(8)*mp+rz(9)*mp+rz(10)*mp+rz(11)*mp+rz(12)*mp+rz(13)*mne+rz(14)*mne+rz(15)*mne+rz(16)*mne+rz(17)*mne+rz(18)*mne)/&
        (6.0*mp+6.0*mne)

                !distance between nucleus and each electron
                rnew=sqrt((comx-rx(a))**2+(comy-ry(a))**2+(comz-rz(a))**2)  !radius for electron new position
                rold=sqrt((comx-rxold)**2+(comy-ryold)**2+(comz-rzold)**2)  !radius for electron old position

                
                p1snew=(4*z**3*rnew**2)*(exp(-2.0*z*rnew/ao))/(ao**3)
                p1sold=(4*z**3*rold**2)*(exp(-2.0*z*rold/ao))/(ao**3)

                p2snew=(z**3*rnew**2)*(2.0-(z*rnew)/ao)**2*exp(-1.0*z*rnew/ao)/(8*ao)
                p2sold=(z**3*rold**2)*(2.0-(z*rold)/ao)**2*exp(-1.0*z*rold/ao)/(8*ao)

                p2pnew=(z**5*rnew**4)*(exp(-1.0*z*rnew/ao))/(8*ao**5)
                p2pold=(z**5*rold**4)*(exp(-1.0*z*rold/ao))/(8*ao**5)



                u(1)=exp(-1*(E2-E1)/534)*p1snew/p1sold
                u(2)=exp(-1*(E2-E1)/534)*p2snew/p2sold
                u(3)=exp(-1*(E2-E1)/534)*p2pnew/p2pold

                !u=p1snew/p1sold

print*,"1s orbital, probability(u(1))=", u(1),"E2=",E2,"E1=",E1
print*,"2s orbital, probability(u(2))=", u(2)
print*,"2p orbital, probability(u(3))=", u(3)

                do i=1,4

                    if (i==1) then

                        if (count1s<2) then
                            if(u(1)>1) then
                                count1s=count1s+1
                                if (u(1)>1) exit

                            end if

                            if (u(1)<1) then
                                if(rdm(4)<u(1)) then
                                    count1s=count1s+1
                                    if (rdm(4)<u(1)) exit
                                end if

                            end if
                        end if

                    end if



                    if (i==2) then

                        if (count2s<2) then
                            if (u(2)>1) then
                                count2s=count2s+1
                                if (u(2)>1) exit
                            end if

                            if (u(2)<1) then
                                if (rdm(5)<u(1)) then
                                    count2s=count2s+1
                                    if (rdm(5)<u(1)) exit
                                end if

                            end if
                        end if

                    end if



                    if (i==3) then
                      
                        if (count2p<2) then
                            if (u(3)>1) then
                                count2p=count2p+1
                                if (u(3)>1) exit
                            end if

                            if (u(3)<1) then
                                if (rdm(6)<u(1)) then
                                    count2p=count2p+1
                                    if (rdm(6)<u(1)) exit
                                end if

                            end if
                        end if

                    end if


                    if (i==4) then  !!!if all three cases trail move failed, we assign old values
        
                        energy(t)=E1
                        rx(a)=rxold  !assigning old position into xx,xy,xz array for cal energypot
                        ry(a)=ryold
                        rz(a)=rzold

                    end if

                end do


        end if


end do

    !print*,"u=",u,"rdm=",rdm(4)
end subroutine trailchange
end module trailchangemodule
	  
