bigread.f90

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

subroutine BigRead( InputFile, InputConf, Nham, Nljgrs, Niongrs, Nsp, &
					MaxSp, MaxSteps, MaxBeads, MaxInt, MaxReal, MaxEE, &
					nmoves, Uwidth, BETA, ZETA, NAMElj, MASSlj, &
					EPS, SIG, CP, ALP, RMAX, KAPPA, LAMDA, GAMMA, &
					NAMEion, MASSion, CHARGE, NAMEsp, &
					NSTEPS, LENLJ, LENION, BEADTYPE, GROUPTYPE, &
					BEADDAMP, NTRIALS, STEPSTART, STEPLENGTH, &
					METHOD, INTPARAM , REALPARAM , BONDSAPART, &
					ERSTEPS, ERSTART, EREND, EESTEPS, EEW, &
					MASSsp, FromDisk, Nmol, StartConf, Volume, &
					PROB_MOVE, PROB_DR, PROB_SP_CD, PROB_SP_RG, &
					Nequil, Nprod, Ncollect, Neew, &
					NinCycle, Nstorage, Nadjust, histtype, &
					Nresmol, Nrefresh, Nrefrmoves, &
					DXYZ, DROT, tag_val, tag_mol, &
					Ubins, Umin, Umax, Nmin, Nmax, LNW, &
					Nstages, NST_TOT, NST_WT, STOP_H, Seed )

implicit none

! InputFile contains the parameters for the run.
! InputConf has the information for continuing a previous run.

character*30, intent(in)								:: InputFile
character*30, intent(in)								:: InputConf

! Nham is the number of hamiltonians.
! Nljgrs is the number of Lennard-Jones groups in the simulation.
! Niongrs is the number of ionic groups in the simulation.
! Nsp is the number of species in the simulation.

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

! MaxSp is the maximum number of species in the simulation.
! MaxSteps is the maximum number of CB steps for a molecule.
! MaxBeads is the maximum number of LJ and ionic beads in a molecule.
! MaxInt is the maximum number of integer parameters for a CB method.
! MaxReal is the maximum number of real parameters for a CB method.

integer, intent(in)										:: MaxSp
integer, intent(in)										:: MaxSteps
integer, intent(in)										:: MaxBeads
integer, intent(in)										:: MaxInt
integer, intent(in)										:: MaxReal
integer, intent(in)										:: MaxEE

integer, intent(in)										:: nmoves

! Uwidth is the energy width of the histogram bins.
! beta is the reciprical temperature.
! mu is the chemical potential term

real, intent(out)										:: Uwidth
real, dimension(Nham), intent(out)						:: BETA
real, dimension(MaxSp,Nham), intent(out)				:: ZETA

! NAMElj containes the names of the LJ groups.
! MASSlj contains the mass of the LJ groups.
! EPS contains the eps_ij parameters for each hamiltonian.
! SIG contains the sigma_ij parameters for each hamiltonian.
! CP contains the C_ij parameters for each hamiltonian.
! ALP contains the alpha_ij parameters for each hamiltonian.
! RMAX contains the Rmax parameters for each hamiltonian.
! KAPPA contains the correction to the epsilom_ij cross parameters.
! LAMDA contains the correction to the sigma_ij cross parameters.

character*15, dimension(Nljgrs), intent(out)			:: NAMElj
real, dimension(Nljgrs), intent(out)					:: MASSlj
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: EPS
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: SIG
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: CP
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: ALP
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: RMAX
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: KAPPA
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: LAMDA
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: GAMMA

! NAMEion containes the names of the ionic groups.
! MASSion contains the mass of the ionic groups.
! CHARGE contains the charge of a bead for each hamiltonian.

character*15, dimension(Niongrs), intent(out)			:: NAMEion
real, dimension(Niongrs), intent(out)					:: MASSion
real, dimension(Niongrs,Nham), intent(out)				:: CHARGE

