FOTRAN program for Displacement of particle and Acceptance of Move using Monte Carlo simulation

Average enrgy calculation involves both Displacement of the selected particle and LJ potential calculation

// Displacement of particle and Acceptance of Move

subroutine displace(success, ihs)

Use sp
Use xyz
Use randmod
Use moves
Use tmmod

implicit none

logical            :: success
integer            :: ihs,i

integer            :: mol

real,dimension(1:5)    :: xold, yold, zold
real,dimension(1:5)    :: xnew, ynew, znew
real            :: dx, dy, dz, xo,yo,zo,xn,yn,zn
real            :: emffnew, emffold, ioneffold, ioneffnew, ioneffold1, ioneffnew1
real            :: emwfnew, emwfold, ionewfnew, ionewfold, ionewfnew1, ionewfold1
real            :: de, dlnpsi

success = .false.
emffnew = 0
emffold = 0
emwfnew = 0
emwfold = 0
ioneffold = 0
ioneffold1 = 0
ioneffnew = 0
ioneffnew1 = 0
ionewfold = 0
ionewfold1 = 0
ionewfnew = 0
ionewfnew1 = 0
pacc = 0.0

if( Nmol == 0 ) return
ihs = 1
if( ran2(Seed) > 0.33 ) ihs = 2

if( ihs == 1 ) then

        dx = 0.0
        dy = 0.0
        dz = ( 2.0 * ran2(Seed) - 1.0 ) * dh * bh

else

        dx = ( 2.0 * ran2(Seed) - 1.0 ) * ds * bs
        dy = ( 2.0 * ran2(Seed) - 1.0 ) * ds * bs
        dz = 0.0

end if

mol = int( Nmol * ran2(Seed) ) + 1

if (mol .eq. Nmol) then
             
    xo = Xr
    yo = yr
    zo = zr

    xn = xo + dx
    yn = yo + dy
    zn = zo + dz

    if( xn > bs )  xn = xn - bs * aint( xn / bs )
        if( yn > bs )  yn = yn - bs * aint( yn / bs )
        if( nwalls == 0 .and. zn > bh ) zn = zn - bh * aint( zn / bh )
      
    if( nwalls > 0  .and. zn > bh )  return
    if( xn < 0.0 )  xn = xn - bs * aint( xn / bs - 1.0 )
    if( yn < 0.0 )  yn = yn - bs * aint( yn / bs - 1.0 )
    if( nwalls == 0 .and. zn < 0.0 )  zn = zn - bh * aint( zn / bh - 1.0 )
    if( nwalls > 0  .and. zn < 0.0 )  return

        call ljmoleculeion(ioneffold,ionewfold)
        call ljmoleculeion1(ioneffold1,ionewfold1)

        ioneffold = (1-lambda) * ioneffold + (lambda * ioneffold1)

        Xr = xn
        Yr = yn
        Zr = zn

        call ljmoleculeion(ioneffnew,ionewfnew)
        call ljmoleculeion1(ioneffnew1,ionewfnew1)

        ioneffnew = (1-lambda) * ioneffnew + (lambda * ioneffnew1)

    de = (ioneffnew + ionewfnew ) - (ioneffold + ionewfold)
        
else
    xold(:) = X(mol,:)
    yold(:) = Y(mol,:)
    zold(:) = Z(mol,:)

    xnew = xold + dx
    ynew = yold + dy
    znew = zold + dz

    do i=1,5
        if( xnew(i) > bs )  xnew(i) = xnew(i) - bs * aint( xnew(i) / bs )
        if( ynew(i) > bs )  ynew(i) = ynew(i) - bs * aint( ynew(i) / bs )
        if( nwalls == 0 .and. znew(i) > bh ) znew(i) = znew(i) - bh * aint( znew(i) / bh )
        if( nwalls > 0  .and. znew(i) > bh )  return
                          
        if( xnew(i) < 0.0 )  xnew(i) = xnew(i) - bs * aint( xnew(i) / bs - 1.0 )
        if( ynew(i) < 0.0 )  ynew(i) = ynew(i) - bs * aint( ynew(i) / bs - 1.0 )
        if( nwalls == 0 .and. znew(i) < 0.0 )  znew(i) = znew(i) - bh * aint( znew(i) / bh - 1.0 )
        if( nwalls > 0  .and. znew(i) < 0.0 )  return

    end do
                      
    call ljmolecule(mol, emffold, emwfold)
      

    X(mol,:) = xnew(:)
    Y(mol,:) = ynew(:)
    Z(mol,:) = znew(:)

    call ljmolecule(mol, emffnew, emwfnew)
     
    de = ( emffnew + emwfnew ) - ( emffold + emwfold )
     
endif


dlnpsi = - beta * de

pacc = -beta * de

if ( de > 2160.0) then

 !write(*,*)' problem is in displaces stop stop stop stopstop*************************',mol
  write(*,*)' change',ihs,dx,dy,dz,de
  write(*,*) '  energies', emffnew,emffold

!  write(*,*)' xold', xold
 ! write(*,*)' yold', yold
  !write(*,*)' zold', zold


  !write(*,*)' xnew', xnew
  !write(*,*)' ynew', ynew
  !write(*,*)' znew', znew
 
 end if
if( pacc < 0.0 ) then

    pacc = exp(pacc)

else

    pacc = 1.0

end if

if( log( ran2(Seed) ) < dlnpsi ) then  
      
    success = .True.
    !ebff = ebff + emffnew - emffold + ioneffnew - ioneffold
      
    ebff=ebff+de
    ebwf = ebwf + emwfnew - emwfold + ionewfnew - ionewfold
    
else
    if (mol .eq. Nmol) then
        xr = xo
        yr = yo
        zr = zo
    else
        X(mol,:) = xold(:)
        Y(mol,:) = yold(:)
        Z(mol,:) = zold(:)
    endif
  
end if

return

end subroutine displace