predict.f90

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

subroutine predict( Nmin, Nmax, Ubins, Umin, Uwidth, Nham, &
					MaxMol, UN_HIST, BETA, ZETA )

implicit none

integer, intent(in)						:: Nmin, Nmax, Ubins, Nham, MaxMol

real, intent(in)						:: Umin, Uwidth

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

real, dimension(Nham)					:: BETA, ZETA

! Local Stuff

integer, parameter						:: ndim = 2
integer									:: h
integer									:: iter

real									:: lamda_b = 0.0001
real									:: lamda_z = 0.0001
real									:: ftol = 1.0e-6

real, dimension(ndim+1,ndim)			:: p
real, dimension(ndim+1)					:: y

real, dimension(2)						:: x1
real, dimension(Nmin:Nmax)				:: f1

real, external							:: funkk


x1(1) = BETA(1)
x1(2) = ZETA(1)

call fdist( f1, x1, Nmin, Nmax, Ubins, Umin, Uwidth, Nham, 1, MaxMol, &
			UN_HIST, BETA, ZETA )

write(*,*) 
write(*,*) ' The new temperatures and activities are: '
write(*,*) ' Ham.    Temp      Activity '

do h = 2, Nham

	p(1,1) = BETA(h)
	p(1,2) = ZETA(h)

	p(2,1) = BETA(h) + lamda_b
	p(2,2) = ZETA(h)

	p(3,1) = BETA(h)
	p(3,2) = ZETA(h) + lamda_z

	y(1) = funkk( p(1,:), Nmin, Nmax, Ubins, Umin, Uwidth, Nham, h, MaxMol, &
				  UN_HIST, BETA, ZETA, f1 )
	y(2) = funkk( p(2,:), Nmin, Nmax, Ubins, Umin, Uwidth, Nham, h, MaxMol, &
				  UN_HIST, BETA, ZETA, f1 )
	y(3) = funkk( p(3,:), Nmin, Nmax, Ubins, Umin, Uwidth, Nham, h, MaxMol, &
				  UN_HIST, BETA, ZETA, f1 )

	call amoeba( p, y, ndim+1, ndim, ndim, ftol, funkk, iter, &
				 Nmin, Nmax, Ubins, Umin, Uwidth, Nham, h, MaxMol, &
				 UN_HIST, BETA, ZETA, f1 )

	BETA(h) = p(1,1)
	ZETA(h) = p(1,2)

	write(*,"(1x,I3,4x,F8.4,4x,F8.4)") h, 1.0/BETA(h), log(ZETA(h))

end do

return

end subroutine predict

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine fdist( f, x, Nmin, Nmax, Ubins, Umin, Uwidth, Nham, ham, MaxMol, &
				  UN_HIST, BETA, ZETA )

implicit none

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

real									:: Umin, Uwidth

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

real, dimension(Nmin:Nmax)				:: f
real, dimension(2)						:: x

! Local

integer									:: i, j

real									:: energy, Xi, lnS

lnS = 0.0
Xi = 0.0
f = 0.0

do i = Nmin, Nmax

	do j = 1, Ubins

		if( UN_HIST(j,i,ham) > 0.0 ) then

			energy = Umin + ( real(j) - 0.5 ) * Uwidth

			lnS = log( UN_HIST(j,i,ham) ) + ( BETA(ham) - x(1) ) * energy - &
					real(i) * ( log(ZETA(ham)) - log(x(2)) )

			lnS = exp( lnS )

			f(i) = f(i) + lnS 

			Xi = Xi + lnS

		end if

	end do

end do

f = f / Xi

return

end subroutine fdist

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function funkk( x, Nmin, Nmax, Ubins, Umin, Uwidth, Nham, ham, MaxMol, &
				UN_HIST, BETA, ZETA, f1 )

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

real									:: Umin, Uwidth

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

real, dimension(Nmin:Nmax)				:: f1

real									:: funkk
real, dimension(2)						:: x

! Local

real, dimension(Nmin:Nmax)				:: f

integer									:: i		

if( x(2) < 0.0 ) then

	funkk = 1.0e100

else

	call fdist( f, x, Nmin, Nmax, Ubins, Umin, Uwidth, Nham, ham, MaxMol, &
				UN_HIST, BETA, ZETA )

	funkk = 0.0

	do i = Nmin, Nmax

		funkk = funkk + abs( f(i) - f1(i) )

	end do

end if

return

end function funkk