methods.f90

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

subroutine BendTor( lenlj, lenion, &
					Xlj_new, Ylj_new, Zlj_new, &
					Xion_new, Yion_new, Zion_new, &
					LJBEAD, IONBEAD, &
					BEADTYPE, MaxInt, MaxReal, &
					INTPARAM, REALPARAM, &
					BoxSize, Nham, BETA, LNW, &
					FxSt, newold, &
					ncount, nDoOver, &
					xtr, ytr, ztr, &
					uvib, ubend, utor, Seed )  

implicit none

integer									:: lenlj, lenion

real, dimension(lenlj)					:: Xlj_new, Ylj_new, Zlj_new
real, dimension(lenion)					:: Xion_new, Yion_new, Zion_new

integer, dimension(lenlj+lenion)		:: LJBEAD, IONBEAD

character*5, dimension(lenlj+lenion)	:: BEADTYPE

integer									:: MaxInt, MaxReal

integer, dimension(MaxInt)				:: INTPARAM
real, dimension(MaxReal)				:: REALPARAM

real									:: BoxSize

integer									:: Nham

real, dimension(Nham)					:: BETA, LNW

integer									:: FxSt, newold, ncount, nDoOver

real									:: xtr, ytr, ztr
real									:: uvib, ubend, utor

integer									:: Seed

real, external							:: ran2

! Local Stuff

integer									:: m
integer									:: ntemp1, ntemp2
integer									:: noutfold, nbendp

logical									:: success

real									:: req, rtrial, utemp
real									:: Largest, Ratio
real, dimension(Nham)					:: dU
real, dimension(30)						:: Xc, Yc, Zc
real, dimension(1)						:: ZERO2 = 0.0



if( BEADTYPE( INTPARAM(4) ) == 'LJ' )	then

	Xc(3) = Xlj_new( LJBEAD( INTPARAM(4) ) )
	Yc(3) = Ylj_new( LJBEAD( INTPARAM(4) ) )
	Zc(3) = Zlj_new( LJBEAD( INTPARAM(4) ) )

else if( BEADTYPE( INTPARAM(4) ) == 'ION' ) then

	Xc(3) = Xion_new( IONBEAD( INTPARAM(4) ) )
	Yc(3) = Yion_new( IONBEAD( INTPARAM(4) ) )
	Zc(3) = Zion_new( IONBEAD( INTPARAM(4) ) )

end if

do m = 1, INTPARAM(1) + 2 * INTPARAM(2)

	if( BEADTYPE( INTPARAM(m+4) ) == 'LJ' )	then

		Xc(m+4) = Xlj_new( LJBEAD( INTPARAM(m+4) ) )
		Yc(m+4) = Ylj_new( LJBEAD( INTPARAM(m+4) ) )
		Zc(m+4) = Zlj_new( LJBEAD( INTPARAM(m+4) ) )

	else if( BEADTYPE( INTPARAM(m+4) ) == 'ION' ) then

		Xc(m+4) = Xion_new( IONBEAD( INTPARAM(m+4) ) )
		Yc(m+4) = Yion_new( IONBEAD( INTPARAM(m+4) ) )
		Zc(m+4) = Zion_new( IONBEAD( INTPARAM(m+4) ) )

	end if

end do

if( newold == 1 ) then

	Xc(4) = Xc(3)
	Yc(4) = Yc(3)
	Zc(4) = Zc(3)

else if( newold == 2 ) then

	Xc(4) = xtr
	Yc(4) = ytr
	Zc(4) = ztr

end if

noutfold = 2 + INTPARAM(1) + 2 * INTPARAM(2)

call Outfold( noutfold, 0, BoxSize, &
			  Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 )

if( FxSt == 1 ) then

	rtrial = REALPARAM(1)

	uvib = 0.0

else if( FxSt == 2 ) then

	req = REALPARAM(2)

	if( newold == 1 ) then

		call bondl( REALPARAM(1), req, BETA(1), Seed, rtrial )

		uvib = 0.5 * REALPARAM(1) * ( rtrial - req ) ** 2.0

	else if( newold == 2 ) then

		rtrial = sqrt( ( Xc(4) - Xc(3) ) ** 2 + ( Yc(4) - Yc(3) ) ** 2 + &
					   ( Zc(4) - Zc(3) ) ** 2 )

		uvib = 0.5 * REALPARAM(1) * ( rtrial - req ) ** 2.0

	end if

end if

success = .false.

ncount = 0

