displace.f90

蒙特卡罗的一个程序分析 与大家分享 共同研究

		
		case ( 1 )			! Rotation in y-z plane.

			do i = 1, lenlj
				
				T(1) = Ylj_old(i) - ycom
				T(2) = Zlj_old(i) - zcom

				T = matmul( M, T )

				Ylj_new(i) = T(1) + ycom
				Zlj_new(i) = T(2) + zcom
				Xlj_new(i) = Xlj_old(i)

			end do

			if( lenion > 0 ) then

				do i = 1, lenion
					
					T(1) = Yion_old(i) - ycom
					T(2) = Zion_old(i) - zcom

					T = matmul( M, T )

					Yion_new(i) = T(1) + ycom
					Zion_new(i) = T(2) + zcom
					Xion_new(i) = Xion_old(i)
			
				end do
		
			end if

		case ( 2 )			! Rotation in x-z plane.

			do i = 1, lenlj
				
				T(1) = Zlj_old(i) - zcom
				T(2) = Xlj_old(i) - xcom

				T = matmul( M, T )

				Zlj_new(i) = T(1) + zcom
				Xlj_new(i) = T(2) + xcom
				Ylj_new(i) = Ylj_old(i)

			end do

			if( lenion > 0 ) then

				do i = 1, lenion

					T(1) = Zion_old(i) - zcom
					T(2) = Xion_old(i) - xcom

					T = matmul( M, T )

					Zion_new(i) = T(1) + zcom
					Xion_new(i) = T(2) + xcom
					Yion_new(i) = Yion_old(i)
			
				end do
		
			end if

		case ( 3 )			! Rotation in x-y plane.

			do i = 1, lenlj

				T(1) = Xlj_old(i) - xcom
				T(2) = Ylj_old(i) - ycom

				T = matmul( M, T )

				Xlj_new(i) = T(1) + xcom
				Ylj_new(i) = T(2) + ycom
				Zlj_new(i) = Zlj_old(i)

			end do

			if( lenion > 0 ) then

				do i = 1, lenion

					T(1) = Xion_old(i) - xcom
					T(2) = Yion_old(i) - ycom

					T = matmul( M, T )

					Xion_new(i) = T(1) + xcom
					Yion_new(i) = T(2) + ycom
					Zion_new(i) = Zion_old(i)
			
				end do
		
			end if
	
	end select

end if


do i = 1, lenlj
	
	if( Xlj_new(i) > BoxSize )  Xlj_new(i) = Xlj_new(i) - &
								BoxSize * aint( Xlj_new(i) / BoxSize )
	if( Ylj_new(i) > BoxSize )  Ylj_new(i) = Ylj_new(i) - &
								BoxSize * aint( Ylj_new(i) / BoxSize )
	if( Zlj_new(i) > BoxSize )  Zlj_new(i) = Zlj_new(i) - &
								BoxSize * aint( Zlj_new(i) / BoxSize )

	if( Xlj_new(i) < 0.0 )  Xlj_new(i) = Xlj_new(i) - &
							BoxSize * aint( Xlj_new(i) / BoxSize - 1 )
	if( Ylj_new(i) < 0.0 )  Ylj_new(i) = Ylj_new(i) - &
							BoxSize * aint( Ylj_new(i) / BoxSize - 1 )
	if( Zlj_new(i) < 0.0 )  Zlj_new(i) = Zlj_new(i) - &
							BoxSize * aint( Zlj_new(i) / BoxSize - 1 )

end do


if( stlj == 1 )	then

	if( Nmol(0) == 1 ) then

		ULJ_new = 0.0

	else

		call e6molecule( lenlj, Xlj_new, Ylj_new, Zlj_new, &
						 TYPElj(1:lenlj), DAMPlj2(1:lenlj), &
						 DAMPlj3(1:lenlj), &
						 Nlj - lenlj, Xlj(endlj+1:Nlj), &
						 Ylj(endlj+1:Nlj), Zlj(endlj+1:Nlj), &
						 TYPElj(endlj+1:Nlj), DAMPlj2(endlj+1:Nlj), &
						 DAMPlj3(endlj+1:Nlj), &
						 Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
						 BoxSize, ULJ_new	)

	end if

else if( stlj + lenlj - 1 == Nlj ) then

	call e6molecule( lenlj, Xlj_new, Ylj_new, Zlj_new, &
					 TYPElj(stlj:endlj), DAMPlj2(stlj:endlj), &
					 DAMPlj3(stlj:endlj), &
					 Nlj - lenlj, Xlj(1:stlj-1), Ylj(1:stlj-1), &
					 Zlj(1:stlj-1), TYPElj(1:stlj-1), &
					 DAMPlj2(1:stlj-1), DAMPlj3(1:stlj-1), &
					 Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					 BoxSize, ULJ_new )

