bigread.f90

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

	read(20,*)

	bend = 0

	do while ( bend < LENLJ(i) + LENION(i) )

		read(20,*) bst, METHOD(bst,i)

		bend = bst

		select case ( METHOD(bst,i) )

			case( 'Random' )

			case( 'Sphere' )

				read(20,*) INTPARAM(1,bst,i), REALPARAM(1,bst,i)

			case( 'Cone' )

				read(20,*) INTPARAM(1:3,bst,i), &
						   REALPARAM( 1:3*INTPARAM(3,bst,i), bst, i )

				REALPARAM( 1:2*INTPARAM(3,bst,i), bst, i ) = &
				REALPARAM( 1:2*INTPARAM(3,bst,i), bst, i ) * Pi / 180.0

			case( 'FxBend' ) 

				read(20,*) INTPARAM(1:2,bst,i)
				read(20,*) INTPARAM(3:2+INTPARAM(1,bst,i),bst,i), &
						   REALPARAM(1:1+2*INTPARAM(1,bst,i),bst,i)

				do j = 1, INTPARAM(1,bst,i)

					REALPARAM( 1+2*j, bst, i ) = REALPARAM( 1+2*j, bst, i ) * Pi / 180.0

				end do

			case( 'FxBendTor' )

				read(20,*) INTPARAM(1:4,bst,i)

				if( INTPARAM(3,bst,i) == 1 ) then

					read(20,*) INTPARAM(5:4+INTPARAM(1,bst,i)+2*INTPARAM(2,bst,i),bst,i), &
							   REALPARAM(1:1+2*INTPARAM(1,bst,i)+6*INTPARAM(2,bst,i),bst,i)

				else if( INTPARAM(3,bst,i) == 2 ) then

					read(20,*) INTPARAM(5:4+INTPARAM(1,bst,i)+2*INTPARAM(2,bst,i),bst,i), &
							   REALPARAM(1:1+2*INTPARAM(1,bst,i)+4*INTPARAM(2,bst,i),bst,i)

				end if
				
				do j = 1, INTPARAM(1,bst,i)

					REALPARAM( 1+2*j, bst, i ) = REALPARAM( 1+2*j, bst, i ) * Pi / 180.0

				end do


			case( 'Stretch' )

				read(20,*) INTPARAM(1,bst,i), REALPARAM(1:2,bst,i)

			case( 'StBend' ) 

				read(20,*) INTPARAM(1:2,bst,i)
				read(20,*) INTPARAM(3:2+INTPARAM(1,bst,i),bst,i), &
						   REALPARAM(1:2+2*INTPARAM(1,bst,i),bst,i)

				do j = 1, INTPARAM(1,bst,i)

					REALPARAM( 2+2*j, bst, i ) = REALPARAM( 2+2*j, bst, i ) * Pi / 180.0

				end do

			case( 'StBendTor' )

				read(20,*) INTPARAM(1:4,bst,i)

				if( INTPARAM(3,bst,i) == 1 ) then

					read(20,*) INTPARAM(5:4+INTPARAM(1,bst,i)+2*INTPARAM(2,bst,i),bst,i), &
							   REALPARAM(1:2+2*INTPARAM(1,bst,i)+6*INTPARAM(2,bst,i),bst,i)

				else if( INTPARAM(3,bst,i) == 2 ) then

					read(20,*) INTPARAM(5:4+INTPARAM(1,bst,i)+2*INTPARAM(2,bst,i),bst,i), &
							   REALPARAM(1:2+2*INTPARAM(1,bst,i)+4*INTPARAM(2,bst,i),bst,i)

				end if
				
				do j = 1, INTPARAM(1,bst,i)

					REALPARAM( 2+2*j, bst, i ) = REALPARAM( 2+2*j, bst, i ) * Pi / 180.0

				end do

			case( 'ResRand' )

				read(20,*) INTPARAM(1:1,bst,i)

			case( 'ResRot' )

			case( 'ResCir' )

			case( 'ResCR4' )

			case( 'ResMatch' )

			case( 'Match' )

				read(20,*) INTPARAM(1,bst,i)

			case( 'Complete' )

		end select

	end do

end do

read(20,*)
read(20,*) FromDisk
read(20,*)

if( FromDisk ) open( 30, file = InputConf )

do i = 1, Nsp

	if( FromDisk ) then

		read(30,*) Species, Nmol1
		read(20,*)

	else 

		read(20,*) Species, Nmol1

	end if
	
	do j = 1, Nsp

		if( Species == NAMEsp(j) ) then

			Nmol(j) = Nmol1

		end if

	end do

end do

Nmol(0) = sum( Nmol(1:Nsp) )

if( FromDisk ) then

	read(30,*)
	read(30,*) StartConf, Seed
	read(30,*)
	read(30,*) tag_val(1:Nsp), tag_mol(1:Nsp)

	read(20,*)
	read(20,*)
	read(20,*)

else

	StartConf = 0
	
	read(20,*)
	read(20,*) Seed
	read(20,*)

	tag_val = 0
	tag_mol	= 0
 
end if

read(20,*) PROB_MOVE(1:nmoves)

do i = 2, nmoves

	PROB_MOVE(i) = PROB_MOVE(i) + PROB_MOVE(i-1)

end do