do while( .NOT. success )

	ncount = ncount + 1

	if( ncount > nDoOver ) return

	if( newold == 1 ) then

		call Sphere( Xc(3), Yc(3), Zc(3), rtrial, &
					 Xc(4), Yc(4), Zc(4), Seed )

	end if

	ubend = 0.0
	utor = 0.0
	
	do m = 1, INTPARAM(1)

		Xc(2) = Xc(m+4)
		Yc(2) = Yc(m+4)
		Zc(2) = Zc(m+4)

		if( FxSt == 1 ) then

			ntemp1 = 2 * m
			ntemp2 = 1 + 2 * m
			
		else if( FxSt == 2 ) then
		
			ntemp1 = 1 + 2 * m
			ntemp2 = 2 + 2 * m 
		 
		end if

		call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1), &
						REALPARAM(ntemp2), utemp )

		ubend = ubend + utemp

	end do

	nbendp = INTPARAM(1)

	do m = 1, INTPARAM(2)

		Xc(1) = Xc(2*m+nbendp+3)
		Yc(1) = Yc(2*m+nbendp+3)
		Zc(1) = Zc(2*m+nbendp+3)

		Xc(2) = Xc(2*m+nbendp+4)
		Yc(2) = Yc(2*m+nbendp+4)
		Zc(2) = Zc(2*m+nbendp+4)

		if( INTPARAM(3) == 1 ) then

			if( FxSt == 1 ) then

				ntemp1 = 2*nbendp-4+6*m
				ntemp2 = 2*nbendp+1+6*m
				
			else if( FxSt == 2 ) then
			
				ntemp1 = 2*nbendp-3+6*m
				ntemp2 = 2*nbendp+2+6*m
		
			end if

			call IntraTorsion( Xc, Yc, Zc, 1, &
							   REALPARAM(ntemp1:ntemp2), &
							   utemp )

		else if( INTPARAM(3) == 2 ) then

			if( FxSt == 1 ) then

				ntemp1 = 2*nbendp-2+4*m
				ntemp2 = 2*nbendp+1+4*m
				
			else if( FxSt == 2 ) then
			
				ntemp1 = 2*nbendp-1+4*m
				ntemp2 = 2*nbendp+2+4*m
			 
			end if

			call IntraTorsion( Xc, Yc, Zc, 2, &
							   REALPARAM(ntemp1:ntemp2), &
							   utemp )

		end if

		utor = utor + utemp

	end do

	if( newold == 1 ) then

		dU = - ( ubend + utor ) * BETA

		Largest = maxval( LNW + dU )

		Ratio = log( sum( exp( LNW + dU - Largest ) ) ) + Largest

		if( log( Ran2(Seed) ) < Ratio ) then
		
			success	= .true.

			ncount = 0

			xtr = Xc(4)
			ytr = Yc(4)
			ztr = Zc(4)

		end if

	else if( newold == 2 ) then
	
		success = .true.

	end if

end do

return

end subroutine Bendtor

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

subroutine Bend( lenlj, lenion, &
					Xlj_new, Ylj_new, Zlj_new, &
					Xion_new, Yion_new, Zion_new, &
					LJBEAD, IONBEAD, &
					BEADTYPE, MaxInt, MaxReal, &
					INTPARAM, REALPARAM, &
					BoxSize, Nham, BETA, LNW, &
					FxSt, newold, &
					ncount, nDoOver, &
					xtr, ytr, ztr, &
					uvib, ubend, Seed )  

implicit none

integer									:: lenlj, lenion

real, dimension(lenlj)					:: Xlj_new, Ylj_new, Zlj_new
real, dimension(lenion)					:: Xion_new, Yion_new, Zion_new

integer, dimension(lenlj+lenion)		:: LJBEAD, IONBEAD

character*5, dimension(lenlj+lenion)	:: BEADTYPE

integer									:: MaxInt, MaxReal

integer, dimension(MaxInt)				:: INTPARAM
real, dimension(MaxReal)				:: REALPARAM

real									:: BoxSize

integer									:: Nham

real, dimension(Nham)					:: BETA, LNW

integer									:: FxSt, newold, ncount, nDoOver

real									:: xtr, ytr, ztr
real									:: uvib, ubend

integer									:: Seed

real, external							:: ran2

! Local Stuff

integer									:: m
integer									:: ntemp1, ntemp2
integer									:: noutfold

logical									:: success

real									:: req, rtrial, utemp
real									:: Largest, Ratio
real, dimension(Nham)					:: dU
real, dimension(30)						:: Xc, Yc, Zc
real, dimension(1)						:: ZERO2 = 0.0


! INTPARAM(1) = # of bending potentials evaluated
! INTPARAM(2) = bead # from which the current bead is to be grown

if( BEADTYPE( INTPARAM(2) ) == 'LJ' )	then

	Xc(2) = Xlj_new( LJBEAD( INTPARAM(2) ) )
	Yc(2) = Ylj_new( LJBEAD( INTPARAM(2) ) )
	Zc(2) = Zlj_new( LJBEAD( INTPARAM(2) ) )