else

	call e6molecule( lenlj, Xlj_new, Ylj_new, Zlj_new, &
					 TYPElj(stlj:endlj), DAMPlj2(stlj:endlj), &
					 DAMPlj3(stlj:endlj), &
					 Nlj - lenlj, temp1, temp2, temp3, temp4, &
					 temp20, temp21, &
					 Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					 BoxSize, ULJ_new )

end if

dULJ = ULJ_new - ULJ_old

if( lenion > 0 ) then

CoulCombo = ec * ec * 1.0e10 / ( 4.0 * Pi * eps0 * kB )

	do i = 1, lenion
	
		if( Xion_new(i) > BoxSize )	Xion_new(i) = Xion_new(i) - &
									BoxSize * aint( Xion_new(i) / BoxSize )
		if( Yion_new(i) > BoxSize ) Yion_new(i) = Yion_new(i) - &
									BoxSize * aint( Yion_new(i) / BoxSize )
		if( Zion_new(i) > BoxSize ) Zion_new(i) = Zion_new(i) - &
									BoxSize * aint( Zion_new(i) / BoxSize )

		if( Xion_new(i) < 0.0 ) Xion_new(i) = Xion_new(i) - &
								BoxSize * aint( Xion_new(i) / BoxSize - 1 )
		if( Yion_new(i) < 0.0 ) Yion_new(i) = Yion_new(i) - &
								BoxSize * aint( Yion_new(i) / BoxSize - 1 )
		if( Zion_new(i) < 0.0 ) Zion_new(i) = Zion_new(i) - &
								BoxSize * aint( Zion_new(i) / BoxSize - 1 )
	
	end do

	
	do i = 1, lenion
	
		if( Xion_old(i) > BoxSize ) Xion_old(i) = Xion_old(i) - &
									BoxSize * aint( Xion_old(i) / BoxSize )
		if( Yion_old(i) > BoxSize ) Yion_old(i) = Yion_old(i) - &
									BoxSize * aint( Yion_old(i) / BoxSize )
		if( Zion_old(i) > BoxSize ) Zion_old(i) = Zion_old(i) - &
									BoxSize * aint( Zion_old(i) / BoxSize )

		if( Xion_old(i) < 0.0 ) Xion_old(i) = Xion_old(i) - &
								BoxSize * aint( Xion_old(i) / BoxSize - 1 )
		if( Yion_old(i) < 0.0 ) Yion_old(i) = Yion_old(i) - &
								BoxSize * aint( Yion_old(i) / BoxSize - 1 )
		if( Zion_old(i) < 0.0 ) Zion_old(i) = Zion_old(i) - &
								BoxSize * aint( Zion_old(i) / BoxSize - 1 )
	
	end do

	if( stion == 1 ) then

		if( Nmol(0) == 1 ) then

			UREAL_old = 0.0

		else
	
			call RealMolecule( lenion, Xion_old, Yion_old, Zion_old, &
							   TYPEion(1:lenion), DAMPion(1:lenion), &
							   Nion - lenion, Xion(endion+1:Nion), &
							   Yion(endion+1:Nion), Zion(endion+1:Nion), &
							   TYPEion(endion+1:Nion), DAMPion(endion+1:Nion), &
							   Nham, Niongrs, CHARGE, &
							   BoxSize, Alpha, UREAL_old )

		end if

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

		call RealMolecule( lenion, Xion_old, Yion_old, Zion_old, &
						   TYPEion(stion:endion), DAMPion(stion:endion), &
						   Nion - lenion, Xion(1:stion-1), &
						   Yion(1:stion-1), Zion(1:stion-1), &
						   TYPEion(1:stion-1), DAMPion(1:stion-1), &
						   Nham, Niongrs, CHARGE, &
						   BoxSize, Alpha, UREAL_old )
	
	else

		temp5( 1:stion-1 ) = Xion( 1:stion-1 )
		temp5( stion:Nion-lenion ) = Xion( endion+1:Nion )
		
		temp6( 1:stion-1 ) = Yion( 1:stion-1 )
		temp6( stion:Nion-lenion ) = Yion( endion+1:Nion )
		
		temp7( 1:stion-1 ) = Zion( 1:stion-1 )
		temp7( stion:Nion-lenion ) = Zion( endion+1:Nion )
		
		temp8( 1:stion-1 ) = TYPEion( 1:stion-1 )
		temp8( stion:Nion-lenion ) = TYPEion( endion+1:Nion )

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

		call RealMolecule( lenion, Xion_old, Yion_old, Zion_old, &
						   TYPEion(stion:endion), DAMPion(stion:endion), &
						   Nion - lenion, temp5, temp6, temp7, temp8, temp22, &
						   Nham, Niongrs, CHARGE, BoxSize, Alpha, UREAL_old )

	end if
	
	if( stion == 1 ) then

		if( Nmol(0) == 1 ) then

			UREAL_new = 0.0

		else
	
			call RealMolecule( lenion, Xion_new, Yion_new, Zion_new, &
							   TYPEion(1:lenion), DAMPion(1:lenion), &
							   Nion - lenion, Xion(endion+1:Nion), &
							   Yion(endion+1:Nion), Zion(endion+1:Nion), &
							   TYPEion(endion+1:Nion), DAMPion(endion+1:Nion), &
							   Nham, Niongrs, CHARGE, &
							   BoxSize, Alpha, UREAL_new )

		end if

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

		call RealMolecule( lenion, Xion_new, Yion_new, Zion_new, &
						   TYPEion(stion:endion), DAMPion(stion:endion), &
						   Nion - lenion, Xion(1:stion-1), &
						   Yion(1:stion-1), Zion(1:stion-1), &
						   TYPEion(1:stion-1), DAMPion(1:stion-1), &
						   Nham, Niongrs, CHARGE, &
						   BoxSize, Alpha, UREAL_new )

	else
	
		call RealMolecule( lenion, Xion_new, Yion_new, Zion_new, &
						   TYPEion(stion:endion), DAMPion(stion:endion), &
						   Nion - lenion, temp5, temp6, temp7, temp8, temp22, &
						   Nham, Niongrs, CHARGE, BoxSize, Alpha, UREAL_new )

	end if

	dUREAL = UREAL_new - UREAL_old

	DELTAX = Xion_new - Xion_old
	DELTAY = Yion_new - Yion_old
	DELTAZ = Zion_new - Zion_old

	call Surf_Move( lenion, DELTAX, DELTAY, DELTAZ, &
					TYPEion(stion:endion), DAMPion(stion:endion), &
					Nham, Niongrs, CHARGE, BoxSize, &
					SUMQX, SUMQY, SUMQZ, &
					SUMQX_NEW, SUMQY_NEW, SUMQZ_NEW, dUSURF )

	call Fourier_Move( lenion, Xion_new, Yion_new, Zion_new, &
					   TYPEion(stion:endion), DAMPion(stion:endion), &
					   Nham, Niongrs, CHARGE, BoxSize, &
					   Kmax, Nkvec, KX, KY, KZ, CONST, &
					   EXPX(:,stion:endion), EXPY(:,stion:endion), &
					   EXPZ(:,stion:endion), EXPX_NEW(:,1:lenion), &
					   EXPY_NEW(:,1:lenion), EXPZ_NEW(:,1:lenion), &
					   SUMQEXPV, SUMQEXPV_NEW, dUFOURIER )

	dUREAL = dUREAL	* CoulCombo
	dUSURF = dUSURF	* CoulCombo
	dUFOURIER = dUFOURIER * CoulCombo

