add.f90

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

Success = .False.

tag_val = 0
tag_mol = 0

llj = LENLJ(SpID)

lion = LENION(SpID)

allocate( Xlj_new(llj) )
allocate( Ylj_new(llj) )
allocate( Zlj_new(llj) )

allocate( TYPElj_new(llj) )

allocate( DAMPlj2_new(llj) )
allocate( DAMPlj3_new(llj) )

Xlj_new = 0.0
Ylj_new = 0.0
Zlj_new = 0.0

TYPElj_new( 1:llj ) = GROUPTYPE_LJ( 1:llj, SpID)

ULJBOX_new = ULJBOX
ULJLR_new = ULJLR
ULJ_MOL_new = 0.0

Uintra_new = 0.0

if( lion > 0 ) then

	allocate( Xion_new(lion) )
	allocate( Yion_new(lion) )
	allocate( Zion_new(lion) )

	allocate( TYPEion_new(lion) )

	allocate( EXPX_new(0:Kmax, lion ) )
	allocate( EXPY_new(-Kmax:Kmax, lion ) )
	allocate( EXPZ_new(-Kmax:Kmax, lion ) )

	allocate( DAMPion_new(lion) )

	Xion_new = 0.0
	Yion_new = 0.0
	Zion_new = 0.0

	TYPEion_new( 1:lion ) = GROUPTYPE_ION( 1:lion, SpID)

end if

UREAL_new = UREAL
USURF_new = USURF
UFOURIER_new = UFOURIER
USELF_CH_new = USELF_CH
USELF_MOL_new = 0.0
UION_MOL_new = 0.0

SUMQX_new = SUMQX
SUMQY_new = SUMQY
SUMQZ_new = SUMQZ

SUMQEXPV_new = SUMQEXPV

TMP = LNPSI

LnPi_old = LnPi

ljbk = 0
ionbk = 0

do i = 1, lion + llj
	
	if(	BEADTYPE(i,SpID) == 'LJ' ) then

		ljbk = ljbk + 1

		DAMPlj2_new( ljbk ) = sqrt( BEADDAMP( i, 1, SpID ) )
		DAMPlj3_new( ljbk ) = BEADDAMP( i, 1, SpID ) ** (1.0/3.0)

	else if( BEADTYPE(i,SpID) == 'ION' ) then

		ionbk = ionbk + 1

		DAMPion_new( ionbk ) = sqrt( BEADDAMP( i, 1, SpID ) )

	end if

end do

New = .True.

ZETA_tmp(:) = ZETA(SpID,:)