! NAMEsp contains the name of each species.
! NSTEPS is the number of CB steps needed to grow each species.
! LENLJ is the LJ length of each species.
! LENION is the ionic length of each species.
! BEADTYPE indicates whether a bead is LJ or ionic.
! GROUPTYPE indicates the group identity of each bead.
! MASSsp contains the mass of the species.
! NTRIALS is the number of trials for each CB step.
! STEPSTART is the starting bead for each CB step.
! STEPLENGTH is the bead length of each CB step.
! METHOD is the method used for each CB step.
! INTPARAM contains the integer parameters for each CB step.
! REALPARAM contains the real parameters for each CB step.
! BONDAPART indicates how many bonds separate two beads.

character*15, dimension(Nsp), intent(out)				:: NAMEsp
integer, dimension(MaxSp), intent(out)					:: NSTEPS
integer, dimension(MaxSp), intent(out)					:: LENLJ 
integer, dimension(MaxSp), intent(out)					:: LENION
character*5, dimension(MaxBeads,MaxSp), intent(out)		:: BEADTYPE
integer, dimension(MaxBeads,MaxSp), intent(out)			:: GROUPTYPE 
real, dimension(MaxBeads,MaxEE,MaxSp), intent(out)		:: BEADDAMP
integer, dimension(MaxSteps,MaxSp), intent(out)			:: NTRIALS
integer, dimension(MaxSteps,MaxSp), intent(out)			:: STEPSTART 
integer, dimension(MaxSteps,MaxSp), intent(out)			:: STEPLENGTH
character*10, dimension(MaxBeads,MaxSp), intent(out)	:: METHOD
integer, dimension(MaxInt,MaxBeads,MaxSp), intent(out)	:: INTPARAM	
real, dimension(MaxReal,MaxBeads,MaxSp), intent(out)	:: REALPARAM
integer, dimension(MaxBeads,MaxBeads,MaxSp), intent(out):: BONDSAPART

integer, dimension(MaxSp), intent(out)					:: ERSTEPS
integer, dimension(MaxSteps,MaxSp), intent(out)			:: ERSTART
integer, dimension(MaxSteps,MaxSp), intent(out)			:: EREND
integer, dimension(MaxSp), intent(out)					:: EESTEPS
real, dimension(MaxEE,MaxSp), intent(out)				:: EEW

real, dimension(MaxSp), intent(out)						:: MASSsp

! FromDisk indicates whether or not to read the initial configuration from disk.

logical, intent(out)									:: FromDisk

! Nmol is the total and per species number of molecules in the simulation.
 
integer, dimension(0:MaxSp), intent(out)				:: Nmol

! StartConf is the starting configuration number.

integer, intent(out)									:: StartConf

! BoxSize contains the length of the simulation boxes.

real, intent(out)										:: Volume

! PROB_MOVE is the accumulative probability of performing a 
! Displacement / Rotation, Volume change, transfer, or regrowth.
! PROB_DR is the accumulative probability of performing a Displacement or Rotation.
! PROB_SP_CD is the accumulative probability of selecting a given species for creation/destruction.
! PROB_SP_RG is the accumulative probability of selecting a given species for regrowth.

real, dimension(nmoves), intent(out)					:: PROB_MOVE 
real, dimension(2), intent(out)							:: PROB_DR 
real, dimension(MaxSp), intent(out)						:: PROB_SP_CD
real, dimension(MaxSp), intent(out)						:: PROB_SP_RG

! Nequil is the number of equilibration steps.
! Nprod is the number of production steps.
! Nblocks is the number of MC steps before histogram data collection.
! NinCycle is the number of MC steps per cycle.
! Nstorage is the number of MC steps before storing the current configuration.
! Nadjust is the number of MC steps before adjusting the maximum displacement, etc.

integer, intent(out)									:: Nequil
integer, intent(out)									:: Nprod
integer, intent(out)									:: Ncollect
integer, intent(out)									:: Neew
integer, intent(out)									:: NinCycle
integer, intent(out)									:: Nstorage
integer, intent(out)									:: Nadjust
integer, intent(out)									:: Nresmol
integer, intent(out)									:: Nrefresh
integer, intent(out)									:: Nrefrmoves

character*10, intent(out)								:: histtype

