regrow.f90

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

! the lines that are followed by ! n-alkane

if( ran2(seed) < 0.5 ) then						! n-alkane
												
	do i = 1, lenlj								! n-alkane

		Xlj_grow1(i) = Xlj( endlj - i + 1 )		! n-alkane
		Ylj_grow1(i) = Ylj( endlj - i + 1 )		! n-alkane
		Zlj_grow1(i) = Zlj( endlj - i + 1 )		! n-alkane

		Xlj_grow2(i) = Xlj( endlj - i + 1 )		! n-alkane
		Ylj_grow2(i) = Ylj( endlj - i + 1 )		! n-alkane
		Zlj_grow2(i) = Zlj( endlj - i + 1 )		! n-alkane

	end do 									    ! n-alkane

else 											! n-alkane

	Xlj_grow1( 1:lenlj ) = Xlj( stlj:endlj )
	Ylj_grow1( 1:lenlj ) = Ylj( stlj:endlj )
	Zlj_grow1( 1:lenlj ) = Zlj( stlj:endlj )

	Xlj_grow2( 1:lenlj ) = Xlj( stlj:endlj )
	Ylj_grow2( 1:lenlj ) = Ylj( stlj:endlj )
	Zlj_grow2( 1:lenlj ) = Zlj( stlj:endlj )
												
end if											! n-alkane

TYPElj_grow( 1:lenlj ) = TYPElj( stlj:endlj )

DAMPlj2_grow( 1:lenlj ) = DAMPlj2( stlj:endlj )
DAMPlj3_grow( 1:lenlj ) = DAMPlj3( stlj:endlj )

StartGrowthStep = int( ran2(Seed) * ( NSTEPS(SpID) - 1 ) ) + 2

! For n-alkanes only

StartGrowthStep = NSTEPS(SpID) / 2 + int( ran2(Seed) * ( NSTEPS(SpID) / 2 ) ) + 1  ! n-alkane

lenljst = 0
lenionst = 0

do i = StartGrowthStep, NSTEPS(SpID)

	do k = STEPSTART(i,SpID), STEPSTART(i,SpID) + STEPLENGTH(i,SpID) - 1
	
		if(	BEADTYPE(k,SpID) == 'LJ' ) then

			lenljst = lenljst + 1
		
		else if( BEADTYPE(k,SpID) == 'ION' ) then

			lenionst = lenionst + 1

		end if

	end do

end do

ULJLR_grow = ULJLR
ULJLR_new = ULJLR

if( lenljst > 0 ) then

	NGROUPS = 0
		
	do k = 1, Nlj-lenlj

		NGROUPS( TYPElj_new(k) ) = NGROUPS( TYPElj_new(k) ) + 1

	end do

	do k = 1, lenlj - lenljst

		NGROUPS( TYPElj_grow(k) ) = NGROUPS( TYPElj_grow(k) ) + 1

	end do

	do h = 1, Nham

		if( CP(1,1,h) > 0.0 ) then

			EPS_CP(:,:,h) = EPS(:,:,h) * CP(:,:,h) * ALP(:,:,h) / &
							( ALP(:,:,h) - 6.0 ) / 4.0

		else

			EPS_CP(:,:,h) = EPS(:,:,h) * ALP(:,:,h) / ( ALP(:,:,h) - 6.0 ) * &
							( ALP(:,:,h) / 6.0 ) ** ( 6.0 / ( ALP(:,:,h) - 6.0 ) ) / 4.0

		end if

	end do

	call lrcorr( Nljgrs, Nham, BoxSize, NGROUPS, EPS_CP, SIG, XLRCORR, &
				 ELRCORR, ULJLR_new )

	ULJLR_grow = ULJLR_new

end if

ULJBOX_new = 0.0
ULJBOX_grow = 0.0

ULJ_MOL_new = 0.0
ULJ_MOL_grow = 0.0

UION_MOL_new = 0.0
UION_MOL_grow = 0.0

Uintra_new = 0.0
Uintra_grow = 0.0

