bondl.f90

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

! Taken from Frenkel and Smit.
! bondl ==> Algorithm 41, page 410
! gauss ==> Algorithm 42, page 411

subroutine bondl( Kr, req, beta, Seed, rtrial )

implicit none

real, intent(in)		:: Kr
real, intent(in)		:: req
real, intent(in)		:: beta

integer					:: Seed

real, intent(out)		:: rtrial

real, external			:: ran2

! Local stuff

real					:: sigma, a
logical					:: ready


!sigma = sqrt( 1.0 / ( 2.0 * beta * Kr ) )
sigma = sqrt( 1.0 / ( beta * Kr ) )

a = ( req + 3.0 * sigma ) ** 2.0

ready = .false.

do while ( .NOT. ready )

	call gauss( sigma, req, Seed, rtrial )

	if( ran2(Seed) < rtrial*rtrial/a ) ready = .true.

end do

return

end

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

subroutine angle( Kth, theq, beta, Seed, thtr )

implicit none

real, intent(in)		:: Kth
real, intent(in)		:: theq
real, intent(in)		:: beta

integer					:: Seed

real, intent(out)		:: thtr

real, external			:: ran2

! Local stuff

real					:: sigma

!sigma = sqrt( 1.0 / ( 2.0 * beta * Kr ) )
sigma = sqrt( 1.0 / ( beta * Kth ) )

! if in 3-d, sin(theta) must be taken into account.

call gauss( sigma, theq, Seed, thtr )

return

end

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

subroutine gauss( sigma, req, Seed, rtrial )

implicit none

real, intent(in)		:: sigma
real, intent(in)		:: req

integer					:: Seed

real, intent(out)		:: rtrial

real, external			:: ran2

! Local stuff

real					:: r, v1, v2

r = 2.0

do while( r > 1.0 )

	v1 = 2.0 * ran2(Seed) - 1.0
	v2 = 2.0 * ran2(Seed) - 1.0
	r = v1 * v1 + v2 * v2

end do

rtrial = v1 * sqrt( -2.0 * log(r) / r )
rtrial = req + sigma * rtrial

return

end