else if( BEADTYPE( INTPARAM(2) ) == 'ION' ) then

	Xc(2) = Xion_new( IONBEAD( INTPARAM(2) ) )
	Yc(2) = Yion_new( IONBEAD( INTPARAM(2) ) )
	Zc(2) = Zion_new( IONBEAD( INTPARAM(2) ) )

end if

do m = 1, INTPARAM(1)

	if( BEADTYPE( INTPARAM(m+2) ) == 'LJ' )	then

		Xc(m+3) = Xlj_new( LJBEAD( INTPARAM(m+2) ) )
		Yc(m+3) = Ylj_new( LJBEAD( INTPARAM(m+2) ) )
		Zc(m+3) = Zlj_new( LJBEAD( INTPARAM(m+2) ) )

	else if( BEADTYPE( INTPARAM(m+2) ) == 'ION' ) then

		Xc(m+3) = Xion_new( IONBEAD( INTPARAM(m+2) ) )
		Yc(m+3) = Yion_new( IONBEAD( INTPARAM(m+2) ) )
		Zc(m+3) = Zion_new( IONBEAD( INTPARAM(m+2) ) )

	end if

end do

if( newold == 1 ) then

	Xc(3) = Xc(2)
	Yc(3) = Yc(2)
	Zc(3) = Zc(2)

else if( newold == 2 ) then

	Xc(3) = xtr
	Yc(3) = ytr
	Zc(3) = ztr

end if

noutfold = INTPARAM(1) + 2

call Outfold( noutfold, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), &
			  ZERO2, ZERO2, ZERO2 )

if( FxSt == 1 ) then

	rtrial = REALPARAM(1)

	uvib = 0.0

else if( FxSt == 2 ) then

	req = REALPARAM(2)

	if( newold == 1 ) then

		call bondl( REALPARAM(1), req, BETA(1), Seed, rtrial )

		uvib = 0.5 * REALPARAM(1) * ( rtrial - req ) ** 2.0

	else if( newold == 2 ) then

		rtrial = sqrt( ( Xc(3) - Xc(2) ) ** 2 + ( Yc(3) - Yc(2) ) ** 2 + &
					   ( Zc(3) - Zc(2) ) ** 2 )

		uvib = 0.5 * REALPARAM(1) * ( rtrial - req ) ** 2.0

	end if

end if

success = .false.

ncount = 0

do while( .NOT. Success )

	ncount = ncount + 1

	if( ncount > nDoOver ) return

	if( newold == 1 ) then

		call Sphere( Xc(2), Yc(2), Zc(2), rtrial, &
					 Xc(3), Yc(3), Zc(3), Seed )

	end if

	ubend = 0.0
	
	do m = 1, INTPARAM(1)

		Xc(1) = Xc(m+3)
		Yc(1) = Yc(m+3)
		Zc(1) = Zc(m+3)

		if( FxSt == 1 ) then

			ntemp1 = 2 * m
			ntemp2 = 1 + 2 * m
			
		else if( FxSt == 2 ) then
		
			ntemp1 = 1 + 2 * m
			ntemp2 = 2 + 2 * m 
		 
		end if

		call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1), &
						REALPARAM(ntemp2), utemp )

		ubend = ubend + utemp

	end do

	if( newold == 1 ) then

		dU = - ubend * BETA

		Largest = maxval( LNW + dU )

		Ratio = log( sum( exp( LNW + dU - Largest ) ) ) + Largest

		if( log( ran2(Seed) ) < Ratio ) then
		
			success	= .true.

			ncount = 0

			xtr = Xc(3)
			ytr = Yc(3)
			ztr = Zc(3)

		end if

	else if( newold == 2 ) then

		success = .true.

	end if

end do

return

end subroutine Bend

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

subroutine Stretch( lenlj, lenion, &
					Xlj_new, Ylj_new, Zlj_new, &
					Xion_new, Yion_new, Zion_new, &
					LJBEAD, IONBEAD, &
					BEADTYPE, MaxInt, MaxReal, &
					INTPARAM, REALPARAM, &
					BoxSize, Nham, BETA, &
					FxSt, newold, &
					xtr, ytr, ztr, &
					uvib, Seed )  

implicit none

integer									:: lenlj, lenion

real, dimension(lenlj)					:: Xlj_new, Ylj_new, Zlj_new
real, dimension(lenion)					:: Xion_new, Yion_new, Zion_new

integer, dimension(lenlj+lenion)		:: LJBEAD, IONBEAD

character*5, dimension(lenlj+lenion)	:: BEADTYPE

integer									:: MaxInt, MaxReal

integer, dimension(MaxInt)				:: INTPARAM
real, dimension(MaxReal)				:: REALPARAM

real									:: BoxSize

integer									:: Nham

real, dimension(Nham)					:: BETA

integer									:: FxSt, newold

real									:: xtr, ytr, ztr
real									:: uvib

integer									:: Seed

! Local Stuff

