create.f90

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

subroutine Create( InputConf, Nham, Nljgrs, Niongrs, Nsp, &
				   MaxSp, MaxSteps, MaxBeads, &
				   MaxInt, MaxReal, MaxEE, &
				   BETA, ZETA, EPS, SIG, CP, ALP, RMAX, CHARGE, &
				   NSTEPS, LENLJ, LENION, BEADTYPE, GROUPTYPE, &
				   BEADDAMP, NTRIALS, STEPSTART, STEPLENGTH, &
				   METHOD, INTPARAM , REALPARAM , BONDSAPART, &
				   EESTEPS, FromDisk, BoxSize,  &
				   Nlj, Nion, MaxMol, MaxLJ, MaxIon, Nmol, &
				   Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
				   Xion, Yion, Zion, TYPEion, DAMPion, &
				   SPECIES, LENGTHlj, LENGTHion, &
				   STARTlj, STARTion,	&
				   LNW, LNPSI, LnPi, XLRCORR, ELRCORR, &
				   Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
				   EXPX, EXPY, EXPZ, SUMQEXPV, &
				   SUMQX, SUMQY, SUMQZ, &
				   ULJBOX, ULJLR, UREAL, UFOURIER, &
				   USURF, USELF_CH, USELF_MOL, &
				   ULJ_MOL, UION_MOL, UINTRA, &
				   tag_val, tag_mol, &
				   Nres, Nresmol, reslen, Xresmols, &
				   Yresmols, Zresmols, Uint_resm, &
				   Ulj_resm, Uion_resm, &
				   f14_lj, f14_ion, Seed )				  
				  
implicit none

! InputConf has the information for continuing a previous run.

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

! beta is the reciprocal temperature.

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

! zeta is a gcmc parameter = f(mu, T)
! zeta = q(T) * exp(beta * mu_B) from Frenkel + Smit

real, dimension(MaxSp,Nham), intent(in)					:: ZETA

! 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_ij parameters for each hamiltonian.

real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: EPS
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: SIG
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: CP
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: ALP
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: RMAX

! CHARGE contains the charge of a bead for each hamiltonian.

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

! 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.
! 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.
! BONDSAPART indicates how many bonds separate two beads.

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

integer, dimension(MaxSp), intent(in)					:: EESTEPS

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

logical, intent(in)										:: FromDisk

! BoxSize contains the length of the simulation boxes.

real, intent(in)										:: BoxSize

! Nlj is the number of LJ beads in the simulation boxes.
! Nion is the number of ionic beads in the simulation boxes.

integer, intent(out)									:: Nlj
integer, intent(out)									:: Nion

! MaxMol is the maximum number of molecules per box.
! MaxLJ is the maximum number of LJ beads per box.
! MaxIon is the maximum number of ionic beads per box.

integer, intent(in)										:: MaxMol
integer, intent(in)										:: MaxLJ
integer, intent(in)										:: MaxIon

! Nmol is the number of molecules in the simulation box.

integer, dimension(0:MaxSp), intent(in)					:: Nmol

! Xlj, Ylj, and Zlj are the coordinates of the LJ beads.
! Xion, Yion, and Zion are the coordinates of the ionic beads.

real, dimension( MaxLJ ), intent(out)					:: Xlj, Ylj, Zlj
real, dimension( MaxIon ), intent(out)					:: Xion, Yion, Zion

! TYPElj contains the group identity of each LJ bead.
! TYPEion contains the group identity of each ionic bead.

integer, dimension( MaxLJ ), intent(out)				:: TYPElj
integer, dimension( MaxIon ), intent(out)				:: TYPEion

real, dimension( MaxLJ ), intent(out)					:: DAMPlj2, DAMPlj3
real, dimension( MaxIon ), intent(out)					:: DAMPion

! SPECIES contains the species identity of each molecule.
! LENGTHlj contains the number of LJ beads in each molecule.
! LENGTHion	contains the number of ionic beads in each molecule.
! STARTlj contains the starting LJ bead number for each molecule.
! STARTion contains the starting ionic bead number for each molecule.

