conrot.f90

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

				end if

				! Routine to find solution(s).

				gate_tmp = 0

				do j = 1, nb

					nsol = nsol + 1

					a = phi(1)

					if( j == nb ) gate_tmp = gate

					phi_sol(1,nsol) = zbrent4(F3,xb1(j),xb2(j),1.0e-7, &
												i, l(1:4), th(1:3), phi(1:3), &
												q1, q2, qp1, r41, &
												gate_tmp, T1, T2 )

					phi(1) = a

					if( nint( phi_sol(1,nsol) ) == 1000 ) then
					
						nsol = nsol - 1
						cycle
					
					end if
					
					phi_sol(0,nsol) = phi(0)
					phi_sol(2,nsol) = phi(2)
					phi_sol(3,nsol) = phi(3)

					! determine r_sol's

					MM = matmul( matmul( T0lab, T0 ), T1 )
					y = 0.0
					y(1) = l(2)
					r_sol(2,:,nsol) = r_sol(1,:,nsol) + matmul( MM(:,:), y(:) )

					MM = matmul( MM, T2 )
					y = 0.0
					y(1) = l(3)
					r_sol(3,:,nsol) = r_sol(2,:,nsol) + matmul( MM(:,:), y(:) )

				end do							

			end if

			phi(1) = phi(1) + phi_w

			if( F3a2_sol ) phi(1) = phi(1) - phi_inc

			F3a1_num = F3a2_num
			F3a1_sol = F3a2_sol
			dF31 = dF32

		end do

	end do

	!if( mod(nsol,2) > 0 ) then

	!	write(*,*) ' Number of solutions in ConRot is Odd!! '

	!end if

!	if( nsol == 0 ) cycle
	if( nsol == 0 ) return

	!return

	Ubt_m = 0.0
	Ubt_n = 0.0

	Ubend_m = 0.0
	Ubend_n = 0.0

	Utor_m = 0.0
	Utor_n = 0.0

	do j = 1, nsol

		do i = 1, 3

			Xb_new(mapf(i)) = r_sol(i,1,j)
			Yb_new(mapf(i)) = r_sol(i,2,j) 
			Zb_new(mapf(i)) = r_sol(i,3,j)

		end do

		do i = 1, Nbend

			ibd = 0

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

				ibd = i

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

				ibd = i

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

				ibd = i

			end if

			if( ibd == i ) then

				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(j) = Ubend_n(j) + b

				if( IBEND(1,i) == mapf(3) .AND. IBEND(2,i) == mapf(-2) .AND. &
					IBEND(3,i) == mapf(-1)	) then

					Ubt_n(j) = Ubt_n(j) + b

				else if( IBEND(3,i) == mapf(3) .AND. IBEND(2,i) == mapf(-2) .AND. &
					IBEND(1,i) == mapf(-1)	) then

					Ubt_n(j) = Ubt_n(j) + b

				end if

			end if

		end do			

		do i = 1, Ntor

			ibd = 0

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

				ibd = i

			else if( ITOR(2,i) == mapf(1) .OR. ITOR(2,i) == mapf(2) .OR. &
						ITOR(2,i) == mapf(3) )	then

				ibd = i

			else if( ITOR(3,i) == mapf(1) .OR. ITOR(3,i) == mapf(2) .OR. &
						ITOR(3,i) == mapf(3) )	then

				ibd = i

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

				ibd = i

			end if

			if( ibd == i ) then

				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(j) = Utor_n(j) + b

			end if

		end do			

	end do

	boltz(:) = exp( -beta * ( Ubt_n(:) + Utor_n(:) ) )
