cranksh.f90

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

! 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:1+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, c, d
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, u12, nrm, vv, ww
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, 4

!	r(i,1) = Xb_new(i)
!	r(i,2) = Yb_new(i)
!	r(i,3) = Zb_new(i)

!end do

!u(1,:) = r(2,:) - r(1,:)
!u(2,:) = r(3,:) - r(1,:)
!u(3,:) = r(4,:) - r(1,:)

!do i = 1, 3

!	l(i) = sqrt(dot_product(u(i,:),u(i,:)))

!	j = nint( l(i) * 80.0 )

!	th_h(j) = th_h(j) + 1

!	u(i,:) = u(i,:) / l(i)

!end do

! Calculate the old bond angles.

!th_m(1) = dot_product(u(1,:),u(2,:))
!j = int( ( th_m(1) + 1.0 ) / 2.0 * 180.0 )
!th_h(j) = th_h(j) + 1

!th_m(1) = dot_product(u(1,:),u(2,:))
!j = int( ( th_m(1) + 1.0 ) / 2.0 * 60.0 )
!th_h(j) = th_h(j) + 1

!th_m(2) = dot_product(u(1,:),u(3,:))
!j = 60 + int( ( th_m(2) + 1.0 ) / 2.0 * 60.0 )
!th_h(j) = th_h(j) + 1

!th_m(3) = dot_product(u(2,:),u(3,:))
!j = 120 + int( ( th_m(3) + 1.0 ) / 2.0 * 60.0 )
!th_h(j) = th_h(j) + 1


!do i = 1, 4

!	r(i,1) = Xb_new(i)
!	r(i,2) = Yb_new(i)
!	r(i,3) = Zb_new(i)

!end do

!u(1,:) = r(3,:) - r(2,:)
!u(2,:) = r(4,:) - r(3,:)
!u(3,:) = r(2,:) - r(4,:)

!u(1,:) = r(2,:) - r(1,:)
!u(2,:) = r(3,:) - r(1,:)
!u(3,:) = r(4,:) - r(1,:)

!do i = 1, 3

!	l(i) = sqrt(dot_product(u(i,:),u(i,:)))

!	j = nint( l(i) * 80.0 )

!	th_h(j) = th_h(j) + 1

!	u(i,:) = u(i,:) / l(i)

!end do

! Calculate the old bond angles.

!th_m(1) = dot_product(u(1,:),u(2,:))
!th_m(2) = dot_product(u(1,:),u(3,:))
!th_m(3) = dot_product(u(2,:),u(3,:))

!do i = 1, 0
		  
!	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))

call Sphere( 0.0, 0.0, 0.0, 1.0, nrm(1), nrm(2), nrm(3), Seed )

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

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

do i = 1, Nrot

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

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

	bl = sqrt( dot_product( u12, u12 ) )

	u12 = u12 / bl

	ctheta = dot_product( nrm, u12 )

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

	stheta = sqrt( dot_product( vv, vv ) )

	vv = vv / stheta
	
	ww(1) = vv(2) * u12(3) - vv(3) * u12(2)
	ww(2) = vv(3) * u12(1) - vv(1) * u12(3)
	ww(3) = vv(1) * u12(2) - vv(2) * u12(1)

	UU = cphi * u12 + sphi * ww

	Xb_new(mapf(i+1)) = Xc(1) + bl * UU(1)
	Yb_new(mapf(i+1)) = Yc(1) + bl * UU(2)
	Zb_new(mapf(i+1)) = Zc(1) + bl * UU(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+1)

	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 BendSh