integer, dimension(MaxMol), intent(out)					:: SPECIES
integer, dimension(MaxMol), intent(out)					:: LENGTHlj, LENGTHion
integer, dimension(MaxMol), intent(out)					:: STARTlj, STARTion

! LNW contains the weight of each hamiltonian.

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

! LNPSI contains the log(psi) for each hamiltonian.

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

! LnPi contains the log(pi).

real, intent(out)										:: LnPi

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

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

! 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

! EXPX contains the value of exp( i*kx*x ) for a given kx and ion.

complex, dimension(0:Kmax, MaxIon), intent(out)			:: EXPX
complex, dimension(-Kmax:Kmax, MaxIon), intent(out)		:: EXPY, EXPZ

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

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

! SUMQX is the summation of qi * xi for all ions in the box.

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

! ULJ is the LJ energy of the system without the long range correction.
! UFOURIER is the coulombic fourier energy of the system.
! UREAL is the coulombic real energy of the system.
! USURF is the coulombic surface energy of the system.
! USELF_CH is the summation of the square of all the charges.
! USELF_MOL is the self energy of a given molecule.
! ULJ_MOL is the non bonded intramolecular LJ energy of a molecule.
! UION_MOL is the non bonded intramolecular ionic energy of a molecule.
! UINTRA contains the bending and torsional energy of each molecule.

real, dimension(Nham), intent(inout)					:: ULJBOX, ULJLR
real, dimension(Nham), intent(inout)					:: UREAL, UFOURIER
real, dimension(Nham), intent(inout)					:: USURF, USELF_CH
real, dimension(MaxMol,Nham), intent(inout)				:: ULJ_MOL
real, dimension(MaxMol,Nham), intent(inout)				:: USELF_MOL
real, dimension(MaxMol,Nham), intent(inout)				:: UION_MOL
real, dimension(MaxMol), intent(inout)					:: UINTRA

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

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

! 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 stuff

integer											:: h, i, j, k, m, n
integer											:: ljb, ionb
integer											:: mol2, Nbeads
integer											:: mol,  bead, sp
integer											:: stlj, stion
integer											:: llj, lion
integer											:: endlj, endion
integer											:: Count
integer											:: ljbk, ljbm
integer											:: ionbk, ionbm
integer											:: ntemp1, ntemp2
integer											:: nbendp

integer, dimension(MaxBeads, MaxSp)				:: LJBEAD, IONBEAD
integer, dimension(Nljgrs)						:: NGROUPS
integer, dimension(0:MaxSp)						:: Nmol_new
integer, dimension(MaxLJ)						:: TYPElj_tmp
integer, dimension(MaxIon)						:: TYPEion_tmp

integer, allocatable, dimension(:)				:: TYPElj_new
integer, allocatable, dimension(:)				:: TYPEion_new

logical											:: Ions	= .False.
logical											:: New

real											:: x, y, z
real											:: Uintra_new
real											:: Uvib, Ubend, Utor
real											:: Random
real											:: CoulCombo
real											:: Largest
real											:: rtrial
real											:: Utemp

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

real, dimension(20)								:: Xc, Yc, Zc
real, dimension(4)								:: ZERO2 = 0.0
real, dimension(Nham)							:: ULJBOX_new
real, dimension(Nham)							:: ULJLR_new
real, dimension(Nham)							:: ULJ_MOL_new
real, dimension(Nham)							:: UREAL_new
real, dimension(Nham)							:: USURF_new
real, dimension(Nham)							:: UFOURIER_new
real, dimension(Nham)							:: USELF_CH_new
real, dimension(Nham)							:: USELF_MOL_new
real, dimension(Nham)							:: UION_MOL_new
real, dimension(Nham)							:: ULJ
real, dimension(Nham)							:: UION
real, dimension(Nham)							:: UTOT
real, dimension(Nham)							:: dULJBOX
real, dimension(Nham)							:: dUREAL
real, dimension(Nham)							:: dUSURF
real, dimension(Nham)							:: dUFOURIER
real, dimension(Nham)							:: dUSELF_CH
real, dimension(Nham)							:: ULJ_MOL_part
real, dimension(Nham)							:: UION_MOL_part
real, dimension(Nham)							:: SUMQX_new 
real, dimension(Nham)							:: SUMQY_new
real, dimension(Nham)							:: SUMQZ_new
real, dimension(Nham)							:: BIGW_new, BIGW_old
real, dimension(Nsp)							:: PROB_SP
real, dimension(MaxLJ)							:: Xlj_tmp, Ylj_tmp, Zlj_tmp
real, dimension(MaxIon)							:: Xion_tmp, Yion_tmp, Zion_tmp
real, dimension(MaxLJ)							:: DAMPlj2_tmp, DAMPlj3_tmp 
real, dimension(MaxIon)							:: DAMPion_tmp
real, dimension(Nljgrs, Nljgrs, Nham)			:: EPS_CP