do i = 1, ERSTEPS(SpID)
	
	call Grow( EREND(i,SpID), STEPSTART(:,SpID), STEPLENGTH(:,SpID), &
			   NTRIALS(:,SpID), llj, lion, MaxBeads, METHOD(:,SpID), &
			   MaxInt, MaxReal, INTPARAM(:,:,SpID), REALPARAM(:,:,SpID), &
			   BEADTYPE(:,SpID), New, ERSTART(i,SpID), &
			   Xlj_new, Ylj_new, Zlj_new, &
			   TYPElj_new, DAMPlj2_new, DAMPlj3_new, &
			   Xion_new, Yion_new, Zion_new, &
			   TYPEion_new, DAMPion_new, &
			   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_new, EXPY_new, EXPZ_new, SUMQEXPV_new, &
			   BETA, BIGW_new, LNW, &
			   Nlj, Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
			   Nion, Xion, Yion, Zion, TYPEion, DAMPion, &
			   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 )

	do j = ERSTART(i,SpID), EREND(i,SpID)
	
		BIGW_new = BIGW_new / real( NTRIALS(j,SpID) )

	end do

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

		LnPi_tmp = -1.0e10

	else

		TMP = TMP + log( BIGW_new ) + &
				log( ZETA_tmp * ( BoxSize**3.0 ) / real( Nmol(SpID) + 1 ) ) / real( ERSTEPS(SpID) )
		
		Largest = maxval( LNW + TMP ) 

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

		LnPi_tmp = LnPi_tmp + ( EEW(1,SpID) - EEW(EESTEPS(SpID),SpID) ) / real( ERSTEPS(SpID) )

	end if

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

	LnPi_old = LnPi_tmp

	if( i == ERSTEPS(SpID) ) then

		Success = .true.

		dU = ULJBOX_new + ULJLR_new + ULJ_MOL_new + &
			 UREAL_new + USURF_new + UFOURIER_new - &
			 USELF_CH_new - USELF_MOL_new + UION_MOL_new - &
			 ( ULJBOX + ULJLR + 0.0 + &
			   UREAL + USURF + UFOURIER - &
			   USELF_CH - 0.0 + 0.0 )

		dU = -dU * BETA

		LNPSI = LNPSI + dU + &
				log( ZETA_tmp * ( BoxSize**3.0 ) / real( Nmol(SpID) + 1 ) )

		Largest = maxval( LNW + LNPSI )

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


		Xlj( Nlj+1:Nlj+llj ) = Xlj_new( 1:llj )
		Ylj( Nlj+1:Nlj+llj ) = Ylj_new( 1:llj )
		Zlj( Nlj+1:Nlj+llj ) = Zlj_new( 1:llj )

		TYPElj( Nlj+1:Nlj+llj ) = TYPElj_new( 1:llj )

		DAMPlj2( Nlj+1:Nlj+llj ) = DAMPlj2_new( 1:llj )
		DAMPlj3( Nlj+1:Nlj+llj ) = DAMPlj3_new( 1:llj )

		ULJBOX = ULJBOX_new
		ULJLR = ULJLR_new

		if( lion > 0 ) then

			Xion( Nion+1:Nion+lion ) = Xion_new( 1:lion )
			Yion( Nion+1:Nion+lion ) = Yion_new( 1:lion )
			Zion( Nion+1:Nion+lion ) = Zion_new( 1:lion )

			TYPEion( Nion+1:Nion+lion ) = TYPEion_new( 1:lion )

			DAMPion( Nion+1:Nion+lion ) = DAMPion_new( 1:lion )

			UREAL = UREAL_new
			USURF = USURF_new
			UFOURIER = UFOURIER_new
			USELF_CH = USELF_CH_new

			SUMQX = SUMQX_new
			SUMQY = SUMQY_new
			SUMQZ = SUMQZ_new

			SUMQEXPV = SUMQEXPV_new

			EXPX( :, Nion+1:Nion+lion ) = EXPX_new( :, 1:lion )
			EXPY( :, Nion+1:Nion+lion ) = EXPY_new( :, 1:lion )
			EXPZ( :, Nion+1:Nion+lion ) = EXPZ_new( :, 1:lion )

		end if

		tag_val = tag_val + 1
		tag_mol = Nmol(0) + 1

		if( EESTEPS(SpID) == tag_val ) then
			
			tag_val = 0
			tag_mol = 0

		end if

		SPECIES( Nmol(0)+1 ) = SpID

		LENGTHlj( Nmol(0)+1 ) = llj
		LENGTHion( Nmol(0)+1 ) = lion

		if( Nmol(0) == 0 ) then

			STARTlj(1) = 1
			STARTion(1) = 1

		else

			STARTlj( Nmol(0)+1 ) = STARTlj( Nmol(0) ) + LENGTHlj( Nmol(0) )
			STARTion( Nmol(0)+1 ) = STARTion( Nmol(0) ) + LENGTHion( Nmol(0) )

		end if

		USELF_MOL( Nmol(0)+1,: ) = USELF_MOL_new(:)

		ULJ_MOL( Nmol(0)+1,: ) = ULJ_MOL_new(:)

		UION_MOL( Nmol(0)+1,: ) = UION_MOL_new(:)

		UINTRA( Nmol(0)+1 ) = Uintra_new

		Nmol(0) = Nmol(0) + 1
		Nmol(SpID) = Nmol(SpID) + 1

		Nlj = Nlj + llj

		Nion = Nion + lion

	end if

end do

ER_HIST(i, SpID, 1) = ER_HIST(i, SpID, 1) + 1

deallocate( Xlj_new )
deallocate( Ylj_new )
deallocate( Zlj_new )

deallocate( TYPElj_new )

deallocate( DAMPlj2_new )
deallocate( DAMPlj3_new )

if( lion > 0 ) then

	deallocate( Xion_new )
	deallocate( Yion_new )
	deallocate( Zion_new )

	deallocate( TYPEion_new )

	deallocate( DAMPion_new )

	deallocate( EXPX_new )
	deallocate( EXPY_new )
	deallocate( EXPZ_new )

end if

return

end subroutine Add