zbrent.f90

蒙特卡罗的一个程序分析 与大家分享 共同研究

! zbrent2 and funk2 are used to find the value of rmax/rm 

FUNCTION zbrent2(funk2,x1,x2,tol,alpha,cparam)
      INTEGER ITMAX
      REAL zbrent2,tol,x1,x2,funk2,EPS,alpha,cparam
      EXTERNAL funk2
      PARAMETER (ITMAX=100,EPS=3.e-8)
      INTEGER iter
      REAL a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
      a=x1
      b=x2
      fa=funk2(a,alpha,cparam)
      fb=funk2(b,alpha,cparam)

	  ! Start JRE addition

	  do while ( fb < 0.0 .AND. b > 0.1*x2 )

		b = 0.99*b
		fb=funk2(b,alpha,cparam)

	  end do

      if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.)) then  

		write(*,*) ' The Exp-6 potential has no repulsive part. '
		write(*,*) ' The program will stop. '
		stop

	  end if

	  ! End JRE addition

      c=b
      fc=fb
      do 11 iter=1,ITMAX
        if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then
          c=a
          fc=fa

          d=b-a
          e=d
        endif
        if(abs(fc).lt.abs(fb)) then
          a=b
          b=c
          c=a
          fa=fb
          fb=fc
          fc=fa
        endif
        tol1=2.*EPS*abs(b)+0.5*tol
        xm=.5*(c-b)
        if(abs(xm).le.tol1 .or. fb.eq.0.)then
          zbrent2=b
          return
        endif
        if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
          s=fb/fa
          if(a.eq.c) then
            p=2.*xm*s
            q=1.-s
          else

            q=fa/fc
            r=fb/fc
            p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
            q=(q-1.)*(r-1.)*(s-1.)
          endif
          if(p.gt.0.) q=-q
          p=abs(p)
          if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
            e=d
            d=p/q
          else
            d=xm
            e=d
          endif
        else
          d=xm
          e=d
        endif
        a=b
        fa=fb
        if(abs(d) .gt. tol1) then
          b=b+d
        else

          b=b+sign(tol1,xm)
        endif
        fb=funk2(b,alpha,cparam)
11    continue
      pause 'zbrent2 exceeding maximum iterations'
      zbrent2=b
      return
	  END function zbrent2


	function funk2( x, alpha, cparam )

	implicit none

	real				:: x, alpha, cparam, funk2

	funk2 = ( x ** 7.0 ) * exp( alpha * ( 1.0 - x ) ) - cparam

	end function funk2


! zbrent3 and funk3 are used to find the value of sigma/rm 

FUNCTION zbrent3(funk3,x1,x2,tol,alpha,cparam)
      INTEGER ITMAX
      REAL zbrent3,tol,x1,x2,funk3,EPS,alpha,cparam
      EXTERNAL funk3
      PARAMETER (ITMAX=100,EPS=3.e-8)
      INTEGER iter
      REAL a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
      a=x1
      b=x2
      fa=funk3(a,alpha,cparam)
      fb=funk3(b,alpha,cparam)

	  if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.))pause

      c=b
      fc=fb
      do 11 iter=1,ITMAX
        if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then
          c=a
          fc=fa

          d=b-a
          e=d
        endif
        if(abs(fc).lt.abs(fb)) then
          a=b
          b=c
          c=a
          fa=fb
          fb=fc
          fc=fa
        endif
        tol1=2.*EPS*abs(b)+0.5*tol
        xm=.5*(c-b)
        if(abs(xm).le.tol1 .or. fb.eq.0.)then
          zbrent3=b
          return
        endif
        if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
          s=fb/fa
          if(a.eq.c) then
            p=2.*xm*s
            q=1.-s
          else

            q=fa/fc
            r=fb/fc
            p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
            q=(q-1.)*(r-1.)*(s-1.)
          endif
          if(p.gt.0.) q=-q
          p=abs(p)
          if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
            e=d
            d=p/q
          else
            d=xm
            e=d
          endif
        else
          d=xm
          e=d
        endif
        a=b
        fa=fb
        if(abs(d) .gt. tol1) then
          b=b+d
        else

          b=b+sign(tol1,xm)
        endif
        fb=funk3(b,alpha,cparam)
11    continue
      pause 'zbrent3 exceeding maximum iterations'
      zbrent3=b
      return
	  END function zbrent3


	function funk3( x, alpha, cparam )

	implicit none

	real				:: x, alpha, cparam, funk3

	funk3 = ( x ** 6.0 ) * exp( alpha * ( 1.0 - x ) ) - cparam * alpha / 6.0

	end function funk3