grow.f90

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

Subroutine Grow( NSteps, STEPSTART, STEPLENGTH, NTRIALS, lenlj, lenion, &  
				 MaxBeads, METHOD, MaxInt, MaxReal, INTPARAM, REALPARAM, &
				 BEADTYPE, New, StartGrowthStep, &
				 Xlj_new, Ylj_new, Zlj_new, &
				 TYPElj_new, DAMPlj2_new, DAMPlj3_new, &
				 Xion_new, Yion_new, Zion_new, &
				 TYPEion_new, DAMPion_new, &
				 Nham, ULJBOX, ULJLR, ULJ_MOL, UREAL, USURF, &
				 UFOURIER, USELF_CH, USELF_MOL, UION_MOL, Uintra, &
				 Niongrs, CHARGE, Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
				 SUMQX, SUMQY, SUMQZ, EXPX_new, EXPY_new, EXPZ_new, &
				 SUMQEXPV, BETA, BIGW, LNW, &
				 Nlj, Xlj, Ylj, Zlj, &
				 TYPElj, DAMPlj2, DAMPlj3, &
				 Nion, Xion, Yion, Zion, &
				 TYPEion, DAMPion, &
				 BoxSize, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
				 XLRCORR, ELRCORR, &
				 Nres, Nresmol, reslen, Xresmols, &
				 Yresmols, Zresmols, Uint_resm, &
				 Ulj_resm, Uion_resm, &
				 BONDSAPART, f14_lj, f14_ion, Seed )
				 
implicit none

! This subroutine grows a molecule using configurational bias. 

! NSteps is the number of configurational bias steps it takes to grow the molecule.

integer, intent(in)										:: NSteps

! STEPSTART is the start bead of each CB step.
! STEPLENGTH is the number of beads in each CB step.
! NTRIALS is the number of trial orientations/locations for each CB step.

integer, dimension(NSteps)								:: STEPSTART
integer, dimension(NSteps)								:: STEPLENGTH
integer, dimension(NSteps)								:: NTRIALS

! lenlj is the number of LJ beads in the molecule to be grown.
! lenion is the number of ionic beads in the molecule to be grown.

integer, intent(in)										:: lenlj, lenion

! METHOD is the method used to grow each bead.
! MaxInt is the maximum number of integer parameters.
! MaxReal is the maximum number of real	parameters.
! INTPARAM are the integer parameters needed by the method to grow each bead.
! REALPARAM are the real parameters needed by the method to grow each bead.
! BEADTYPE indicates whether the bead is 'LJ' or 'ION'.

integer, intent(in)										:: MaxBeads
character*10, dimension( MaxBeads ), intent(in)			:: METHOD
integer, intent(in)										:: MaxInt
integer, intent(in)										:: MaxReal
integer, dimension( MaxInt, MaxBeads ), intent(in)		:: INTPARAM
real, dimension( MaxReal, MaxBeads ), intent(in)		:: REALPARAM
character*5, dimension( MaxBeads ), intent(in)			:: BEADTYPE

! New is a logical to indicate if a new molecule is being grown or the
! Rosenbluth weight of an old molecule is being determined.
! New = .True. when a new molecule is to be grown.

logical, intent(in)										:: New

! StartGrowthStep is the starting step number for gegrowing the molecule.
! StartGrowthStep = 1, for regrowing the whole molecule.

integer, intent(in)										:: StartGrowthStep

! Xlj_new, Ylj_new, and Zlj_new contain the new coordinates of the grown
! molecule if New = .True., or the coordinates of an existing molecule if
! New = .False.
! When regrowing part of a molecule Xlj_new, etc. contain the coordinates 
! of the molecule up to StartGrowthStep - 1

real, dimension(lenlj), intent(inout)					:: Xlj_new
real, dimension(lenlj), intent(inout)					:: Ylj_new
real, dimension(lenlj), intent(inout)					:: Zlj_new

! TYPElj_new contains the group identity of the LJ beads to be grown.

integer, dimension(lenlj), intent(in)					:: TYPElj_new

real, dimension(lenlj), intent(inout)					:: DAMPlj2_new
real, dimension(lenlj), intent(inout)					:: DAMPlj3_new

