cranksh.f90

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

Subroutine CrankSh( Nb, 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, &
					Nrot, mapf, BETA, LNW, BoxSize, &
					dUintra, dUlj, dUion, &
					success, Seed, th_h, phi_h )

implicit none

! Nbs is the number of beads in the molecule.

integer, intent(in)										:: Nb

! Xb, Yb, and Zb are the coordinates of the beads.

real, dimension(Nb), intent(inout)						:: Xb, Yb, Zb

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, 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

integer													:: Nrot
integer, dimension(0:2+Nrot)							:: mapf

real, dimension(Nham)									:: BETA

real, dimension(Nham)									:: LNW

! BoxSize is the length of the simulation box.

real, intent(in)										:: BoxSize

real													:: dUintra

real, dimension(Nham)									:: dUlj, dUion

logical													:: success

! Seed is the current random number generator seed value.

integer, intent(inout)									:: Seed

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

real, external											:: ran2

! Local Variables

integer										:: i, j
integer										:: mapf1

logical, dimension(Nbend)					:: cbend
logical, dimension(Ntor)					:: ctor

real, parameter								:: Pi = 3.14159265359
real, parameter								:: ec = 1.60217733e-19
real, parameter								:: eps0 = 8.854187817e-12
real, parameter								:: kB = 1.380658e-23

real										:: dphi_max
real										:: dphi
real										:: sth, a, b
real										:: bl, cphi, sphi, h
real										:: ctheta, stheta
real										:: Ubend_m, Ubend_n
real										:: Utor_m, Utor_n
real										:: PiRatio, Largest
real										:: CoulCombo
											   
real, dimension(1)							:: ZERO2 = 0.0
real, dimension(1)							:: ONE = 1.0
real, dimension(3)							:: w, y, z
real, dimension(3)							:: UU, TT, u21, u23, nrm, vv
real, dimension(4)							:: Xc, Yc, Zc
real, dimension(8)							:: l, th_m, phi
real, dimension(8,3)						:: r, u
real, dimension(Nb)							:: Xb_new, Yb_new, Zb_new
real, dimension(Nb)							:: Xb_old, Yb_old, Zb_old
real, dimension(Nham)						:: Ulj_new, Ulj_old
real, dimension(Nham)						:: Uion_new, Uion_old
real, dimension(Nham)						:: dU, U_part


										
success = .false.

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

dUintra = 0.0
dUlj = 0.0
dUion = 0.0

! The molecule must be unfolded from periodic boundary conditions.

Xb_new = Xb
Yb_new = Yb
Zb_new = Zb

call Outfold( Nb, 0, BoxSize, Xb_new, Yb_new, Zb_new, &
			  ZERO2, ZERO2, ZERO2 )

Xb_old = Xb_new
Yb_old = Yb_new
Zb_old = Zb_new

! For data collection about bending and torsion angles.

!do i = 1, 6
!
!	r(i,1) = Xb_new(i)
!	r(i,2) = Yb_new(i)
!	r(i,3) = Zb_new(i)
!
!end do
!
!do i = 2, 6
!
!	u(i,:) = r(i,:) - r(i-1,:)
!
!	l(i) = sqrt(dot_product(u(i,:),u(i,:)))
!
!	u(i,:) = u(i,:) / l(i)
!
!end do
!
!! Calculate the old bond angles.
!
!do i = 2, 5
!		  
!	th_m(i) = dot_product(u(i,:),u(i+1,:))
!
!	if( th_m(i) > 1.0 ) then
!
!		write(*,*) ' th_m = ', th_m(i)
!		th_m(i) = 1.0
!
!	else if( th_m(i) < -1.0 ) then
!
!		write(*,*) ' th_m = ', th_m(i)
!		th_m(i) = -1.0
!
!	end if
!
!	th_m(i) = dacos(th_m(i))
!
!	j = nint( th_m(i) * 180.0 / pi )
!
!	th_h(j) = th_h(j) + 1
!
!end do
!
!! find phi0_m, phi1_m
!! cos(phi(i)) = ( (u(i-1) x u(i)) / sin(th(i-1)) ) dot ( (u(i) x u(i+1)) / sin(th(i)) )
!
!do i = 3, 5
!
!	sth = dsin(th_m(i-1))
!
!	z(1) = ( u(i-1,2) * u(i,3) - u(i-1,3) * u(i,2) ) / sth 
!	z(2) = ( u(i-1,3) * u(i,1) - u(i-1,1) * u(i,3) ) / sth 
!	z(3) = ( u(i-1,1) * u(i,2) - u(i-1,2) * u(i,1) ) / sth 
!
!	sth = dsin(th_m(i))
!
!	y(1) = ( u(i,2) * u(i+1,3) - u(i,3) * u(i+1,2) ) / sth 
!	y(2) = ( u(i,3) * u(i+1,1) - u(i,1) * u(i+1,3) ) / sth 
!	y(3) = ( u(i,1) * u(i+1,2) - u(i,2) * u(i+1,1) ) / sth
!
!	a = -dot_product( z, y )	! a = cos(phi)
!
!	w(1) = z(2) * y(3) - z(3) * y(2) 
!	w(2) = z(3) * y(1) - z(1) * y(3) 
!	w(3) = z(1) * y(2) - z(2) * y(1)
!	
!	b = -dot_product( w(:), u(i,:) )	! b = sin(phi)
!	
!	if( a > 1.0 ) a = 1.0
!	if( a < -1.0 ) a = -1.0
!
!	if( a > 0 .AND. b > 0 ) then
!	
!		phi(i) = dacos(a)
!		  
!	else if( a > 0 .AND. b < 0 ) then
!	
!		phi(i) = -dacos(a)
!		  
!	else if( a < 0 .AND. b > 0 ) then
!	
!		phi(i) = dacos(a)
!		  
!	else if( a < 0 .AND. b < 0 ) then
!	
!		phi(i) = -dacos(a)
!
!	end if
!
!	j = nint( phi(i) * 180.0 / pi )
!
!	phi_h(j) = phi_h(j) + 1
!
!end do