real									:: req, rtrial
real, dimension(30)						:: Xc, Yc, Zc
real, dimension(1)						:: ZERO2 = 0.0


uvib = 0.0

if( BEADTYPE( INTPARAM(1) ) == 'LJ' )	then

	Xc(1) = Xlj_new( LJBEAD( INTPARAM(1) ) )
	Yc(1) = Ylj_new( LJBEAD( INTPARAM(1) ) )
	Zc(1) = Zlj_new( LJBEAD( INTPARAM(1) ) )

else if( BEADTYPE( INTPARAM(1) ) == 'ION' ) then

	Xc(1) = Xion_new( IONBEAD( INTPARAM(1) ) )
	Yc(1) = Yion_new( IONBEAD( INTPARAM(1) ) )
	Zc(1) = Zion_new( IONBEAD( INTPARAM(1) ) )

end if

if( newold == 1 ) then

	Xc(2) = Xc(1)
	Yc(2) = Yc(1)
	Zc(2) = Zc(1)

else if( newold == 2 ) then

	Xc(2) = xtr
	Yc(2) = ytr
	Zc(2) = ztr

end if

call Outfold( 2, 0, BoxSize, Xc, Yc, Zc, &
			  ZERO2, ZERO2, ZERO2 )

if( FxSt == 1 ) then

	rtrial = REALPARAM(1)

	uvib = 0.0

else if( FxSt == 2 ) then

	req = REALPARAM(2)

	if( newold == 1 ) then

		call bondl( REALPARAM(1), req, BETA(1), Seed, rtrial )

		uvib = 0.5 * REALPARAM(1) * ( rtrial - req ) ** 2.0

	else if( newold == 2 ) then

		rtrial = sqrt( ( Xc(2) - Xc(1) ) ** 2 + ( Yc(2) - Yc(1) ) ** 2 + &
					   ( Zc(2) - Zc(1) ) ** 2 )

		uvib = 0.5 * REALPARAM(1) * ( rtrial - req ) ** 2.0

	end if

end if

if( newold == 1 ) then

	call Sphere( Xc(1), Yc(1), Zc(1), rtrial, &
				 Xc(2), Yc(2), Zc(2), Seed )

end if

xtr = Xc(2)
ytr = Yc(2)
ztr = Zc(2)

return

end subroutine stretch

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

subroutine cone2( lenlj, lenion, nb, &
					Xlj_new, Ylj_new, Zlj_new, &
					Xion_new, Yion_new, Zion_new, &
					LJBEAD, IONBEAD, &
					BEADTYPE, MaxInt, MaxReal, &
					INTPARAM, REALPARAM, &
					BoxSize, &
					xtr, ytr, ztr, &
					Seed )  

implicit none

integer									:: lenlj, lenion, nb

real, dimension(lenlj)					:: Xlj_new, Ylj_new, Zlj_new
real, dimension(lenion)					:: Xion_new, Yion_new, Zion_new

integer, dimension(lenlj+lenion)		:: LJBEAD, IONBEAD

character*5, dimension(lenlj+lenion)	:: BEADTYPE

integer									:: MaxInt, MaxReal

integer, dimension(MaxInt)				:: INTPARAM
real, dimension(MaxReal)				:: REALPARAM

real									:: BoxSize

real, dimension(nb)						:: xtr, ytr, ztr

integer									:: Seed

! Local Stuff

integer									:: m
integer									:: ntemp1, ntemp2
integer									:: ntemp3, ntemp4

real, dimension(30)						:: Xc, Yc, Zc
real, dimension(1)						:: ZERO2 = 0.0



do m = 1, 2

	if( BEADTYPE( INTPARAM(m) ) == 'LJ' )	then

		Xc(m) = Xlj_new( LJBEAD( INTPARAM(m) ) )
		Yc(m) = Ylj_new( LJBEAD( INTPARAM(m) ) )
		Zc(m) = Zlj_new( LJBEAD( INTPARAM(m) ) )

	else if( BEADTYPE( INTPARAM(m) ) == 'ION' ) then

		Xc(m) = Xion_new( IONBEAD( INTPARAM(m) ) )
		Yc(m) = Yion_new( IONBEAD( INTPARAM(m) ) )
		Zc(m) = Zion_new( IONBEAD( INTPARAM(m) ) )

	end if

end do

call Outfold( 2, 0, BoxSize, Xc, Yc, Zc, &
			  ZERO2, ZERO2, ZERO2 )

ntemp1 = Nb+1
ntemp2 = 2*Nb
ntemp3 = 2*Nb+1
ntemp4 = 3*Nb

call Cone( Nb, Xc, Yc, Zc, REALPARAM(1:Nb), &
		   REALPARAM(ntemp1:ntemp2), REALPARAM(ntemp3:ntemp4), &
		   xtr, ytr, ztr, Seed )			 

return

end subroutine cone2