! Xion_new, Yion_new, and Zion_new have a similar definition to that above for
! Xlj_new, Ylj_new, and Zlj_new.

real, dimension(lenion), intent(inout)					:: Xion_new
real, dimension(lenion), intent(inout)					:: Yion_new
real, dimension(lenion), intent(inout)					:: Zion_new

! TYPEion_new contains the group identity of the ionic beads to be grown.

integer, dimension(lenion), intent(in)					:: TYPEion_new

real, dimension(lenion), intent(inout)					:: DAMPion_new

! Nham is the number of hamiltonians.

integer, intent(in)										:: Nham

! ULJBOX is the LJ energy of the phase before and after the molecule is grown.
! ULJLR is the long range correction to the LJ energy.
! ULJ_MOL is the non bonded LJ energy of the molecule to be grown.

real, dimension(Nham), intent(inout)					:: ULJBOX
real, dimension(Nham), intent(inout)					:: ULJLR
real, dimension(Nham), intent(inout)					:: ULJ_MOL

! UREAL is the real coulombic energy.
! USURF is the surface energy.
! UFOURIER is the reciprical space coulombic energy.
! USELF_CH is the self charge energy.
! USELF_MOL is the self energy of the molecule to be grown.

real, dimension(Nham), intent(inout)					:: UREAL
real, dimension(Nham), intent(inout)					:: USURF
real, dimension(Nham), intent(inout)					:: UFOURIER	
real, dimension(Nham), intent(inout)					:: USELF_CH
real, dimension(Nham), intent(inout)					:: USELF_MOL
real, dimension(Nham), intent(inout)					:: UION_MOL

! Uintra is the intramolecular energy of the grown molecule.

real, intent(inout)										:: Uintra

! Niongrs is the number of ion groups.
! CHARGE contains the charge for a given group and hamiltonian.

integer, intent(in)										:: Niongrs
real, dimension(Niongrs, Nham), intent(in)				:: CHARGE

! Alpha is an Ewald sum parameter, Alpha = kappa * L, for kappa in A + T.

real, intent(in)										:: Alpha

! Kmax is an Ewald sum parameter.
! Nkvec is the number of k-vectors used in the Fourier sum.
! KX, KY, KZ contain the vector identity of the Nkvec vectors.
! CONST contains the constant part of the Fourier summation for a given Nkvec.

integer, intent(in)										:: Kmax
integer, intent(in)										:: Nkvec
integer, dimension(Nkvec), intent(in)					:: KX, KY, KZ
real, dimension(Nkvec), intent(in)						:: CONST

! SUMQX is the summation of qi * xi for all ions in the box before and after 
! the growth.

real, dimension(Nham), intent(inout)					:: SUMQX, SUMQY, SUMQZ

! EXPX_new, etc. contain exp( i*kx*x ) for the beads of the grown molecule.

complex, dimension(0:Kmax, lenion), intent(inout)		:: EXPX_new
complex, dimension(-Kmax:Kmax, lenion), intent(inout)	:: EXPY_new, EXPZ_new

! SUMQEXPV contains the summation of qi*exp(i*(kx*x + ky*y + kz*z)) 
! for a given k-vector and hamiltonian before and after the growth.

complex, dimension(Nkvec, Nham), intent(inout)			:: SUMQEXPV

! beta contains the recprical temperatures.

real, dimension(Nham), intent(in)						:: BETA

! BIGW is the Rosenbluth factor of the grown molecule.

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

! LNW contains the weight of each state.

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

! Nlj is the number of LJ beads in the phase the molecule is grown in.
! Xlj, Ylj, and Zlj are the coordinates of the Nlj LJ beads.
! TYPElj contains the group identity of the Nlj LJ beads.

integer, intent(in)										:: Nlj
real, dimension(Nlj), intent(in)						:: Xlj, Ylj, Zlj
integer, dimension(Nlj), intent(in)						:: TYPElj
real, dimension(Nlj), intent(in)						:: DAMPlj2, DAMPlj3

! Nion is the number of ionic beads in the phase the molecule is grown in.
! Xion, Yion, and Zion are the coordinates of the Nion ionic beads.
! TYPEion contains the group identity of the Nion ionic beads.