dphi_max = real( mapf(0) ) / 180.0 * Pi

Xc(1) = Xb_new(mapf(1))
Yc(1) = Yb_new(mapf(1))
Zc(1) = Zb_new(mapf(1))

Xc(2) = Xb_new(mapf(2))
Yc(2) = Yb_new(mapf(2))
Zc(2) = Zb_new(mapf(2))

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

bl = sqrt( dot_product( u21, u21 ) )

u21 = u21 / bl

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

cphi = cos( dphi )
sphi = sin( dphi )

do i = 1, Nrot

	Xc(3) = Xb_new(mapf(i+2))
	Yc(3) = Yb_new(mapf(i+2))
	Zc(3) = Zb_new(mapf(i+2))

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

	bl = sqrt( dot_product( u23, u23 ) )

	u23 = u23 / bl

	ctheta = dot_product( u21, u23 )

	vv = u23 - ctheta * u21

	stheta = sqrt( dot_product( vv, vv ) )

	vv = vv / stheta

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

	nrm = -nrm / sqrt( dot_product( nrm, nrm ) )

	h = ctheta / stheta

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

	UU = cphi * vv + sphi * nrm

	TT(1) = ctheta * u21(1) + stheta * UU(1)
	TT(2) = ctheta * u21(2) + stheta * UU(2)
	TT(3) = ctheta * u21(3) + stheta * UU(3)

!	TT = TT / sqrt( dot_product( TT, TT ) )

	Xb_new(mapf(i+2)) = Xc(2) + bl * TT(1)
	Yb_new(mapf(i+2)) = Yc(2) + bl * TT(2)
	Zb_new(mapf(i+2)) = Zc(2) + bl * TT(3)

end do

Ubend_m = 0.0
Ubend_n = 0.0

Utor_m = 0.0
Utor_n = 0.0

Ulj_new = 0.0
Ulj_old = 0.0

Uion_new = 0.0
Uion_old = 0.0

cbend = .false.
ctor = .false.

do j = 1, Nrot

	mapf1 = mapf(j+2)

	do i = 1, Nbend

		if( IBEND(1,i) == mapf1 .OR. IBEND(2,i) == mapf1 .OR. &
			IBEND(3,i) == mapf1 )	then

			if( cbend(i) ) cycle

			Xc(1:3) = Xb_new( IBEND(1:3,i) )
			Yc(1:3) = Yb_new( IBEND(1:3,i) )
			Zc(1:3) = Zb_new( IBEND(1:3,i) )

			call IntraBend( Xc, Yc, Zc, RBEND(1,i), RBEND(2,i), b )

			Ubend_n = Ubend_n + b

			Xc(1:3) = Xb_old( IBEND(1:3,i) )
			Yc(1:3) = Yb_old( IBEND(1:3,i) )
			Zc(1:3) = Zb_old( IBEND(1:3,i) )

			call IntraBend( Xc, Yc, Zc, RBEND(1,i), RBEND(2,i), b )

			Ubend_m = Ubend_m + b

			cbend(i) = .true.

		end if

	end do			

	do i = 1, Ntor

		if( ITOR(1,i) == mapf1 .OR. ITOR(2,i) == mapf1 .OR. &
			ITOR(3,i) == mapf1 .OR. ITOR(4,i) == mapf1 ) then

			if( ctor(i) ) cycle

			Xc(1:4) = Xb_new( ITOR(1:4,i) )
			Yc(1:4) = Yb_new( ITOR(1:4,i) )
			Zc(1:4) = Zb_new( ITOR(1:4,i) )

			call IntraTorsion( Xc, Yc, Zc, 2, RTOR(1:4,i), b )

			Utor_n = Utor_n + b

			Xc(1:4) = Xb_old( ITOR(1:4,i) )
			Yc(1:4) = Yb_old( ITOR(1:4,i) )
			Zc(1:4) = Zb_old( ITOR(1:4,i) )

			call IntraTorsion( Xc, Yc, Zc, 2, RTOR(1:4,i), b )

			Utor_m = Utor_m + b

			ctor(i) = .true.

		end if

	end do
	
