amoeba.f90

巨正则系综蒙特卡罗算法的源程序；可以用来进行吸附等分子模拟；最大的好处在于可以插入或删除原子

      SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funkk,iter, &
						Nmin, Nmx, Ubins, Umin, Uwidth, Nham, ham, MaxMol, &
						UN_HIST, BETA, ZETA, f1 )

      INTEGER iter,mp,ndim,np,NMAX,ITMAX
      REAL ftol,p(mp,np),y(mp),funkk
      PARAMETER (NMAX=20,ITMAX=5000)
      EXTERNAL funkk

	  integer									:: Nmin, Nmx, Ubins, MaxMol, Nham, ham

	  real										:: Umin, Uwidth

	  real, dimension(Ubins, 0:MaxMol, Nham)	:: UN_HIST
	  real, dimension(Nham)						:: BETA, ZETA

	  real, dimension(Nmin:Nmx)					:: f1

!     USES amotry,funkk
      INTEGER i,ihi,ilo,inhi,j,m,n
      REAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotry
      iter=0
1     do 12 n=1,ndim
        sum=0.
        do 11 m=1,ndim+1
          sum=sum+p(m,n)
11      continue
        psum(n)=sum
12    continue
2     ilo=1
      if (y(1).gt.y(2)) then
        ihi=1
        inhi=2
      else
        ihi=2

        inhi=1
      endif
      do 13 i=1,ndim+1
        if(y(i).le.y(ilo)) ilo=i
        if(y(i).gt.y(ihi)) then
          inhi=ihi
          ihi=i
        else if(y(i).gt.y(inhi)) then
          if(i.ne.ihi) inhi=i
        endif
13    continue
      rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
      if (rtol.lt.ftol) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue

        return
      endif
      if (iter.ge.ITMAX) pause 'ITMAX exceeded in amoeba'
      iter=iter+2
      ytry=amotry(p,y,psum,mp,np,ndim,funkk,ihi,-1.0, &
					  Nmin, Nmx, Ubins, Umin, Uwidth, Nham, ham, MaxMol, &
					  UN_HIST, BETA, ZETA, f1 )
      if (ytry.le.y(ilo)) then
        ytry=amotry(p,y,psum,mp,np,ndim,funkk,ihi,2.0, &
					  Nmin, Nmx, Ubins, Umin, Uwidth, Nham, ham, MaxMol, &
					  UN_HIST, BETA, ZETA, f1 )
      else if (ytry.ge.y(inhi)) then
        ysave=y(ihi)
        ytry=amotry(p,y,psum,mp,np,ndim,funkk,ihi,0.5, &
					  Nmin, Nmx, Ubins, Umin, Uwidth, Nham, ham, MaxMol, &
					  UN_HIST, BETA, ZETA, f1 )
        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
              y(i)=funkk(psum, Nmin, Nmx, Ubins, Umin, Uwidth, Nham, ham, MaxMol, &
						 UN_HIST, BETA, ZETA, f1 )
            endif
16        continue
          iter=iter+ndim
          goto 1

        endif
      else
        iter=iter-1
      endif
      goto 2
      END