integer, intent(in)										:: Nion
real, dimension(Nion), intent(in)						:: Xion, Yion, Zion
integer, dimension(Nion), intent(in)					:: TYPEion
real, dimension(Nion), intent(in)						:: DAMPion

! BoxSize is the length of the simulation box.

real, intent(in)										:: BoxSize

! Nljgrs is the number of LJ groups in the system.
! EPS is a rank 3 array containing the eps_ij parameters for each hamiltonian.
! SIG is a rank 3 array containing the sigma_ij parameters for each hamiltonian.
! CP is a rank 3 array containing the C_ij parameters for each hamiltonian.
! ALP is a rank 3 array containing the alpha_ij parameters for each hamiltonian.
! RMAX is a rank 3 array containing the Rmax_ij parameters for each hamiltonian.
									
integer, intent(in)										:: Nljgrs
real, dimension(Nljgrs, Nljgrs, Nham), intent(in)		:: EPS, SIG, CP, ALP, RMAX

! XLRCORR and ELRCORR contain correction factors for the long range LJ energy.

real, dimension(605), intent(in)						:: XLRCORR  
real, dimension(605), intent(in)						:: ELRCORR

integer, intent(in)										:: Nres, Nresmol

integer, dimension(Nres), intent(in)					:: reslen

real, dimension(Nresmol, MaxBeads, Nres), intent(in)	:: Xresmols, Yresmols, Zresmols
real, dimension(Nresmol, Nres), intent(in)				:: Uint_resm
real, dimension(Nresmol, Nham, Nres), intent(in)		:: Ulj_resm, Uion_resm

! BONDSAPART gives the bondlengths any two LJ beads are away from each other.

integer, dimension(MaxBeads,MaxBeads)					:: BONDSAPART

! f14_lj and f14_ion are damping factors for the 1-4 intramolecular interactions

real, intent(in)										:: f14_lj
real, intent(in)										:: f14_ion

! Seed is the current random number generator seed value.

integer, intent(inout)									:: Seed

real, external											:: ran2

! Local Variables

integer													:: h, i, j, k, m, n
integer													:: ljb, ionb
integer													:: ljb_temp, ionb_temp
integer													:: lenljst, lenionst
integer													:: Ntry, Nb
integer													:: Selection
integer													:: ncount
integer													:: FxSt, newold
integer													:: nDoOver = 1000
integer													:: resl, resstart
integer													:: resmol
integer, dimension(Nljgrs)								:: NLJGROUPS
integer, dimension(lenlj+lenion)						:: LJBEAD, IONBEAD

real													:: CoulCombo
real													:: Random
real													:: xtr, ytr, ztr
real, dimension(Nham)									:: Utemp
real, dimension(Nham)									:: dULJ_MOL_part
real, dimension(Nham)									:: dUION_MOL_part
real, dimension(Nham)									:: dU
real, dimension(Nham)									:: W
real, dimension(Nham,NSteps)							:: SMALLW
real, dimension(Nljgrs, Nljgrs, Nham)					:: EPS_CP
real, dimension(MaxBeads)								:: Xres, Yres, Zres

real, allocatable, dimension(:,:)						:: BOLTZ
real, allocatable, dimension(:)							:: PROB
real, allocatable, dimension(:)							:: TempX, TempY, TempZ
real, allocatable, dimension(:,:)						:: Xlj_tr, Ylj_tr, Zlj_tr
real, allocatable, dimension(:,:)						:: Xion_tr, Yion_tr, Zion_tr
real, allocatable, dimension(:,:)						:: Xres_tr, Yres_tr, Zres_tr
real, allocatable, dimension(:,:)						:: SUMQX_tr, SUMQY_tr, SUMQZ_tr
real, allocatable, dimension(:,:)						:: dULJBOX_tr
real, allocatable, dimension(:,:)						:: dULJLR_tr
real, allocatable, dimension(:,:)						:: dULJ_MOL_tr
real, allocatable, dimension(:,:)						:: dUREAL_tr
real, allocatable, dimension(:,:)						:: dUSURF_tr
real, allocatable, dimension(:,:)						:: dUFOURIER_tr
real, allocatable, dimension(:,:)						:: dUSELF_CH_tr
real, allocatable, dimension(:,:)						:: dUSELF_MOL_tr
real, allocatable, dimension(:,:)						:: dUION_MOL_tr
real, allocatable, dimension(:,:)						:: UVIB
real, allocatable, dimension(:,:)						:: UBEND
real, allocatable, dimension(:,:)						:: UTOR

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