real, allocatable, dimension(:)					:: Xlj_new, Ylj_new, Zlj_new
real, allocatable, dimension(:)					:: Xion_new, Yion_new, Zion_new
real, allocatable, dimension(:)					:: DAMPlj2_new, DAMPlj3_new
real, allocatable, dimension(:)					:: DAMPion_new

complex, dimension(0:Kmax, MaxIon)				:: EXPX_tmp
complex, dimension(-Kmax:Kmax, MaxIon)			:: EXPY_tmp, EXPZ_tmp
complex, dimension(Nkvec, Nham)					:: SUMQEXPV_new

complex, allocatable, dimension(:,:)			:: EXPX_new, EXPY_new, EXPZ_new
complex, allocatable, dimension(:,:)			:: CZERO1, CZERO2


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

do i = 1, Nsp

	if( LENION(i) > 0 ) Ions = .True.

end do

LJBEAD = 0
IONBEAD = 0

ljb = 0
ionb = 0

do j = 1, Nsp

	ljb = 0
	ionb = 0

	do i = 1, LENLJ(j) + LENION(j)
	
		if(	BEADTYPE(i,j) == 'LJ' ) then

			ljb = ljb + 1
		
			LJBEAD(i,j) = ljb
		
			IONBEAD(i,j) = ionb
		
		else if( BEADTYPE(i,j) == 'ION' ) then

			ionb = ionb + 1

			IONBEAD(i,j) = ionb

			LJBEAD(i,j) = ljb

		end if

	end do

end do

DAMPlj2 = 1.0
DAMPlj3 = 1.0
DAMPion = 1.0