! If the intramolecular nonbonded interactions are included
! as part of the reservoir energy, then uncomment the 
! following lines.  Also change the appropriate lines in 
! resread.f90, alpha_ch.f90 and grow.f90. 

!	if( BEADTYPE( mapf1 ) == 'LJ' ) then
!
!		do i = 1, Nb
!
!			if( BEADTYPE(i) == 'LJ' .AND. BONDSAPART(mapf1,i) >= 3 ) then
!
!				call e6molecule( 1, Xb_new(mapf1), Yb_new(mapf1), &
!								 Zb_new(mapf1), GROUPTYPE(mapf1), &
!								 ONE, ONE, &
!								 1, Xb_new(i), Yb_new(i), &
!								 Zb_new(i), GROUPTYPE(i), &
!								 ONE, ONE, &
!								 Nham, Nljgrs, EPS, SIG, CP, &
!								 ALP, RMAX, BoxSize, U_part )
!
!				if( BONDSAPART(mapf1,i) == 3 ) U_part = f14_lj * U_part 
!
!				Ulj_new = Ulj_new + U_part
!
!				call e6molecule( 1, Xb_old(mapf1), Yb_old(mapf1), &
!								 Zb_old(mapf1), GROUPTYPE(mapf1), &
!								 ONE, ONE, &
!								 1, Xb_old(i), Yb_old(i), &
!								 Zb_old(i), GROUPTYPE(i), &
!								 ONE, ONE, &
!								 Nham, Nljgrs, EPS, SIG, CP, &
!								 ALP, RMAX, BoxSize, U_part )
!
!				if( BONDSAPART(mapf1,i) == 3 ) U_part = f14_lj * U_part 
!
!				Ulj_old = Ulj_old + U_part
!
!			end if
!
!		end do
!
!	end if
!
!	if( BEADTYPE( mapf1 ) == 'ION' ) then
!
!		do i = 1, Nb
!
!			if( BEADTYPE(i) == 'ION' .AND. BONDSAPART(mapf1,i) >= 3 ) then
!
!				call IonMolecule( 1, Xb_new(mapf1), Yb_new(mapf1), &
!								  Zb_new(mapf1), GROUPTYPE(mapf1), ONE, &
!								  1, Xb_new(i), Yb_new(i), &
!								  Zb_new(i), GROUPTYPE(i), ONE, &
!								  Nham, Niongrs, CHARGE, &
!								  BoxSize, U_part )
!
!				if( BONDSAPART(mapf1,i) == 3 ) U_part = f14_ion * U_part 
!
!				Uion_new = Uion_new + U_part * CoulCombo
!
!				call IonMolecule( 1, Xb_old(mapf1), Yb_old(mapf1), &
!								  Zb_old(mapf1), GROUPTYPE(mapf1), ONE, &
!								  1, Xb_old(i), Yb_old(i), &
!								  Zb_old(i), GROUPTYPE(i), ONE, &
!								  Nham, Niongrs, CHARGE, &
!								  BoxSize, U_part )
!
!				if( BONDSAPART(mapf1,i) == 3 ) U_part = f14_ion * U_part 
!
!				Uion_old = Uion_old + U_part * CoulCombo
!
!			end if
!
!		end do
!
!	end if

end do

dU = Ubend_n - Ubend_m + Utor_n - Utor_m

!dU = dU + Ulj_new - Ulj_old + Uion_new - Uion_old

dU = -BETA * dU

Largest = maxval( LNW + dU )

PiRatio = log( sum( exp( LNW + dU - Largest ) ) ) + Largest

if( log(ran2(Seed)) < PiRatio ) then

	success = .true.

	dUintra = Ubend_n - Ubend_m + Utor_n - Utor_m
	dUlj = Ulj_new - Ulj_old 
	dUion = Uion_new - Uion_old

	Xb = Xb_new
	Yb = Yb_new
	Zb = Zb_new

end if

return

end subroutine CrankSh


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

Subroutine BendSh( Nb, 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, &
				   Nrot, mapf, BETA, LNW, BoxSize, &
				   dUintra, dUlj, dUion, &
				   success, Seed, th_h, phi_h )

implicit none