complex, allocatable, dimension(:,:,:)					:: EXPX_tr
complex, allocatable, dimension(:,:,:)					:: EXPY_tr
complex, allocatable, dimension(:,:,:)					:: EXPZ_tr
complex, allocatable, dimension(:,:,:)					:: SUMQEXPV_tr
complex, allocatable, dimension(:,:)					:: CZERO1, CZERO2


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

W = exp( LNW )

LJBEAD = 0
IONBEAD = 0

ljb = 0
ionb = 0

do i = 1, lenion + lenlj
	
	if(	BEADTYPE(i) == 'LJ' ) then

		ljb = ljb + 1
		
		LJBEAD(i) = ljb

		IONBEAD(i) = ionb
		
	else if( BEADTYPE(i) == 'ION' ) then

		ionb = ionb + 1

		IONBEAD(i) = ionb

		LJBEAD(i) = ljb

	end if

end do

ljb = 0
ionb = 0

! NLJGROUPS, for calculation of ULJLR.

NLJGROUPS = 0

! NLJGROUPS for complete molecules. 
	
do i = 1, Nlj

	NLJGROUPS( TYPElj(i) ) = NLJGROUPS( TYPElj(i) ) + 1

end do

! NLJGROUPS for molecule being grown. 

do i = 1, StartGrowthStep - 1

	do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1
	
		ljb = LJBEAD(k)
		ionb = IONBEAD(k)

		if( BEADTYPE(k) == 'LJ' ) then
		
			NLJGROUPS( TYPElj_new(ljb) ) = NLJGROUPS( TYPElj_new(ljb) ) + 1

		end if

	end do

end do

resstart = 0

BIGW = 1.0

