weights.f90

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

subroutine Weights(	Nham, stop_h, LNW, PSI_PI )

implicit none

integer, intent(in)							:: Nham
integer, intent(in)							:: stop_h

real, dimension(Nham), intent(inout)		:: LNW
real, dimension(Nham), intent(in)			:: PSI_PI

! Local stuff

integer							:: i, hbase
real							:: psi_pi_crit, w_crit
real, dimension(Nham)			:: W

W = exp( LNW )

hbase = stop_h

if( hbase > Nham ) then

	hbase = 0

	psi_pi_crit = PSI_PI(1)
	w_crit = W(1)

else

	psi_pi_crit = PSI_PI(hbase)
	w_crit = W(hbase)

end if

do i = 1, Nham

	if( PSI_PI(i) > psi_pi_crit .OR. hbase == 0 ) then
		
		W(i) = psi_pi_crit / PSI_PI(i) * w_crit

	end if

end do

W = W / sum( W )

LNW = log( W )

return

end subroutine Weights


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


subroutine EEWeights( Nees, EEW, ALP_PROB )

implicit none

integer, intent(in)							:: Nees

real, dimension(Nees), intent(inout)		:: EEW
integer, dimension(0:Nees)					:: ALP_PROB

! Local stuff

integer							:: i, ref
real							:: eew_ref, sum
real, dimension(Nees)			:: W, dF

W = 0.0

ref = 0
i = 0

ALP_PROB(Nees) = ALP_PROB(0)

do while( ref == 0 )

	i = i + 1

	if( ALP_PROB(i) /= 0 ) ref = i

end do

eew_ref = EEW(ref)

dF = 0.0

do i = 1, Nees

	if( ALP_PROB(i) /= 0 ) then

		dF(i) = -log( real(ALP_PROB(i)) / real(ALP_PROB(ref)) ) + eew(i) - eew_ref

	end if

end do

do i = 1, Nees

	if( ALP_PROB(i) /= 0 ) then

		W(i) = exp( dF(i) + eew_ref )

	end if

end do

do i = 1, Nees

	if( ALP_PROB(i) == 0 ) then

		W(i) = maxval(W)

	end if

end do

EEW = log( W )

sum = 0.0

do i = 1, Nees

	sum = sum + EEW(i)

end do

sum = sum / real(Nees)

do i = 1, Nees

	EEW(i) = EEW(i) - sum

end do


return

end subroutine EEWeights