else

	dUREAL = 0.0
	dUSURF = 0.0
	dUFOURIER = 0.0

end if

dU = dULJ + dUREAL + dUSURF + dUFOURIER

LNPSI_new = LNPSI - BETA * dU

Largest = maxval( LNW + LNPSI_new )

LnPi_new = log( sum( exp( LNW + LNPSI_new - Largest ) ) ) + Largest

if( log( ran2(Seed) ) < LnPi_new - LnPi ) then		

	Success = .True.

	LnPi = LnPi_new 

	LNPSI = LNPSI_new

	Xlj( stlj:endlj ) = Xlj_new( 1:lenlj )
	Ylj( stlj:endlj ) = Ylj_new( 1:lenlj )
	Zlj( stlj:endlj ) = Zlj_new( 1:lenlj )

	ULJ = ULJ + dULJ

	if( lenion > 0 ) then

		Xion( stion:endion ) = Xion_new( 1:lenion )
		Yion( stion:endion ) = Yion_new( 1:lenion )
		Zion( stion:endion ) = Zion_new( 1:lenion )

		EXPX(:,stion:endion) = EXPX_NEW(:,1:lenion)
		EXPY(:,stion:endion) = EXPY_NEW(:,1:lenion)
		EXPZ(:,stion:endion) = EXPZ_NEW(:,1:lenion)

		SUMQEXPV = SUMQEXPV_NEW

		SUMQX = SUMQX_NEW
		SUMQY = SUMQY_NEW
		SUMQZ = SUMQZ_NEW

		UFOURIER = UFOURIER + dUFOURIER
		UREAL = UREAL + dUREAL
		USURF = USURF + dUSURF

	end if

end if

deallocate( Xlj_old )
deallocate( Ylj_old )
deallocate( Zlj_old )

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

if( lenion > 0 ) then

	deallocate( Xion_old )
	deallocate( Yion_old )
	deallocate( Zion_old )

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

	deallocate( DELTAX )
	deallocate( DELTAY )
	deallocate( DELTAZ )

	deallocate( EXPX_NEW )
	deallocate( EXPY_NEW )
	deallocate( EXPZ_NEW )

end if

return

end subroutine Disp_Rot