! DXYZ is the maximum displacement per species.
! DROT is the maximum rotation per species.
! CDV is the constant volume change.

real, dimension(MaxSp), intent(out)					:: DXYZ
real, dimension(MaxSp), intent(out)					:: DROT

integer, dimension(MaxSp), intent(inout)			:: tag_val
integer, dimension(MaxSp), intent(inout)			:: tag_mol

! Ubins is the number of energy histogram bins.
! Umin is the minimum energy in the energy histogram.
! Umax is the maximum energy in the energy histogram.

integer, intent(inout)								:: Ubins
real, intent(out)									:: Umin
real, intent(out)									:: Umax

! Nmin is the minimum number of particles in the previous simulation.
! Nmax is the maximum number of particles in the previous simulation.

integer, dimension(MaxSp), intent(out)				:: Nmin
integer, dimension(MaxSp), intent(out)				:: Nmax

real, dimension(Nham), intent(out)					:: LNW

! Nstages is the number of stages for determining the weights
! NST_TOT is the number of MC steps in each stage.
! NST_WT is the number of MC steps before averaging in each stage.
! STOP_H is the decision hamiltonian number for which the weights will 
!	be changed after that stage.

integer													:: Nstages
integer, dimension(Nstages), intent(out)				:: NST_TOT
integer, dimension(Nstages), intent(out)				:: NST_WT
integer, dimension(Nstages), intent(out)				:: STOP_H

! Seed is the current random number generator seed value.

integer, intent(out)								:: Seed


! Local Stuff
integer										:: h, i, j, k, m
integer										:: bst, bend
integer										:: group
integer										:: step, repeat
integer										:: ns, nb, nh
integer										:: Nmol1
integer										:: ref1, ref2
integer										:: len
integer, dimension(MaxBeads)				:: REFBEAD
integer, dimension(MaxBeads)				:: BASE
integer, dimension(MaxBeads)				:: BONDS

character*5									:: btype
character*15								:: Species
character*30								:: InFile

real, parameter								:: Pi = 3.14159265359
real, parameter								:: Nav = 6.02214e23

real, external								:: zbrent2, funk2, zbrent3, funk3


open( 20, file = InputFile )

read(20,*)
read(20,*)
read(20,*)
read(20,*) BETA(1:Nham) 

BETA = 1.0 / BETA

read(20,*)
read(20,*) ZETA(1:Nsp, 1:Nham)
read(20,*)
read(20,*) Volume
read(20,*)
read(20,*) Uwidth
read(20,*)
read(20,*)
read(20,*)
read(20,*)
read(20,*)

do i = 1, Nljgrs

	read(20,*) j, NAMElj(j), MASSlj(j)

end do

read(20,*)

! Sigma is read in from input file.

do i = 1, Nljgrs * Nham

	read(20,*) j, k, EPS(k,k,j), SIG(k,k,j), CP(k,k,j), ALP(k,k,j)

end do

do h = 1, Nham

	read(20,*)
	read(20,*)
	read(20,*) ( ( KAPPA( i, j, h ), j = 1, Nljgrs ), i = 1, Nljgrs )
	read(20,*)
	read(20,*) ( ( LAMDA( i, j, h ), j = 1, Nljgrs ), i = 1, Nljgrs )
	read(20,*)
	read(20,*) ( ( GAMMA( i, j, h ), j = 1, Nljgrs ), i = 1, Nljgrs )

end do

! The following combining rules are used
! eps_ij = sqrt( eps_ii * eps_jj )
! sig_ij = 0.5*( sig_ii + sig_jj )
! alp_ij = sqrt( alp_ii * alp_jj )

