add.f90

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

Subroutine Add( Nlj, Nion, MaxMol, MaxSp, MaxBeads, MaxEE, Nmol, &
				Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
				Xion, Yion, Zion, TYPEion, DAMPion, &
				MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
				METHOD, MaxInt, MaxReal, INTPARAM, &
				REALPARAM, BEADTYPE, BEADDAMP, BONDSAPART, &
				ERSTEPS, ERSTART, EREND, EESTEPS, EEW, &
				BoxSize, SPECIES, LENGTHlj, LENGTHion, &
				STARTlj, STARTion, &
				LENLJ, LENION, GROUPTYPE_LJ, GROUPTYPE_ION, &
				Nham, BETA, ZETA, LNW, LNPSI, LnPi, XLRCORR, ELRCORR, &
				Nljgrs, EPS, SIG, CP, ALP, RMAX, Niongrs, CHARGE, &
				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, SpID, &
				Success, tag_val, tag_mol, &
				Nres, Nresmol, reslen, Xresmols, &
				Yresmols, Zresmols, Uint_resm, &
				Ulj_resm, Uion_resm, &
				f14_lj, f14_ion, ER_HIST, Seed )

implicit none

! This subroutine regrows a flexible molecule from StartGrowthStep to the 
! last step and then accepts or rejects the move.

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

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

! MaxMol is the maximum number of molecules per box.
! MaxSp is the maximum number of Species in the simulation.
! MaxBeads in the maximum number of beads per molecule.

integer, intent(in)										:: MaxMol
integer, intent(in)										:: MaxSp
integer, intent(in)										:: MaxBeads
integer, intent(in)										:: MaxEE

! Nmol is the number of molecules in the simulation boxes.
! Nmol(0, iphase) ==> total number of molecules in that phase.
! Nmol(1:MaxSp, iphase) ==> number of molecules of that species in that phase.

integer, dimension(0:MaxSp), intent(inout)				:: 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(Nlj+MaxBeads), intent(inout)			:: Xlj, Ylj, Zlj
real, dimension(Nion+MaxBeads), intent(inout)			:: Xion, Yion, Zion

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

integer, dimension(Nlj+MaxBeads), intent(inout)			:: TYPElj
integer, dimension(Nion+MaxBeads), intent(inout)		:: TYPEion

real, dimension(Nlj+MaxBeads), intent(inout)			:: DAMPlj2, DAMPlj3
real, dimension(Nion+MaxBeads), intent(inout)			:: DAMPion

! MaxSteps is the maximum number of CB steps.

integer, intent(in)										:: MaxSteps

! 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(MaxSteps, MaxSp), intent(in)			:: STEPSTART
integer, dimension(MaxSteps, MaxSp), intent(in)			:: STEPLENGTH
integer, dimension(MaxSteps, MaxSp), intent(in)			:: NTRIALS

! 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'.
! BONDSAPART gives the bondlengths any two LJ beads are away from each other.

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

integer, dimension(MaxSp), intent(in)					:: ERSTEPS
integer, dimension(MaxSteps,MaxSp), intent(in)			:: ERSTART
integer, dimension(MaxSteps,MaxSp), intent(in)			:: EREND

integer, dimension(MaxSp), intent(in)					:: EESTEPS
real, dimension(MaxEE, MaxSp), intent(in)				:: EEW

! BoxSize is the length of the simulation box.

real, intent(in)										:: BoxSize

! 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(Nmol(0)+1), intent(inout)			:: SPECIES
integer, dimension(Nmol(0)+1), intent(inout)			:: LENGTHlj, LENGTHion
integer, dimension(Nmol(0)+1), intent(inout)			:: STARTlj, STARTion

! LENLJ contains the number of LJ beads of each species.
! LENION contains the number of ionic beads of each species.

integer, dimension(MaxSp), intent(in)					:: LENLJ, LENION

! GROUPTYPE_LJ contains the group identity of the LJ beads for a given species.
! GROUPTYPE_ION contains the group identity of the ionic beads for a given species.

integer, dimension(MaxBeads,MaxSp), intent(in)			:: GROUPTYPE_LJ, GROUPTYPE_ION

! Nham is the number of hamiltonians being used.

integer, intent(in)										:: Nham

! beta contains 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

! LNW contains the weight of each hamiltonian.

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

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

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

! LnPi contains the log(pi).

real, intent(inout)										:: LnPi

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

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

! 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

! Niongrs is the number of ionic groups in the system.
! CHARGE is a rank 2 array containing the charge of group i for each 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

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

complex, dimension(0:Kmax, Nion+MaxBeads), intent(inout)		:: EXPX
complex, dimension(-Kmax:Kmax, Nion+MaxBeads), intent(inout)	:: 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)				:: USELF_MOL
real, dimension(MaxMol,Nham), intent(inout)				:: ULJ_MOL
real, dimension(MaxMol,Nham), intent(inout)				:: UION_MOL
real, dimension(MaxMol), intent(inout)					:: UINTRA

! SpID indicates which type of species was used in the move.

integer													:: SpID

! Success is a logical indicating whether the move was successful or not.

logical, intent(out)									:: Success

integer													:: tag_val
integer													:: 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

integer, dimension(MaxSteps+1,MaxSp,2), intent(inout)	:: ER_HIST

! Seed is the current random number generator seed value.

integer, intent(inout)									:: Seed

real, external											:: ran2

! Local Variables

integer										:: i, j
integer										:: llj, lion
integer										:: ljbk, ionbk
integer, allocatable, dimension(:) 			:: TYPElj_new
integer, allocatable, dimension(:) 			:: TYPEion_new

logical										:: New

real										:: Largest
real										:: LnPi_tmp, LnPi_old
real										:: Uintra_new
real										:: CoulCombo
real, dimension(Nham)						:: BIGW_new
real, dimension(Nham)						:: TMP, dU
real, dimension(Nham)						:: ULJBOX_new, ULJLR_new
real, dimension(Nham)						:: UREAL_new, UFOURIER_new
real, dimension(Nham)						:: USURF_new, USELF_CH_new
real, dimension(Nham)						:: USELF_MOL_new
real, dimension(Nham)						:: ULJ_MOL_new
real, dimension(Nham)						:: UION_MOL_new
real, dimension(Nham)						:: SUMQX_new, SUMQY_new, SUMQZ_new
real, dimension(Nham)						:: ZETA_tmp

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

real, parameter								:: Pi = 3.14159265359
real, parameter								:: ec = 1.60217733e-19
real, parameter								:: eps0 = 8.854187817e-12
real, parameter								:: kB = 1.380658e-23
										
complex, dimension(Nkvec, Nham)				:: SUMQEXPV_new
complex, allocatable, dimension(:,:)		:: EXPX_new, EXPY_new, EXPZ_new


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