!	boltz(:) = exp( -beta * Utor_n(:) )

	sboltz_n = sum(boltz(1:nsol))

	PROB(1) = boltz(1) / sboltz_n

	do j = 2, nsol

		PROB(j) = PROB(j-1) + boltz(j) / sboltz_n

	end do

	sol = 0
	Random = ran2(Seed)
	j = 1

	do while ( sol == 0 )

		if( Random < PROB(j) ) sol = j

		j = j + 1

	end do

	sol_n = sol

	do i = 1, 3	  

		Xb_new(mapf(i)) = r_sol(i,1,sol)
		Yb_new(mapf(i)) = r_sol(i,2,sol) 
		Zb_new(mapf(i)) = r_sol(i,3,sol)

	end do

	! Set the r values for beads 1-3 to the new r values.

	r(1,:) = r_sol(1,:,sol)
	r(2,:) = r_sol(2,:,sol)
	r(3,:) = r_sol(3,:,sol)

	! Find the new u's

	do i = -1, 4

		u(i,:) = r(i,:) - r(i-1,:)

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

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

	end do

	! Form matrix for the new Jacobian.

	BB_n = 0.0

	BB_n = 0.0

	! BB_m(1:3,i) = ui x ( r4-ri )

	do i = 1, 3

		y(:) = r(4,:) - r(i,:)

		BB_n(1,i) = u(i,2) * y(3) - u(i,3) * y(2)  
		BB_n(2,i) = u(i,3) * y(1) - u(i,1) * y(3)  
		BB_n(3,i) = u(i,1) * y(2) - u(i,2) * y(1)

	end do  

	! Periodic boundary conditions.

	Xb_tmp = Xb_new
	Yb_tmp = Yb_new
	Zb_tmp = Zb_new

	do i = 1, Nbs
		
		if( Xb_tmp(i) > BoxSize )  Xb_tmp(i) = Xb_tmp(i) - &
									BoxSize * aint( Xb_tmp(i) / BoxSize )
		if( Yb_tmp(i) > BoxSize )  Yb_tmp(i) = Yb_tmp(i) - &
									BoxSize * aint( Yb_tmp(i) / BoxSize )
		if( Zb_tmp(i) > BoxSize )  Zb_tmp(i) = Zb_tmp(i) - &
									BoxSize * aint( Zb_tmp(i) / BoxSize )

		if( Xb_tmp(i) < 0.0 )  Xb_tmp(i) = Xb_tmp(i) - &
								BoxSize * aint( Xb_tmp(i) / BoxSize - 1 )
		if( Yb_tmp(i) < 0.0 )  Yb_tmp(i) = Yb_tmp(i) - &
								BoxSize * aint( Yb_tmp(i) / BoxSize - 1 )
		if( Zb_tmp(i) < 0.0 )  Zb_tmp(i) = Zb_tmp(i) - &
								BoxSize * aint( Zb_tmp(i) / BoxSize - 1 )

	end do



	! *********************************************
	! Now, the problem must be solved in reverse. *
	! *********************************************

	!phi0_m = phi_sol(0,sol) 

	phi1_m = phi_sol(1,sol) 

	th = th_m

	! T0lab does not need to be recalculated.

	! Set phi0 = phi0_m

	phi(0) = phi0_m

	! form T0, a 3x3 matrix.
	! Ti =	cos(th)				sin(th)				0
	!		sin(th)cos(ph)		-cos(th)cos(ph)		sin(ph)
	!		sin(th)sin(ph)		-cos(th)sin(ph)		-cos(ph)

	cth = dcos(th(0))
	sth = dsin(th(0))
	cphi = dcos(phi(0))
	sphi = dsin(phi(0))

	T0(1,1) = cth
	T0(1,2) = sth
	T0(1,3) = 0.0

	T0(2,1) = sth * cphi
	T0(2,2) = -cth * cphi
	T0(2,3) = sphi

	T0(3,1) = sth * sphi
	T0(3,2) = -cth * sphi
	T0(3,3) = -cphi

	! Find r1

	MM = matmul(T0lab, T0)

	y = 0.0
	y(1) = l(1)

	z(:) = r(0,:) + matmul(MM(:,:), y(:))

	r_sol(1,1,:) = z(1)
	r_sol(1,2,:) = z(2)
	r_sol(1,3,:) = z(3)

	! Find the new u1

	u(1,:) = z(:) - r(0,:)

	l(1) = sqrt(dot_product(u(1,:),u(1,:)))

	u(1,:) = u(1,:) / l(1)

	! form T1lab

	! x0 = u1
	! z0 = u1 x u0 / sin(th0) ! says sin(th1) in paper !!
	! y0 = z0 x x0 = z0 x u1

	!sth = sin(th(0))

	z(1) = ( u(1,2) * u(0,3) - u(1,3) * u(0,2) ) / sth 
	z(2) = ( u(1,3) * u(0,1) - u(1,1) * u(0,3) ) / sth 
	z(3) = ( u(1,1) * u(0,2) - u(1,2) * u(0,1) ) / sth 

	y(1) = ( z(2) * u(1,3) - z(3) * u(1,2) )  
	y(2) = ( z(3) * u(1,1) - z(1) * u(1,3) )  
	y(3) = ( z(1) * u(1,2) - z(2) * u(1,1) )  

	T1lab(:,1) = u(1,:)

	T1lab(1,2) = y(1)
	T1lab(2,2) = y(2)
	T1lab(3,2) = y(3)

	T1lab(1,3) = z(1)
	T1lab(2,3) = z(2)
	T1lab(3,3) = z(3)

	r41(:) = matmul( transpose(T1lab(:,:)), r(4,:)-r(0,:) )
	!u61(:) = matmul( transpose(T1lab(:,:)), u(6,:) )

	cth = dcos(th(2)) 
	sth = dsin(th(2)) 

	q1(1) = l(3) + l(2)*cth
	q1(2) = l(2)*sth
	q1(3) = 0.0

	qp1(1) = l(2) + l(3)*cth
	qp1(2) = l(3)*sth
	qp1(3) = 0.0

	!cth = dcos(th(4)) 
	!sth = dsin(th(4)) 

	q2(1) = l(4)
	q2(2) = 0.0
	q2(3) = 0.0

	!l1 = 0.0
	!l1(1) = l(1)

	! Search for solutions to F3.
	! Start by looping over the 2 possible solution branches.

	nsol = 0

	!open(14, file='f5.dat')

	do i = 1, 2

		F3a1_sol = .false.
		F3a2_sol = .false.

		phi(1) = phi1_m	- phi_span

		do while ( phi(1) < phi1_m + phi_span + phi_inc )

			nb = 0

			gate = 0

			F3a2_num = F3( i, l(1:4), th(1:3), phi(1:3), q1, q2, qp1, r41, &
							gate, T1, T2, F3a2_sol )

			if( F3a2_sol ) then

				gate = 0
						
				phi(1) = phi(1) + phi_inc

				F3b2_num = F3( i, l(1:4), th(1:3), phi(1:3), &
								q1, q2, qp1, r41, &
								gate, T1, T2, F3b2_sol )

				dF32 = ( F3b2_num - F3a2_num ) / phi_inc

			end if

			if( F3a1_sol .OR. F3a2_sol ) then

				both = .false.

				if( F3a1_sol .AND. F3a2_sol	) both = .true.

				if( both .AND. F3a1_num * F3a2_num < 0.0 ) then
					
					! Routine to make sure there is only one solution 
					! in the interval.
		
					x1 = phi(1) - phi_inc - phi_w
					x2 = phi(1) - phi_inc 

					nb = 10

					gate = 0

					a = phi(1)

					call zbrak(F3,x1,x2,nblocks,xb1,xb2,nb, &
								i, l(1:4), th(1:3), phi(1:3), &
								q1, q2, qp1, r41, gate)

					phi(1) = a
					
				else if( both .AND. dF31 * dF32 < 0 .AND. f3a1_num * dF31 < 0 ) then

					! Use a smaller range to test for solution.

					x1 = phi(1) - phi_inc - phi_w
					x2 = phi(1) - phi_inc 

					if( F3b2_sol ) x2 = x2 + phi_inc

					nb = 10

					gate = 0

					a = phi(1)

					call zbrak(F3,x1,x2,nblocks*2,xb1,xb2,nb, &
								i, l(1:4), th(1:3), phi(1:3), &
								q1, q2, qp1, r41, gate)

					phi(1) = a

				else if( .NOT. both .AND. F3a2_sol ) then

					x1 = phi(1) - phi_inc 
					x2 = phi(1)	- phi_w	- phi_inc

					if( x2 < phi1_m - phi_span ) then
					
						phi(1) = phi(1) + phi_w

						if( F3a2_sol ) phi(1) = phi(1) - phi_inc

						F3a1_num = F3a2_num
						F3a1_sol = F3a2_sol
						dF31 = dF32
					
						cycle

					end if

					a = phi(1)

					F3l_sol = .false.

					F3_lim = F3a1_num

					j = 0

					do while( .NOT. F3l_sol .AND. j < 20 )

						gate = -nint(F3_lim)

						if( j > 0 .AND. gate == 2 ) gate = 3

						b = zbrent4(F3,x1,x2,1.0e-9, &
									i, l(1:4), th(1:3), phi(1:3), &
									q1, q2, qp1, r41, gate, &
									T1, T2 )

						if( nint( b ) == 1000 ) exit

						if( gate == 3 ) gate = 2

						gate = -gate

						phi(1) = b

						F3_lim = F3( i, l(1:4), th(1:3), phi(1:3), &
									 q1, q2, qp1, r41, &
									 gate, T1, T2, F3l_sol )

						x2 = phi(1)

						j = j + 1

						if( j > 19 ) write(*,*) ' Caught in loop '

					end do

					if( nint( b ) == 1000 ) then
					
						phi(1) = a
						phi(1) = phi(1) + phi_w

						if( F3a2_sol ) phi(1) = phi(1) - phi_inc

						F3a1_num = F3a2_num
						F3a1_sol = F3a2_sol
						dF31 = dF32
					
						cycle

					end if
					
					x2 = phi(1)

					nb = 10
					
					call zbrak(F3,x1,x2,nblocks,xb1,xb2,nb, &
								i, l(1:4), th(1:3), phi(1:3), &
								q1, q2, qp1, r41, gate)

					phi(1) = a

				else if ( .NOT. both .AND. F3a1_sol ) then
				
					x1 = phi(1) - phi_w
					x2 = phi(1)

					if( x2 > phi1_m + phi_span + phi_inc ) then

						phi(1) = phi(1) + phi_w

						if( F3a2_sol ) phi(1) = phi(1) - phi_inc

						F3a1_num = F3a2_num
						F3a1_sol = F3a2_sol
						dF31 = dF32
					
						cycle

					end if

					a = phi(1)

					F3l_sol = .false.

					F3_lim = F3a2_num

					j = 0

					do while( .NOT. F3l_sol .AND. j < 20 )

						gate = -nint(F3_lim)

						if( j > 0 .AND. gate == 2 ) gate = 3

						b = zbrent4(F3,x1,x2,1.0e-9, &
									i, l(1:4), th(1:3), phi(1:3), &
									q1, q2, qp1, r41, &
									gate, T1, T2 )

						if( nint( b ) == 1000 ) exit

						if( gate == 3 ) gate = 2

						gate = -gate

						phi(1) = b

						F3_lim = F3( i, l(1:4), th(1:3), phi(1:3), &
									 q1, q2, qp1, r41, &
									 gate, T1, T2, F3l_sol )

						x2 = phi(1)

						j = j + 1

						if( j > 19 ) write(*,*) ' Caught in loop '

					end do

					if( nint( b ) == 1000 ) then
					
						phi(1) = a
						phi(1) = phi(1) + phi_w

						if( F3a2_sol ) phi(1) = phi(1) - phi_inc

						F3a1_num = F3a2_num
						F3a1_sol = F3a2_sol
						dF31 = dF32
					
						cycle

					end if

					x2 = phi(1)

					nb = 10

					call zbrak(F3,x1,x2,nblocks,xb1,xb2,nb, &
								i, l(1:4), th(1:3), phi(1:3), &
								q1, q2, qp1, r41, gate)

					phi(1) = a

				end if

				! Routine to find solution(s).

				gate_tmp = 0

				do j = 1, nb

					nsol = nsol + 1

					a = phi(1)

					if( j == nb ) gate_tmp = gate

					phi_sol(1,nsol) = zbrent4(F3,xb1(j),xb2(j),1.0e-7, &
												i, l(1:4), th(1:3), phi(1:3), &
												q1, q2, qp1, r41, &
												gate_tmp, T1, T2 )

					phi(1) = a

					if( nint( phi_sol(1,nsol) ) == 1000 ) then
					
						nsol = nsol - 1
						cycle
					
					end if
					
					phi_sol(0,nsol) = phi(0)
					phi_sol(2,nsol) = phi(2)