if( FromDisk ) then

	open( 30, file = InputConf )

	do i = 1, 15 + Nsp + Nham + 2*Nsp + sum(EESTEPS(1:NSp))
		
		read(30,*)

	end do

	Nlj = 0
	Nion = 0

	mol2 = 0

	Nbeads = sum( Nmol(1:Nsp) * ( LENLJ(1:Nsp) + LENION(1:Nsp) ) )
	
	do i = 1, Nbeads

		read(30,*) x, y, z, mol, bead, sp

		if( mol /= mol2 ) then

			SPECIES(mol) = sp
			
			STARTlj(mol) = Nlj + 1
			STARTion(mol) = Nion + 1

			LENGTHlj(mol) = LENLJ(sp)
			LENGTHion(mol) = LENION(sp)

			mol2 = mol

		end if

		if( BEADTYPE(bead,sp) == 'LJ' ) then

			Nlj = Nlj + 1

			Xlj( Nlj ) = x
			Ylj( Nlj ) = y
			Zlj( Nlj ) = z

			TYPElj( Nlj ) = GROUPTYPE(bead, sp)

			if( mol == tag_mol(sp) ) then

				DAMPlj2( Nlj ) = sqrt( BEADDAMP(bead, tag_val(sp), sp) )
				DAMPlj3( Nlj ) = BEADDAMP(bead, tag_val(sp), sp) ** (1.0/3.0)

			end if

		else if( BEADTYPE(bead,sp) == 'ION' ) then

			Nion = Nion + 1

			Xion( Nion ) = x
			Yion( Nion ) = y
			Zion( Nion ) = z

			TYPEion( Nion ) = GROUPTYPE(bead, sp)

			if( mol == tag_mol(sp) ) then

				DAMPion( Nion ) = sqrt( BEADDAMP(bead, tag_val(sp), sp) )

			end if

		end if

	end do

	SUMQX = 0.0
	SUMQY = 0.0
	SUMQZ = 0.0

	ULJBOX = 0.0
	ULJLR = 0.0
	ULJ_MOL = 0.0
	UINTRA = 0.0

	UREAL = 0.0
	UFOURIER = 0.0
	USELF_CH = 0.0
	USURF = 0.0
	USELF_MOL = 0.0
	UION_MOL = 0.0

	EXPX = ( 0.0, 0.0 )
	EXPY = ( 0.0, 0.0 )
	EXPZ = ( 0.0, 0.0 )

	SUMQEXPV = ( 0.0, 0.0 )

	call e6interact( Nlj, Xlj, Ylj, Zlj,  &
					 TYPElj, DAMPlj2, DAMPlj3, &
					 Nmol(0), LENGTHlj, &
					 Nham, Nljgrs, EPS, SIG, CP, &
					 ALP, RMAX, BoxSize, ULJBOX )

	NGROUPS = 0
	
	do j = 1, Nlj

		NGROUPS( TYPElj(j) ) = NGROUPS( TYPElj(j) ) + 1

	end do

	do h = 1, Nham

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

			EPS_CP(:,:,h) = EPS(:,:,h) * CP(:,:,h) * ALP(:,:,h) / &
							( ALP(:,:,h) - 6.0 )	/ 4.0

		else

			EPS_CP(:,:,h) = EPS(:,:,h) * ALP(:,:,h) / ( ALP(:,:,h) - 6.0 ) * &
							( ALP(:,:,h) / 6.0 ) ** ( 6.0 / ( ALP(:,:,h) - 6.0 ) ) / 4.0

		end if

	end do

	call lrcorr( Nljgrs, Nham, BoxSize, NGROUPS, EPS_CP, SIG, XLRCORR, &
				 ELRCORR, ULJLR )

	do j = 1, Nmol(0)

		sp = SPECIES(j)
		stlj = STARTlj(j)

		do k = 1, LENGTHlj(j) + LENGTHion(j) - 1
				
			do m = k + 1, LENGTHlj(j) + LENGTHion(j)

				if( BONDSAPART(k,m,sp) >= 3 .AND. BEADTYPE(k,sp) == 'LJ' .AND. &
					BEADTYPE(m,sp) == 'LJ') then

					ljbk = stlj + LJBEAD(k,sp) - 1
					ljbm = stlj + LJBEAD(m,sp) - 1

					call e6molecule( 1, Xlj(ljbm), Ylj(ljbm), &
									 Zlj(ljbm), TYPElj(ljbm), &
									 DAMPlj2(ljbm), DAMPlj3(ljbm), &
									 1, Xlj(ljbk), Ylj(ljbk), &
									 Zlj(ljbk), TYPElj(ljbk), &
									 DAMPlj2(ljbk), DAMPlj3(ljbk), &
									 Nham, Nljgrs, EPS, SIG, CP, &
									 ALP, RMAX, BoxSize, &
									 ULJ_MOL_part )

					if( BONDSAPART(k,m,sp) == 3 ) ULJ_MOL_part = f14_lj * ULJ_MOL_part 

					ULJ_MOL(j,:) = ULJ_MOL(j,:) + ULJ_MOL_part(:)

				end if

			end do

		end do

	end do

	if( Ions ) then

		call Surf_Move( Nion, Xion, Yion, Zion, &
					    TYPEion, DAMPion, &
						Nham, Niongrs, CHARGE, BoxSize, &
					    SUMQX, SUMQY, SUMQZ, &
						SUMQX_new, SUMQY_new, SUMQZ_new, &
						USURF)

		USURF = USURF * CoulCombo

		SUMQX = SUMQX_new
		SUMQY = SUMQY_new
		SUMQZ = SUMQZ_new