if( lenionst > 0 ) then

	allocate( Xion_grow1(lenion) )
	allocate( Yion_grow1(lenion) )
	allocate( Zion_grow1(lenion) )

	allocate( Xion_grow2(lenion) )
	allocate( Yion_grow2(lenion) )
	allocate( Zion_grow2(lenion) )

	allocate( TYPEion_grow(lenion) )

	allocate( DAMPion_grow(lenion) )

	Xion_grow1( 1:lenion ) = Xion( stion:endion )
	Yion_grow1( 1:lenion ) = Yion( stion:endion )
	Zion_grow1( 1:lenion ) = Zion( stion:endion )

	Xion_grow2( 1:lenion ) = Xion( stion:endion )
	Yion_grow2( 1:lenion ) = Yion( stion:endion )
	Zion_grow2( 1:lenion ) = Zion( stion:endion )

	TYPEion_grow( 1:lenion ) = TYPEion( stion:endion )

	DAMPion_grow( 1:lenion ) = DAMPion( stion:endion )

	UREAL_new = 0.0
	UREAL_grow = 0.0

	Xion_grow1 = -Xion_grow1
	Yion_grow1 = -Yion_grow1
	Zion_grow1 = -Zion_grow1

	call Surf_Move( lenionst, Xion_grow1(lenion-lenionst+1:lenion), &
					Yion_grow1(lenion-lenionst+1:lenion), &
					Zion_grow1(lenion-lenionst+1:lenion), &
					TYPEion_grow(lenion-lenionst+1:lenion), &
					DAMPion_grow(lenion-lenionst+1:lenion), &
					Nham, Niongrs, CHARGE, &
					BoxSize, SUMQX, SUMQY, SUMQZ, &
					SUMQX_new, SUMQY_new, SUMQZ_new, dUSURF )

	Xion_grow1 = -Xion_grow1
	Yion_grow1 = -Yion_grow1
	Zion_grow1 = -Zion_grow1

	SUMQX_grow = SUMQX_new
	SUMQY_grow = SUMQY_new
	SUMQZ_grow = SUMQZ_new


	USURF_new = USURF + dUSURF * CoulCombo

	USURF_grow = USURF_new

	allocate( EXPX_grow1(0:Kmax, lenion ) )
	allocate( EXPY_grow1(-Kmax:Kmax, lenion ) )
	allocate( EXPZ_grow1(-Kmax:Kmax, lenion ) )

	allocate( EXPX_grow2(0:Kmax, lenion ) )
	allocate( EXPY_grow2(-Kmax:Kmax, lenion ) )
	allocate( EXPZ_grow2(-Kmax:Kmax, lenion ) )

	allocate( CZERO1(0:Kmax,lenion) )
	allocate( CZERO2(-Kmax:Kmax,lenion) )

	EXPX_grow1(:, 1:lenion ) = EXPX( :, stion:endion )
	EXPY_grow1(:, 1:lenion ) = EXPY( :, stion:endion )
	EXPZ_grow1(:, 1:lenion ) = EXPZ( :, stion:endion )

	EXPX_grow2 = EXPX_grow1
	EXPY_grow2 = EXPY_grow1
	EXPZ_grow2 = EXPZ_grow1

	CZERO1 = ( 0.0, 0.0 )
	CZERO2 = ( 0.0, 0.0 )

	CHARGE = -CHARGE

	call Fourier_Move( lenionst, Xion_grow1(lenion-lenionst+1:lenion), &
					   Yion_grow1(lenion-lenionst+1:lenion), &
					   Zion_grow1(lenion-lenionst+1:lenion), &
					   TYPEion_grow(lenion-lenionst+1:lenion), &
					   DAMPion_grow(lenion-lenionst+1:lenion), &
					   Nham, Niongrs, CHARGE, BoxSize, &
					   Kmax, Nkvec, KX, KY, KZ, &
					   CONST, CZERO1, CZERO2, CZERO2, &
					   EXPX_grow1(:,lenion-lenionst+1:lenion), &
					   EXPY_grow1(:,lenion-lenionst+1:lenion), &
					   EXPZ_grow1(:,lenion-lenionst+1:lenion), &
					   SUMQEXPV, SUMQEXPV_new, dUFOURIER )

	CHARGE = -CHARGE

	SUMQEXPV_grow = SUMQEXPV_new

	UFOURIER_new = UFOURIER + dUFOURIER * CoulCombo

	UFOURIER_grow = UFOURIER_new

	deallocate( CZERO1 )
	deallocate( CZERO2 )

	USELF_CH_new = 0.0
	USELF_CH_grow = 0.0

	call SelfMolecule( lenion-lenionst, Xion_grow1, &
					   Yion_grow1, Zion_grow1, &
					   TYPEion_grow, DAMPion_grow, &
					   Nham, Niongrs, CHARGE, BoxSize, &
					   Alpha, USELF_MOL_new )

	USELF_MOL_new = USELF_MOL_new * CoulCombo
	
	USELF_MOL_grow = USELF_MOL_new
	
else

	UREAL_new = 0.0
	UREAL_grow = 0.0
	USURF_new = 0.0
	USURF_grow = 0.0
	UFOURIER_new = 0.0
	UFOURIER_grow = 0.0
	USELF_CH_new = 0.0
	USELF_CH_grow = 0.0
	USELF_MOL_new = 0.0
	USELF_MOL_grow = 0.0
	
end if

TMP = LNPSI

LnPi_old = LnPi