PROB_MOVE = PROB_MOVE / PROB_MOVE(nmoves)

read(20,*)
read(20,*) PROB_DR(1:2)

PROB_DR(2) = PROB_DR(1) + PROB_DR(2)

PROB_DR = PROB_DR / PROB_DR(2)

read(20,*)
read(20,*) PROB_SP_CD(1:Nsp)

do i = 2, Nsp

	PROB_SP_CD(i) = PROB_SP_CD(i) + PROB_SP_CD(i-1)

end do

PROB_SP_CD = PROB_SP_CD / PROB_SP_CD(Nsp)

read(20,*)
read(20,*) PROB_SP_RG(1:Nsp)

do i = 2, Nsp

	PROB_SP_RG(i) = PROB_SP_RG(i) + PROB_SP_RG(i-1)

end do

PROB_SP_RG = PROB_SP_RG / PROB_SP_RG(Nsp)

read(20,*)
read(20,*) Nequil, Nprod
read(20,*)
read(20,*) NinCycle
read(20,*)
read(20,*) Nstorage
read(20,*)
read(20,*) Nadjust
read(20,*)
read(20,*) Neew
read(20,*) 
read(20,*) Nrefresh
read(20,*) 
read(20,*) Nrefrmoves
read(20,*) 
read(20,*) Ncollect
read(20,*) 
read(20,*) histtype
read(20,*)
read(20,*) 

i = 1

do while( i <= Nstages )

	read(20,*) ns, nb, nh, repeat

	do j = i, i + repeat

		NST_TOT(j) = ns
		NST_WT(j) = nb
		STOP_H(j) = nh

	end do								  

	i = i + 1 + repeat

end do

if( FromDisk ) then

	read(30,*)
	read(30,*)
	read(30,*) DXYZ(1:Nsp)
	read(30,*) DROT(1:Nsp)
	read(30,*)
	read(30,*)
	read(30,*) Ubins, Umin, Umax
	read(30,*)
	read(30,*) Nmin(1), Nmax(1)

else

	DXYZ = 0.1
	DROT = 0.2

end if

if( FromDisk ) then

	read(30,*)
	read(30,*) 

	do i = 1, Nham

		read(30,*) LNW(i)

	end do

	LNW = log(LNW)

	do j = 1, Nsp

		read(30,*)
		read(30,*) 

		do i = 1, EESTEPS(j)

			read(30,*) EEW(i,j)

		end do

	end do

else

!	LNW = -log( real( Nham ) )

	LNW = -1.0e50
	LNW(1) = 0.0

	EEW = 0.0

end if

close(20)
close(30)


BONDSAPART = 0

do i = 1, Nsp

	InFile = 'ba_'//char(48+i)//'.dat'

	open( 33, file = InFile )

	k = LENLJ(i) + LENION(i)

	do j = 1, k

		read(33,*) BONDSAPART(j,1:k,i)

	end do

	close(33)

end do

return

! Code below is no longer used.

do i = 1, Nsp

	BASE(1) = 1
	BONDS(1) = 0

	do j = 1, LENLJ(i) + LENION(i)

		select case ( METHOD(j,i) )

			case( 'Sphere' )

				BASE(j) = INTPARAM(1,j,i)
				BONDS(j) = 1

			case( 'Cone' )

				do k = j, j + INTPARAM(3,j,i) - 1

					BASE(k) = INTPARAM(2,j,i)
					BONDS(k) = 1
				
				end do

			case( 'ConeTor' )

				BASE(j) = INTPARAM(3,j,i)
				BONDS(j) = 1

			case( 'FxBend' )

				BASE(j) = INTPARAM(2,j,i)
				BONDS(j) = 1

			case( 'FxBendTor' )

				BASE(j) = INTPARAM(4,j,i)
				BONDS(j) = 1

			case( 'Stretch' )

				BASE(j) = INTPARAM(1,j,i)
				BONDS(j) = 1

			case( 'StBend' )

				BASE(j) = INTPARAM(2,j,i)
				BONDS(j) = 1

			case( 'StBendTor' )

				BASE(j) = INTPARAM(4,j,i)
				BONDS(j) = 1

			case( 'Match' )

				BASE(j) = INTPARAM(1,j,i)
				BONDS(j) = 0

			case( 'Complete' )

		end select

	end do

	do j = 1, LENLJ(i) + LENION(i) - 1
	
		do k = j + 1, LENLJ(i) + LENION(i)
		
			ref1 = BASE(j)
			ref2 = BASE(k)

			len = BONDS(k)

			m = 0
			
			do while ( m == 0 )

				if( ref1 == ref2 ) then

					len = len + BONDS(j)

					m = 1
			
				else if( ref1 > ref2 ) then

					if(	ref2 /= j .AND. ref1 /= k ) then

						len = len + BONDS(ref1)
						ref1 = BASE(ref1)
						
					else

						m = 1	

					end if
				
				else if( ref2 > ref1 ) then

					if(	ref2 /= j .AND. ref1 /= k ) then

						len = len + BONDS(ref2)
						ref2 = BASE(ref2)

					else

						m = 1

					end if

				end if

			end do
				
			BONDSAPART(j,k,i) = len
			BONDSAPART(k,j,i) = len
			
		end do
		
	end do
	
end do						 

return

end subroutine BigRead