remove.f90

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

		DAMPion_new( 1:stion-1 ) = DAMPion( 1:stion-1 )
		DAMPion_new( stion:Nion-lenion ) = DAMPion( endion+1:Nion )

	end if

end if

allocate( Xlj_grow(lenlj) )
allocate( Ylj_grow(lenlj) )
allocate( Zlj_grow(lenlj) )

allocate( TYPElj_grow(lenlj) )

allocate( DAMPlj2_grow(lenlj) )
allocate( DAMPlj3_grow(lenlj) )

Xlj_grow( 1:lenlj ) = Xlj( stlj:endlj )
Ylj_grow( 1:lenlj ) = Ylj( stlj:endlj )
Zlj_grow( 1:lenlj ) = Zlj( stlj:endlj )

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

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

ULJBOX_grow = 0.0

NGROUPS = 0
	
do k = 1, Nlj-lenlj

	NGROUPS( TYPElj_new(k) ) = NGROUPS( TYPElj_new(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

ULJ_MOL_grow = 0.0

UION_MOL_grow = 0.0

Uintra_grow = 0.0

if( lenion > 0 ) then

	allocate( Xion_grow(lenion) )
	allocate( Yion_grow(lenion) )
	allocate( Zion_grow(lenion) )

	allocate( TYPEion_grow(lenion) )

	allocate( DAMPion_grow(lenion) )

	Xion_grow( 1:lenion ) = Xion( stion:endion )
	Yion_grow( 1:lenion ) = Yion( stion:endion )
	Zion_grow( 1:lenion ) = Zion( stion:endion )

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

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

	UREAL_grow = 0.0

	Xion_grow = -Xion_grow
	Yion_grow = -Yion_grow
	Zion_grow = -Zion_grow

	call Surf_Move( lenion, Xion_grow, Yion_grow, Zion_grow, &
					TYPEion_grow, DAMPion_grow, &
					Nham, Niongrs, CHARGE, &
					BoxSize, SUMQX, SUMQY, SUMQZ, &
					SUMQX_new, SUMQY_new, SUMQZ_new, dUSURF )

	Xion_grow = -Xion_grow
	Yion_grow = -Yion_grow
	Zion_grow = -Zion_grow

	SUMQX_grow = SUMQX_new
	SUMQY_grow = SUMQY_new
	SUMQZ_grow = SUMQZ_new


	USURF_new = USURF + dUSURF * CoulCombo

	USURF_grow = USURF_new

	allocate( EXPX_grow(0:Kmax, lenion ) )
	allocate( EXPY_grow(-Kmax:Kmax, lenion ) )
	allocate( EXPZ_grow(-Kmax:Kmax, lenion ) )

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

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

	CHARGE = -CHARGE

	call Fourier_Move( lenion, Xion_grow, Yion_grow, Zion_grow, &
					   TYPEion_grow, DAMPion_grow, &
					   Nham, Niongrs, CHARGE, BoxSize, &
					   Kmax, Nkvec, KX, KY, KZ, &
					   CONST, CZERO1, CZERO2, CZERO2, &
					   EXPX_grow, EXPY_grow, EXPZ_grow, &
					   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_grow = 0.0
	USELF_MOL_grow = 0.0

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_grow = 0.0
	
end if

TMP = LNPSI

LnPi_old = LnPi

New = .false.

ZETA_tmp(:) = ZETA(SpID,:)

do i = 1, ERSTEPS(SpID)

	call Grow( EREND(i,SpID), STEPSTART(:,SpID), STEPLENGTH(:,SpID), &
			   NTRIALS(:,SpID), lenlj, lenion, MaxBeads, METHOD(:,SpID), &
			   MaxInt, MaxReal, INTPARAM(:,:,SpID), REALPARAM(:,:,SpID), &
			   BEADTYPE(:,SpID), New, ERSTART(i,SpID), &
			   Xlj_grow, Ylj_grow, Zlj_grow, &
			   TYPElj_grow, DAMPlj2_grow, DAMPlj3_grow, &
			   Xion_grow, Yion_grow, Zion_grow, &
			   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_grow, EXPY_grow, EXPZ_grow, SUMQEXPV_grow, &
			   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 )


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

	end do

	TMP = TMP - log( BIGW_old ) - &
			log( ZETA_tmp * ( BoxSize**3.0 ) / real( Nmol(SpID) ) ) / 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) )

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

	LnPi_old = LnPi_tmp

	if( i == ERSTEPS(SpID) ) then

		Success = .true.

		ULJBOX_new = ULJBOX - ULJBOX_grow

		UREAL_new = UREAL - UREAL_grow
		USELF_CH_new = USELF_CH - USELF_CH_grow

		dU(:) =	ULJBOX_new(:) + ULJLR_new(:) + 0.0 + &
				UREAL_new(:) + USURF_new(:) + UFOURIER_new(:) - &
				USELF_CH_new(:) - 0.0 + 0.0 - &
				( ULJBOX(:) + ULJLR(:) + ULJ_MOL(tag_mol,:) + &
				  UREAL(:) + USURF(:) + UFOURIER(:) - &
				  USELF_CH(:) - USELF_MOL(tag_mol,:) + UION_MOL(tag_mol,:) )

		dU = -dU * BETA

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

		Largest = maxval( LNW + LNPSI )

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

		Xlj = Xlj_new
		Ylj = Ylj_new
		Zlj = Zlj_new

		TYPElj = TYPElj_new

		DAMPlj2 = DAMPlj2_new
		DAMPlj3 = DAMPlj3_new

		ULJBOX = ULJBOX_new
		ULJLR = ULJLR_new

		if( lenion > 0 ) then

			Xion = Xion_new
			Yion = Yion_new
			Zion = Zion_new

			TYPEion = TYPEion_new

			DAMPion = DAMPion_new

			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

			if( stion == 1 ) then

				if( Nmol(0) == 1 ) then

					EXPX = ( 0.0, 0.0 )
					EXPY = ( 0.0, 0.0 )
					EXPZ = ( 0.0, 0.0 )

				else

					EXPX( :, 1:Nion-lenion ) = EXPX( :, endion+1:Nion )
					EXPY( :, 1:Nion-lenion ) = EXPY( :, endion+1:Nion )
					EXPZ( :, 1:Nion-lenion ) = EXPZ( :, endion+1:Nion )

				end if

			else if( stion + lenion - 1 == Nion ) then

				EXPX( :, stion:Nion ) = ( 0.0, 0.0 )
				EXPY( :, stion:Nion ) = ( 0.0, 0.0 )
				EXPZ( :, stion:Nion ) = ( 0.0, 0.0 )

			else

				EXPX( :, stion:Nion-lenion ) = EXPX( :, endion+1:Nion )
				EXPY( :, stion:Nion-lenion ) = EXPY( :, endion+1:Nion )
				EXPZ( :, stion:Nion-lenion ) = EXPZ( :, endion+1:Nion )

			end if

		end if


		if( tag_mol == 1 ) then

			if( Nmol(0) == 1 ) then

				SPECIES = 0

				LENGTHlj = 0
				LENGTHion = 0

				USELF_MOL = 0.0

				ULJ_MOL = 0.0 

				UION_MOL = 0.0 

				UINTRA = 0.0

				STARTlj = 0
				STARTion  = 0

			else

				SPECIES( 1:Nmol(0)-1 ) = SPECIES( 2:Nmol(0) )

				LENGTHlj( 1:Nmol(0)-1 ) = LENGTHlj( 2:Nmol(0) )
				LENGTHion( 1:Nmol(0)-1 ) = LENGTHion( 2:Nmol(0) )

				USELF_MOL( 1:Nmol(0)-1, : ) = USELF_MOL( 2:Nmol(0), : )

				ULJ_MOL( 1:Nmol(0)-1, : ) = ULJ_MOL( 2:Nmol(0), : )

				UION_MOL( 1:Nmol(0)-1, : ) = UION_MOL( 2:Nmol(0), : )

				UINTRA( 1:Nmol(0)-1 ) = UINTRA( 2:Nmol(0) )

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

				do k = 2, Nmol(0) - 1

					STARTlj(k) = STARTlj(k-1)	+ LENGTHlj(k-1)
					STARTion(k) = STARTion(k-1)	+ LENGTHion(k-1)

				end do

			end if

		else if( tag_mol == Nmol(0) ) then

			SPECIES( tag_mol ) = 0

			LENGTHlj( tag_mol ) = 0
			LENGTHion( tag_mol ) = 0

			USELF_MOL( tag_mol, : ) = 0.0

			ULJ_MOL( tag_mol, : ) = 0.0

			UION_MOL( tag_mol, : ) = 0.0

			UINTRA( tag_mol ) = 0.0

			STARTlj( tag_mol ) = 0
			STARTion( tag_mol ) = 0

		else

			SPECIES( tag_mol:Nmol(0)-1 ) = SPECIES( tag_mol+1:Nmol(0) )

			LENGTHlj( tag_mol:Nmol(0)-1 ) = LENGTHlj( tag_mol+1:Nmol(0) )
			LENGTHion( tag_mol:Nmol(0)-1 ) = LENGTHion( tag_mol+1:Nmol(0) )

			USELF_MOL( tag_mol:Nmol(0)-1, : ) = USELF_MOL( tag_mol+1:Nmol(0), : )

			ULJ_MOL( tag_mol:Nmol(0)-1, : ) = ULJ_MOL( tag_mol+1:Nmol(0), : )

			UION_MOL( tag_mol:Nmol(0)-1, : ) = UION_MOL( tag_mol+1:Nmol(0), : )

			UINTRA( tag_mol:Nmol(0)-1 ) = UINTRA( tag_mol+1:Nmol(0) )

			do k = tag_mol, Nmol(0) - 1

				STARTlj(k) = STARTlj(k-1)	+ LENGTHlj(k-1)
				STARTion(k) = STARTion(k-1)	+ LENGTHion(k-1)

			end do

		end if

		Nmol(0) = Nmol(0) - 1
		Nmol(SpID) = Nmol(SpID) - 1

		Nlj = Nlj - lenlj

		Nion = Nion - lenion

		tag_val = 0
		tag_mol = 0

	end if

end do

if( .NOT. success .AND. EESTEPS(SpID) == 1 ) then

	tag_val = 0
	tag_mol = 0

end if

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

deallocate( Xlj_grow )
deallocate( Ylj_grow )
deallocate( Zlj_grow )

deallocate( TYPElj_grow )

deallocate( DAMPlj2_grow )
deallocate( DAMPlj3_grow )

if( lenion > 0 ) then

	deallocate( Xion_grow )
	deallocate( Yion_grow )
	deallocate( Zion_grow )

	deallocate( TYPEion_grow )

	deallocate( DAMPion_grow )

	deallocate( EXPX_grow )
	deallocate( EXPY_grow )
	deallocate( EXPZ_grow )

end if

return

end subroutine Remove