! The Regrowths can be done in two ways:
! 1. With early rejection (ER) used after each 
!    configurational-bias step.
!    For this option, comment ER OFF lines,
!    uncomment ER ON lines.
! 2. With no early rejection.
!    For this option, comment ER ON lines,
!    uncomment ER OFF lines.

do i = StartGrowthStep, NSTEPS(SpID)	 ! ER ON

	New = .False.

	call Grow( i, STEPSTART(:,SpID), STEPLENGTH(:,SpID), &			   ! ER ON
!	call Grow( NSTEPS(SpID), STEPSTART(:,SpID), STEPLENGTH(:,SpID), &  ! ER OFF
			   NTRIALS(:,SpID), lenlj, lenion, MaxBeads, &
			   METHOD(:,SpID), MaxInt, MaxReal, &
			   INTPARAM(:,:,SpID), REALPARAM(:,:,SpID), &
			   BEADTYPE(:,SpID), New, i, &							   ! ER ON
!			   BEADTYPE(:,SpID), New, StartGrowthStep, &			   ! ER OFF
			   Xlj_grow1, Ylj_grow1, Zlj_grow1, &
			   TYPElj_grow, DAMPlj2_grow, DAMPlj3_grow, &
			   Xion_grow1, Yion_grow1, Zion_grow1, &
			   TYPEion_grow, DAMPion_grow, &
			   Nham, ULJBOX_new, ULJLR_new, ULJ_MOL_new, &
			   UREAL_new, USURF_new, UFOURIER_new, USELF_CH_new, &
			   USELF_MOL_new, UION_MOL_new, Uintra_new, &
			   Niongrs, CHARGE, Alpha, &
			   Kmax, Nkvec, KX, KY, KZ, CONST, &
			   SUMQX_new, SUMQY_new, SUMQZ_new, &
			   EXPX_grow1, EXPY_grow1, EXPZ_grow1, &
			   SUMQEXPV_new, BETA, BIGW_old, LNW, &
			   Nlj-lenlj, Xlj_new, Ylj_new, Zlj_new, &
			   TYPElj_new, DAMPlj2_new, DAMPlj3_new, &
			   Nion-lenion, Xion_new, Yion_new, Zion_new, &
			   TYPEion_new, DAMPion_new, &
			   BoxSize, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
			   XLRCORR, ELRCORR, &
			   Nres, Nresmol, reslen, Xresmols, &
			   Yresmols, Zresmols, Uint_resm, &
			   Ulj_resm, Uion_resm, &
			   BONDSAPART(:,:,SpID), f14_lj, f14_ion, Seed )


	New = .True.

	call Grow( i, STEPSTART(:,SpID), STEPLENGTH(:,SpID), &				 ! ER ON
!	call Grow( NSTEPS(SpID), STEPSTART(:,SpID), STEPLENGTH(:,SpID), &	 ! ER OFF
			   NTRIALS(:,SpID), lenlj, lenion, MaxBeads, &
			   METHOD(:,SpID), MaxInt, MaxReal, &
			   INTPARAM(:,:,SpID), REALPARAM(:,:,SpID), &
			   BEADTYPE(:,SpID), New, i, &								 ! ER ON
!			   BEADTYPE(:,SpID), New, StartGrowthStep, &				 ! ER OFF
			   Xlj_grow2, Ylj_grow2, Zlj_grow2, &
			   TYPElj_grow, DAMPlj2_grow, DAMPlj3_grow, &
			   Xion_grow2, Yion_grow2, Zion_grow2, &
			   TYPEion_grow, DAMPion_grow, &
			   Nham, ULJBOX_grow, ULJLR_grow, ULJ_MOL_grow, &
			   UREAL_grow, USURF_grow, UFOURIER_grow, &
			   USELF_CH_grow, USELF_MOL_grow, &
			   UION_MOL_grow, Uintra_grow, &
			   Niongrs, CHARGE, Alpha, &
			   Kmax, Nkvec, KX, KY, KZ, CONST, &
			   SUMQX_grow, SUMQY_grow, SUMQZ_grow, &
			   EXPX_grow2, EXPY_grow2, EXPZ_grow2, &
			   SUMQEXPV_grow, BETA, BIGW_new, LNW, &
			   Nlj-lenlj, Xlj_new, Ylj_new, Zlj_new, &
			   TYPElj_new, DAMPlj2_new, DAMPlj3_new, &
			   Nion-lenion, Xion_new, Yion_new, Zion_new, &
			   TYPEion_new, DAMPion_new, &
			   BoxSize, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
			   XLRCORR, ELRCORR, &
			   Nres, Nresmol, reslen, Xresmols, &
			   Yresmols, Zresmols, Uint_resm, &
			   Ulj_resm, Uion_resm, &
			   BONDSAPART(:,:,SpID), f14_lj, f14_ion, Seed )

	if( minval( BIGW_new ) < 1.0e-250 ) then

		LnPi_tmp = -1.0e10

	else

		TMP = TMP + log( BIGW_new / BIGW_old ) 

		Largest = maxval( LNW + TMP )

		LnPi_tmp = log( sum( exp( LNW + TMP - Largest ) ) ) + Largest

	end if

	if( ran2(Seed) < 1.0 - exp(LnPi_tmp - LnPi_old) ) exit

	LnPi_old = LnPi_tmp

	if( i == NSTEPS(SpID) ) then		 ! ER ON

		Success = .True.

		ULJBOX_grow = ULJBOX + ULJBOX_grow - ULJBOX_new
		ULJ_MOL_grow(:) = ULJ_MOL(Mol,:) + ULJ_MOL_grow(:) - ULJ_MOL_new(:)

		UREAL_grow = UREAL + UREAL_grow - UREAL_new
		USELF_CH_grow = USELF_CH
		UION_MOL_grow(:) = UION_MOL(Mol,:) + UION_MOL_grow(:) - UION_MOL_new(:)

		Uintra_grow	= UINTRA(Mol) + Uintra_grow - Uintra_new

		dU(:) = ULJBOX_grow(:) + ULJLR_grow(:) + ULJ_MOL_grow(:) + &
				UREAL_grow(:) + USURF_grow(:) + UFOURIER_grow(:) - &
				USELF_CH_grow(:) - USELF_MOL_grow(:) + UION_MOL_grow(:) - &
				( ULJBOX(:) + ULJLR(:) + ULJ_MOL(Mol,:) + &
				  UREAL(:) + USURF(:) + UFOURIER(:) - &
				  USELF_CH(:) - USELF_MOL(Mol,:) + UION_MOL(Mol,:) )

 		dU = -dU * BETA

		LNPSI = LNPSI + dU 

		Largest = maxval( LNW + LNPSI )

		LnPi = log( sum( exp( LNW + LNPSI - Largest ) ) ) + Largest

		Xlj( stlj:endlj ) = Xlj_grow2( 1:lenlj )
		Ylj( stlj:endlj ) = Ylj_grow2( 1:lenlj )
		Zlj( stlj:endlj ) = Zlj_grow2( 1:lenlj )

		ULJBOX = ULJBOX_grow
!		ULJLR = ULJLR_grow             ! ULJLR does not change during a regrowth.

		if( lenion > 0 ) then

			Xion( stion:endion ) = Xion_grow2( 1:lenion )
			Yion( stion:endion ) = Yion_grow2( 1:lenion )
			Zion( stion:endion ) = Zion_grow2( 1:lenion )

			UREAL = UREAL_grow
			USURF = USURF_grow
			UFOURIER = UFOURIER_grow
!			USELF_CH = USELF_CH_grow	! USELF_CH does not change during a regrowth.

			SUMQX = SUMQX_grow
			SUMQY = SUMQY_grow
			SUMQZ = SUMQZ_grow

			SUMQEXPV = SUMQEXPV_grow

			EXPX( :, stion:endion ) = EXPX_grow2( :, 1:lenion )
			EXPY( :, stion:endion ) = EXPY_grow2( :, 1:lenion )
			EXPZ( :, stion:endion ) = EXPZ_grow2( :, 1:lenion )

		end if

		USELF_MOL(Mol,:) = USELF_MOL_grow(:)

		ULJ_MOL(Mol,:) = ULJ_MOL_grow(:)

		UION_MOL(Mol,:) = UION_MOL_grow(:)

		UINTRA(Mol) = Uintra_grow

	end if						 ! ER ON

end do							 ! ER ON

deallocate( Xlj_grow1 )
deallocate( Ylj_grow1 )
deallocate( Zlj_grow1 )

deallocate( Xlj_grow2 )
deallocate( Ylj_grow2 )
deallocate( Zlj_grow2 )

deallocate( TYPElj_grow )

deallocate( DAMPlj2_grow )
deallocate( DAMPlj3_grow )

if( lenion > 0 ) then

	deallocate( Xion_grow1 )
	deallocate( Yion_grow1 )
	deallocate( Zion_grow1 )

	deallocate( Xion_grow2 )
	deallocate( Yion_grow2 )
	deallocate( Zion_grow2 )

	deallocate( TYPEion_grow )

	deallocate( DAMPion_grow )

	deallocate( EXPX_grow1 )
	deallocate( EXPY_grow1 )
	deallocate( EXPZ_grow1 )

	deallocate( EXPX_grow2 )
	deallocate( EXPY_grow2 )
	deallocate( EXPZ_grow2 )

end if

return

end subroutine ReGrow