do h = 1, Nham

	do i = 1, Nljgrs - 1

		do j = i + 1, Nljgrs

			EPS(i,j,h) = sqrt( EPS(i,i,h) * EPS(j,j,h) ) 
			EPS(j,i,h) = EPS(i,j,h)	* ( 1 - KAPPA(j,i,h) )
			EPS(i,j,h) = EPS(i,j,h)	* ( 1 - KAPPA(i,j,h) )

			ALP(i,j,h) = sqrt( ALP(i,i,h) * ALP(j,j,h) ) 
			ALP(j,i,h) = ALP(i,j,h)	* ( 1 - GAMMA(j,i,h) )
			ALP(i,j,h) = ALP(i,j,h)	* ( 1 - GAMMA(i,j,h) )

			SIG(i,j,h) = 0.5 * ( SIG(i,i,h) + SIG(j,j,h) ) 
			SIG(j,i,h) = SIG(i,j,h)	* ( 1 - LAMDA(j,i,h) )
			SIG(i,j,h) = SIG(i,j,h)	* ( 1 - LAMDA(i,j,h) )

			if( CP(i,i,h) > 0.0 ) then

				CP(i,j,h) = sqrt( CP(i,i,h) * CP(j,j,h) ) 
				CP(j,i,h) = CP(i,j,h)

				! Convert from sig_ij to rm_ij

				SIG(i,j,h) = SIG(i,j,h) / zbrent3( funk3, 0.5, 1.0, 1.0e-7, &
													ALP(i,j,h), CP(i,j,h) )

				SIG(j,i,h) = SIG(j,i,h) / zbrent3( funk3, 0.5, 1.0, 1.0e-7, &
													ALP(j,i,h), CP(j,i,h) )

				! find ratio of rmax_ij / rm_ij

				RMAX(i,j,h) = zbrent2( funk2, 0.0, 0.8, 1.0e-7, &
										ALP(i,j,h), CP(i,j,h) )

				RMAX(j,i,h) = zbrent2( funk2, 0.0, 0.8, 1.0e-7, &
										ALP(j,i,h), CP(j,i,h) )

			else

				CP(i,j,h) = -1.0 
				CP(j,i,h) = -1.0

			end if

		end do

	end do

end do

do i = 1, Nljgrs 

	do h = 1, Nham

		if( CP(i,i,h) > 0.0 ) then

			! Convert from sig_ii to rm_ii

			SIG(i,i,h) = SIG(i,i,h) / zbrent3( funk3, 0.5, 1.0, 1.0e-7, &
												ALP(i,i,h), CP(i,i,h) )

			! find ratio of rmax_ii / rm_ii

			RMAX(i,i,h) = zbrent2( funk2, 0.0, 0.8, 1.0e-7, &
									ALP(i,i,h), CP(i,i,h) )

			! Water simulations

!			if( i == 1 ) then

!				RMAX(i,i,h) = 0.45

!			end if

		end if

	end do

end do

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

if( Niongrs > 0 ) then

	read(20,*)

	do i = 1, Niongrs

		read(20,*) j, NAMEion(j), MASSion(j)

	end do

	read(20,*)

	do i = 1, Niongrs * Nham

		read(20,*) j, k, CHARGE(k,j)

	end do

end if

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

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

do i = 1, Nsp

	read(20,*)
	read(20,*) NAMEsp(i), ERSTEPS(i), EESTEPS(i), NSTEPS(i), LENLJ(i), LENION(i)
	read(20,*)

	bend = 0

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

		read(20,*) bst, bend, btype, group

		do j = bst, bend

			BEADTYPE(j,i) = btype
			GROUPTYPE(j,i) = group

			if( btype == 'LJ' ) then

				MASSsp(i) = MASSsp(i) + MASSlj(group)

			else if( btype == 'ION' ) then

				MASSsp(i) = MASSsp(i) + MASSion(group)

			end if

		end do

	end do

	read(20,*)

	step = 0

	do while ( step < NSTEPS(i) )

		read(20,*) step, NTRIALS(step,i), STEPSTART(step,i), bend

		STEPLENGTH(step,i) = bend - STEPSTART(step,i) + 1

		do k = STEPSTART(step,i), bend

			REFBEAD(k) = k

		end do

	end do

	read(20,*)

	do j = 1, ERSTEPS(i)

		read(20,*) ERSTART(j,i), ERSTART(j,i), EREND(j,i)

	end do

	! read in alpha values for expanded ensemble

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

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

		read(20,*) step, BEADDAMP(step,1:EESTEPS(i),i)

	end do