resupdate.f90

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

subroutine ResUpdate( reslen, Xb, Yb, Zb, &
					  nresmv, maxresmv, maxresp, &
					  PROB_MV, RESMOVE, RESPARAM, &
					  nvib, nbend, ntor, &
					  maxvib, maxbend, maxtor, &
					  IVIB, IBEND, ITOR, &
					  RVIB, RBEND, RTOR, &
					  RES_ATT, RES_SUC, &
					  MaxBeads, &
					  BEADTYPE, GROUPTYPE, &
					  BONDSAPART, Nham, Nljgrs, &
					  EPS, SIG, CP, ALP, RMAX, &
					  Niongrs, CHARGE, &
					  f14_lj, f14_ion, &
					  BETA, LNW, BoxSize, &
					  Uintra, Ulj, Uion, &
					  Seed, th_h, phi_h )

implicit none

integer, intent(in)						:: reslen

real, dimension(reslen)					:: Xb, Yb, Zb

integer, intent(in)						:: nresmv, maxresmv, maxresp

real, dimension(maxresmv), intent(in)	:: PROB_MV

character*10, dimension(maxresmv)		:: RESMOVE

integer, dimension(maxresp, maxresmv)	:: RESPARAM

integer, intent(in)						:: nvib, nbend, ntor
integer, intent(in)						:: maxvib, maxbend, maxtor

integer, dimension(2, maxvib)			:: IVIB
integer, dimension(3, maxbend)			:: IBEND
integer, dimension(4, maxtor)			:: ITOR

real, dimension(2, maxvib)				:: RVIB
real, dimension(2, maxbend)				:: RBEND
real, dimension(4, maxtor)				:: RTOR

integer, dimension(5)					:: RES_ATT, RES_SUC

integer, intent(in)						:: MaxBeads
character*5, dimension(MaxBeads)		:: BEADTYPE
integer, dimension(MaxBeads)			:: GROUPTYPE 
integer, dimension(MaxBeads,MaxBeads)	:: BONDSAPART

integer, intent(in)						:: Nham
integer, intent(in)						:: Nljgrs

real, dimension(Nljgrs,Nljgrs,Nham)		:: EPS
real, dimension(Nljgrs,Nljgrs,Nham)		:: SIG
real, dimension(Nljgrs,Nljgrs,Nham)		:: CP
real, dimension(Nljgrs,Nljgrs,Nham)		:: ALP
real, dimension(Nljgrs,Nljgrs,Nham)		:: RMAX

integer, intent(in)						:: Niongrs
real, dimension(Niongrs,Nham)			:: CHARGE

real, intent(in)						:: f14_lj
real, intent(in)						:: f14_ion

real, dimension(Nham)					:: BETA

real, dimension(Nham)					:: LNW

! BoxSize is the length of the simulation box.

real, intent(in)						:: BoxSize

real									:: Uintra

real, dimension(Nham)					:: Ulj, Uion

integer									:: Seed

integer, dimension(0:180)				:: th_h
integer, dimension(-180:180)			:: phi_h

real, external							:: ran2

! Local stuff

integer							:: i, j, k, mv
integer							:: axis

integer, dimension(0:20)		:: mapfcs

logical							:: success

real, parameter					:: Pi = 3.14159265359

real							:: r, dUintra
real							:: xcom, ycom, zcom
real							:: xmin, ymin, zmin
real							:: dtheta
real, dimension(2)				:: T
real, dimension(2,2)			:: M
real, dimension(1)				:: ZERO2 = 0.0
real, dimension(Nham)			:: dUlj, dUion


r = ran2(Seed)
mv = 0
i = 1

do while ( mv == 0 )

	if( r < PROB_MV(i) ) mv = i

	i = i + 1

end do

call Outfold( reslen, 0, BoxSize, Xb, Yb, Zb, &
			  ZERO2, ZERO2, ZERO2 )

select case ( RESMOVE(mv) )

	case( 'CrankSh' )

		do i = 1, RESPARAM(1,mv) + 2

			mapfcs(i) = RESPARAM(i+2,mv)

		end do

		mapfcs(0) = RESPARAM(2,mv)

		call CrankSh( reslen, Xb, Yb, Zb, &
						nvib, nbend, ntor, &
						maxvib, maxbend, maxtor, &
						IVIB, IBEND, ITOR, &
						RVIB, RBEND, RTOR, &
						MaxBeads, &
						BEADTYPE, GROUPTYPE, &
						BONDSAPART, Nham, Nljgrs, &
						EPS, SIG, CP, ALP, RMAX, &
						Niongrs, CHARGE, &
						f14_lj, f14_ion, &
						RESPARAM(1,mv), mapfcs, &
						BETA, LNW, BoxSize, &
						dUintra, dUlj, dUion, &
						success, Seed, th_h, phi_h )

		Uintra = Uintra + dUintra
		Ulj = Ulj + dUlj
		Uion = Uion + dUion

		RES_ATT(1) = RES_ATT(1) + 1

		if(success) RES_SUC(1) = RES_SUC(1) + 1

	case( 'BendSh' )

		do i = 1, RESPARAM(1,mv) + 1

			mapfcs(i) = RESPARAM(i+2,mv)

		end do

		mapfcs(0) = RESPARAM(2,mv)

		call BendSh( reslen, Xb, Yb, Zb, &
					 nvib, nbend, ntor, &
					 maxvib, maxbend, maxtor, &
					 IVIB, IBEND, ITOR, &
					 RVIB, RBEND, RTOR, &
					 MaxBeads, &
					 BEADTYPE, GROUPTYPE, &
					 BONDSAPART, Nham, Nljgrs, &
					 EPS, SIG, CP, ALP, RMAX, &
					 Niongrs, CHARGE, &
					 f14_lj, f14_ion, &
					 RESPARAM(1,mv), mapfcs, &
					 BETA, LNW, BoxSize, &
					 dUintra, dUlj, dUion, &
					 success, Seed, th_h, phi_h )

		Uintra = Uintra + dUintra
		Ulj = Ulj + dUlj
		Uion = Uion + dUion

		RES_ATT(2) = RES_ATT(2) + 1

		if(success) RES_SUC(2) = RES_SUC(2) + 1

end select


k = int( ran2(Seed) * reslen ) + 1

xcom = Xb(k)
ycom = Yb(k)
zcom = Zb(k)

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 k = 1, reslen
			
			T(1) = Yb(k) - ycom
			T(2) = Zb(k) - zcom

			T = matmul( M, T )

			Yb(k) = T(1) + ycom
			Zb(k) = T(2) + zcom
			Xb(k) = Xb(k)

		end do

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

		do k = 1, reslen
			
			T(1) = Zb(k) - zcom
			T(2) = Xb(k) - xcom

			T = matmul( M, T )

			Zb(k) = T(1) + zcom
			Xb(k) = T(2) + xcom
			Yb(k) = Yb(k)

		end do

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

		do k = 1, reslen

			T(1) = Xb(k) - xcom
			T(2) = Yb(k) - ycom

			T = matmul( M, T )

			Xb(k) = T(1) + xcom
			Yb(k) = T(2) + ycom
			Zb(k) = Zb(k)

		end do

end select

xmin = minval( Xb )
ymin = minval( Yb )
zmin = minval( Zb )

do j = 1, reslen

	Xb(j) = Xb(j) - xmin + 1.0
	Yb(j) = Yb(j) - ymin + 1.0
	Zb(j) = Zb(j) - zmin + 1.0

end do

return

end subroutine ResUpdate