resmethods.f90

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

subroutine ResRand( reslen, Xres, Yres, Zres, &
					BoxSize, Seed )

implicit none

integer									:: reslen

real, dimension(reslen)					:: Xres, Yres, Zres

real									:: BoxSize

integer									:: Seed

real, external							:: ran2

! Local Stuff

integer									:: i

real									:: xtrans, ytrans, ztrans
real, dimension(1)						:: ZERO2 = 0.0


call Outfold( reslen, 0, BoxSize, Xres, Yres, Zres, &
			  ZERO2, ZERO2, ZERO2 )

xtrans = ran2(Seed) * BoxSize
ytrans = ran2(Seed) * BoxSize
ztrans = ran2(Seed) * BoxSize

xtrans = xtrans - Xres(1)
ytrans = ytrans - Yres(1)
ztrans = ztrans - Zres(1)

do i = 1, reslen

	Xres(i) = Xres(i) + xtrans
	Yres(i) = Yres(i) + ytrans
	Zres(i) = Zres(i) + ztrans

end do


end subroutine ResRand

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

subroutine ResRot( reslen, Xres, Yres, Zres, &
					BoxSize, Seed )

implicit none

integer									:: reslen

real, dimension(reslen)					:: Xres, Yres, Zres

real									:: BoxSize

integer									:: Seed

real, external							:: ran2

! Local Stuff

integer									:: i, axis

real, parameter							:: Pi = 3.14159265359

real									:: xcom, ycom, zcom
real									:: dtheta
real, dimension(1)						:: ZERO2 = 0.0
real, dimension(2)						:: T
real, dimension(2,2)					:: M


call Outfold( reslen, 0, BoxSize, Xres, Yres, Zres, &
			  ZERO2, ZERO2, ZERO2 )

xcom = Xres(1)
ycom = Yres(1)
zcom = Zres(1)

dtheta = ( 2.0 * ran2(Seed) - 1.0 ) * Pi

M(1,1) =  cos( dtheta )
M(2,1) = -sin( dtheta )
M(1,2) =  sin( dtheta )
M(2,2) =  cos( dtheta )

axis = int( 3.0 * ran2(Seed) ) + 1

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

		do i = 1, reslen
			
			T(1) = Yres(i) - ycom
			T(2) = Zres(i) - zcom

			T = matmul( M, T )

			Yres(i) = T(1) + ycom
			Zres(i) = T(2) + zcom
			Xres(i) = Xres(i)

		end do

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

		do i = 1, reslen
			
			T(1) = Zres(i) - zcom
			T(2) = Xres(i) - xcom

			T = matmul( M, T )

			Zres(i) = T(1) + zcom
			Xres(i) = T(2) + xcom
			Yres(i) = Yres(i)

		end do

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

		do i = 1, reslen

			T(1) = Xres(i) - xcom
			T(2) = Yres(i) - ycom

			T = matmul( M, T )

			Xres(i) = T(1) + xcom
			Yres(i) = T(2) + ycom
			Zres(i) = Zres(i)

		end do

end select

return

end subroutine ResRot

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

subroutine ResCir( reslen, Xres, Yres, Zres, &
					BoxSize, Seed )

implicit none

integer									:: reslen

real, dimension(reslen)					:: Xres, Yres, Zres

real									:: BoxSize

integer									:: Seed

real, external							:: ran2

! Local Stuff

integer										:: i

real, parameter								:: Pi = 3.14159265359
real										:: dphi
real										:: bl21, bl23, c, s
real										:: ctheta, stheta

real, dimension(1)							:: ZERO2 = 0.0
real, dimension(3)							:: UU, TT, VV, u21, u23, nrm
real, dimension(3)							:: Xc, Yc, Zc


call Outfold( reslen, 0, BoxSize, Xres, Yres, Zres, &
			  ZERO2, ZERO2, ZERO2 )

Xc(1) = Xres(1)
Yc(1) = Yres(1)
Zc(1) = Zres(1)

Xc(2) = Xres(2)
Yc(2) = Yres(2)
Zc(2) = Zres(2)

Xc(3) = Xres(3)
Yc(3) = Yres(3)
Zc(3) = Zres(3)

u21(1) = Xc(1) - Xc(2)
u21(2) = Yc(1) - Yc(2)
u21(3) = Zc(1) - Zc(2)

bl21 = sqrt( dot_product( u21, u21 ) )

u21 = u21 / bl21

u23(1) = Xc(3) - Xc(2)
u23(2) = Yc(3) - Yc(2)
u23(3) = Zc(3) - Zc(2)

bl23 = sqrt( dot_product( u23, u23 ) )

u23 = u23 / bl23

nrm(1) = u21(2) * u23(3) - u21(3) * u23(2)
nrm(2) = u21(3) * u23(1) - u21(1) * u23(3)
nrm(3) = u21(1) * u23(2) - u21(2) * u23(1)

stheta = sqrt( dot_product( nrm, nrm ) )

nrm = nrm / stheta

dphi = 2.0 * ( ran2(Seed) - 0.5 ) * pi

c = cos( dphi )
s = sin( dphi )

! dphi = pi / 2 takes bead to original position.

do i = 3, reslen

	Xc(3) = Xres(i)
	Yc(3) = Yres(i)
	Zc(3) = Zres(i)

	u23(1) = Xc(3) - Xc(2)
	u23(2) = Yc(3) - Yc(2)
	u23(3) = Zc(3) - Zc(2)

	bl23 = dot_product( u23, u23 )

	if( bl23 > -1.0e-7 .AND. bl23 < +1.0e-7 ) bl23 = 0.0

	bl23 = sqrt( bl23 )

	if( bl23 > 0.0 ) then

		u23 = u23 / bl23

		ctheta = dot_product( u21, u23 )

		if( ctheta < 1.0 + 1.0e-7 .AND. ctheta > 1.0 - 1.0e-7 ) then

			stheta = 0.0

			UU = 0.0

		else

			VV(1) = u21(2) * u23(3) - u21(3) * u23(2)
			VV(2) = u21(3) * u23(1) - u21(1) * u23(3)
			VV(3) = u21(1) * u23(2) - u21(2) * u23(1)

			stheta = sqrt( dot_product( VV, VV ) )

			VV = VV / stheta

			UU(1) = c * VV(1) + s * u21(3) * VV(2) - s * u21(2) * VV(3)
			UU(2) = c * VV(2) + s * u21(1) * VV(3) - s * u21(3) * VV(1)		
			UU(3) = c * VV(3) + s * u21(2) * VV(1) - s * u21(1) * VV(2)

		end if

		TT = ctheta * u21 + stheta * UU

	end if

	Xres(i)	= Xc(2) + bl23 * TT(1)
	Yres(i)	= Yc(2) + bl23 * TT(2)
	Zres(i)	= Zc(2) + bl23 * TT(3)

end do

return

end subroutine ResCir