do i = StartGrowthStep, NSteps

	lenljst = 0
	lenionst = 0

	do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1
	
		if(	BEADTYPE(k) == 'LJ' ) then

			lenljst = lenljst + 1
		
		else if( BEADTYPE(k) == 'ION' ) then

			lenionst = lenionst + 1

		end if

	end do

	Ntry = NTRIALS(i)

	if( lenljst > 0 ) then

		allocate( Xlj_tr( Ntry, ljb + 1 : ljb + lenljst ) )
		allocate( Ylj_tr( Ntry, ljb + 1 : ljb + lenljst ) )
		allocate( Zlj_tr( Ntry, ljb + 1 : ljb + lenljst ) )

	end if

	if( lenionst > 0 ) then

		allocate( Xion_tr( Ntry, ionb + 1 : ionb + lenionst ) )
		allocate( Yion_tr( Ntry, ionb + 1 : ionb + lenionst ) )
		allocate( Zion_tr( Ntry, ionb + 1 : ionb + lenionst ) )

		allocate( EXPX_tr( Ntry, 0:Kmax, ionb + 1 : ionb + lenionst ) )
		allocate( EXPY_tr( Ntry, -Kmax:Kmax, ionb + 1 : ionb + lenionst ) )
		allocate( EXPZ_tr( Ntry, -Kmax:Kmax, ionb + 1 : ionb + lenionst ) )

		allocate( SUMQEXPV_tr( Ntry, Nkvec, Nham ) )

		allocate( SUMQX_tr( Ntry, Nham ) )
		allocate( SUMQY_tr( Ntry, Nham ) )
		allocate( SUMQZ_tr( Ntry, Nham ) )

	end if
	
	allocate( Xres_tr( Ntry, MaxBeads ) )
	allocate( Yres_tr( Ntry, MaxBeads ) )
	allocate( Zres_tr( Ntry, MaxBeads ) )

	allocate( dULJBOX_tr( Ntry, Nham ) )
	allocate( dULJLR_tr( Ntry, Nham ) )
	allocate( dULJ_MOL_tr( Ntry, Nham ) )
	allocate( dUREAL_tr( Ntry, Nham ) )
	allocate( dUSURF_tr( Ntry, Nham ) )
	allocate( dUFOURIER_tr( Ntry, Nham ) )
	allocate( dUSELF_CH_tr( Ntry, Nham ) )
	allocate( dUSELF_MOL_tr( Ntry, Nham ) )
	allocate( dUION_MOL_tr( Ntry, Nham ) )

	allocate( UVIB( Ntry, lenlj + lenion ) )
	allocate( UBEND( Ntry, lenlj + lenion ) )
	allocate( UTOR( Ntry, lenlj + lenion ) )

	dULJBOX_tr = 0.0
	dULJLR_tr = 0.0
	dULJ_MOL_tr = 0.0
	dUREAL_tr = 0.0
	dUSURF_tr = 0.0
	dUFOURIER_tr = 0.0
	dUSELF_CH_tr = 0.0
	dUSELF_MOL_tr = 0.0
	dUION_MOL_tr = 0.0

	UVIB = 0.0
	UBEND = 0.0
	UTOR = 0.0

	allocate( BOLTZ( Nham, Ntry ) )
	allocate( PROB( Ntry ) )

	SMALLW(:,i) = 0.0

	do j = 1, Ntry

		! Find the coordinates of the trial positions.

		if( .NOT. New .AND. j == 1 ) then

			do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1
			
				ljb = LJBEAD(k)
				ionb = IONBEAD(k)

				if( BEADTYPE(k) == 'LJ' ) then
				
					NLJGROUPS( TYPElj_new(ljb) ) = NLJGROUPS( TYPElj_new(ljb) ) + 1

				end if

				if( BEADTYPE(k) == 'LJ' ) then
			
					Xlj_tr(j,ljb) = Xlj_new(ljb)
					Ylj_tr(j,ljb) = Ylj_new(ljb)
					Zlj_tr(j,ljb) = Zlj_new(ljb)
						
				else if( BEADTYPE(k) == 'ION' ) then
					
					Xion_tr(j,ionb) = Xion_new(ionb)
					Yion_tr(j,ionb) = Yion_new(ionb)
					Zion_tr(j,ionb) = Zion_new(ionb)

					EXPX_tr(j,:,ionb) = EXPX_new(:,ionb)
					EXPY_tr(j,:,ionb) = EXPY_new(:,ionb)
					EXPZ_tr(j,:,ionb) = EXPZ_new(:,ionb)

				end if

				if( resstart > 0 ) then
				
					do m = 1, resl

						Xres_tr(j,m) = Xres(m)
						Yres_tr(j,m) = Yres(m)
						Zres_tr(j,m) = Zres(m)

					end do												

				end if

				select case ( METHOD(k) )

					case( 'Stretch' )

						FxSt = 2
						newold = 2

						if( BEADTYPE(k) == 'LJ' ) then
			
							xtr = Xlj_new(ljb)
							ytr = Ylj_new(ljb)
							ztr = Zlj_new(ljb)
						
						else if( BEADTYPE(k) == 'ION' ) then
					
							xtr = Xion_new(ionb)
							ytr = Yion_new(ionb)
							ztr = Zion_new(ionb)

						end if

						call Stretch( lenlj, lenion, &
										Xlj_new, Ylj_new, Zlj_new, &
										Xion_new, Yion_new, Zion_new, &
										LJBEAD, IONBEAD, &
										BEADTYPE, MaxInt, MaxReal, &
										INTPARAM(:,k), REALPARAM(:,k), &
										BoxSize, Nham, BETA, &
										FxSt, newold, &
										xtr, ytr, ztr, &
										UVIB(j,k), Seed )  
					
					case( 'FxBend' )

						FxSt = 1
						newold = 2

						if( BEADTYPE(k) == 'LJ' ) then
			
							xtr = Xlj_new(ljb)
							ytr = Ylj_new(ljb)
							ztr = Zlj_new(ljb)
						
						else if( BEADTYPE(k) == 'ION' ) then