alpha_ch.f90

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

subroutine alpha_ch( MaxSp, MaxMol, MaxBeads, MaxEE, Nmol, &
					 Nlj, Nion, Xlj, Ylj, Zlj, &
					 TYPElj, DAMPlj2, DAMPlj3, &
					 Xion, Yion, Zion, &
					 TYPEion, DAMPion, BoxSize, &
					 LENGTHlj, LENGTHion, SPECIES, &
					 STARTlj, STARTion, &
					 BEADTYPE, BEADDAMP, BONDSAPART, &
					 EESTEPS, EEW, &
					 Nham, BETA, ZETA, LNW, LNPSI, LnPi, &
					 Nljgrs, EPS, SIG, CP, ALP, RMAX, Niongrs, &
					 CHARGE, Alpha, Kmax, Nkvec, KX, KY, KZ, &
					 CONST, EXPX, EXPY, EXPZ, SUMQEXPV, &
					 SUMQX, SUMQY, SUMQZ, ULJBOX, ULJ_MOL, & 
					 UFOURIER, UREAL, USURF, USELF_CH, &
					 USELF_MOL, UION_MOL, &
					 tag_val, tag_mol, updown, &
					 f14_lj, f14_ion, &
					 SpID, Success, Seed )

implicit none

! Nmol is the number of molecules in the simulation box.
! MaxSp is the maximum number of species in the system.
! Nlj is the number of LJ beads in the simulation box.
! Nion is the number of ionic beads in the simulation box.

integer, intent(in)									:: MaxSp
integer, intent(in)									:: MaxMol
integer, intent(in)									:: MaxBeads
integer, intent(in)									:: MaxEE
integer, dimension(0:MaxSp), intent(in)				:: Nmol
integer, intent(in)									:: Nlj
integer, intent(in)									:: Nion

! 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), intent(inout)					:: Xlj, Ylj, Zlj
real, dimension(Nion), 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), intent(in)					:: TYPElj
integer, dimension(Nion), intent(in)				:: TYPEion

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

! BoxSize is the length of the simulation box.

real, intent(in)									:: BoxSize

! LENGTHlj contains the number of LJ beads in each molecule.
! LENGTHion	contains the number of ionic beads in each molecule.
! SPECIES contains the species identity of 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)), intent(in)				:: LENGTHlj, LENGTHion
integer, dimension(Nmol(0)), intent(in)				:: SPECIES
integer, dimension(Nmol(0)), intent(in)				:: STARTlj, STARTion

character*5, dimension(MaxBeads, MaxSp)				:: BEADTYPE
real, dimension(MaxBeads, MaxEE, MaxSp)				:: BEADDAMP

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

integer, dimension(MaxBeads,MaxBeads,MaxSp)			:: BONDSAPART

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

! Nham is the number of hamiltonians being used.

integer, intent(in)									:: Nham

! beta contains the reciprical temperature.

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

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

! 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(inout)				:: LNPSI

! LnPi contains the log(pi).

real, intent(inout)									:: LnPi

! 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), intent(inout)		:: EXPX
complex, dimension(-Kmax:Kmax, Nion), 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

! ULJBOX 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.

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

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

integer												:: tag_val, tag_mol

integer												:: updown

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

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

! SpID gives the identity of the species that was displaced or rotated.

integer												:: SpID

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

logical, intent(out)								:: Success

! Seed is the current random number generator seed value.

integer, intent(inout)								:: Seed 

real, external										:: ran2

! Local Variables.

integer												:: h, i, j, k, m, Count
integer												:: lenlj, stlj, endlj
integer												:: lenion, stion, endion
integer												:: ljbk, ljbm, ionbk, ionbm
integer, dimension( Nlj )							:: temp4
integer, dimension( Nion )							:: temp8
integer, dimension( MaxBeads )						:: LJBEAD, IONBEAD

real												:: CoulCombo
real												:: Largest, LnPi_new
real												:: LnPi_tmp
real												:: max_damp

real, allocatable, dimension( : )					:: Xlj_old, Ylj_old, Zlj_old
real, allocatable, dimension( : )					:: DAMPlj2_old, DAMPlj3_old
real, allocatable, dimension( : )					:: DAMPlj2_new, DAMPlj3_new
real, dimension( Nlj )								:: temp1, temp2, temp3
real, dimension( Nion )								:: temp5, temp6, temp7
real, dimension( Nlj )								:: temp9, temp10
real, dimension( Nion )								:: temp11
real, allocatable, dimension( : )					:: Xion_old, Yion_old, Zion_old
real, allocatable, dimension( : )					:: DAMPion_old, DAMPion_new
real, allocatable, dimension( : )					:: dDAMPion
real, dimension(Nham)								:: LNPSI_new
real, dimension(Nham)								:: dULJBOX
real, dimension(Nham)								:: dULJ_MOL
real, dimension(Nham)								:: dUREAL
real, dimension(Nham)								:: dUSURF
real, dimension(Nham)								:: dUFOURIER
real, dimension(Nham)								:: dUSELF_CH
real, dimension(Nham)								:: dUSELF_MOL
real, dimension(Nham)								:: dUION_MOL
real, dimension(Nham)								:: dU
real, dimension(Nham)								:: U_part
real, dimension(Nham)								:: ZETA_tmp
real, dimension(Nham)								:: SUMQX_NEW, SUMQY_NEW, SUMQZ_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


Success = .false.

if( tag_mol == 0 ) then

	if( Nmol(SpID) == 0 ) then
	
		tag_val = 0
		return

	end if

	tag_mol = int( Nmol(SpID) * ran2(Seed) ) + 1

	i = 0
	Count = 0

	do while ( Count < tag_mol )
		
		i = i + 1
		
		if( SPECIES(i) == SpID ) Count = Count + 1

	end do

	tag_mol = i

end if

lenlj = LENGTHlj( tag_mol )
stlj = STARTlj( tag_mol )
endlj = stlj + lenlj - 1

lenion = LENGTHion( tag_mol )
stion = STARTion( tag_mol )
endion = stion + lenion - 1

allocate( Xlj_old(lenlj) )
allocate( Ylj_old(lenlj) )
allocate( Zlj_old(lenlj) )

allocate( DAMPlj2_old(lenlj) )
allocate( DAMPlj3_old(lenlj) )

allocate( DAMPlj2_new(lenlj) )
allocate( DAMPlj3_new(lenlj) )

Xlj_old( 1:lenlj ) = Xlj( stlj:endlj )
Ylj_old( 1:lenlj ) = Ylj( stlj:endlj )
Zlj_old( 1:lenlj ) = Zlj( stlj:endlj )

DAMPlj2_old( 1:lenlj ) = DAMPlj2( stlj:endlj )
DAMPlj3_old( 1:lenlj ) = DAMPlj3( stlj:endlj )

! Set DAMPlj2_new values later

if( lenion > 0 ) then

	allocate( Xion_old(lenion) )
	allocate( Yion_old(lenion) )
	allocate( Zion_old(lenion) )

	allocate( DAMPion_old(lenion) )
	allocate( DAMPion_new(lenion) )
	allocate( dDAMPion(lenion) )

	Xion_old( 1:lenion ) = Xion( stion:endion )
	Yion_old( 1:lenion ) = Yion( stion:endion )
	Zion_old( 1:lenion ) = Zion( stion:endion )

	DAMPion_old( 1:lenion ) = DAMPion( stion:endion )

	! Set DAMPion_new values later

end if

LJBEAD = 0
IONBEAD = 0

ljbk = 0
ionbk = 0

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

		ljbk = ljbk + 1
		
		LJBEAD(i) = ljbk

		IONBEAD(i) = ionbk
		
		DAMPlj2_new( ljbk ) = sqrt( BEADDAMP( i, tag_val + updown, SpID ) )
		DAMPlj3_new( ljbk ) = BEADDAMP( i, tag_val + updown, SpID )	** (1.0/3.0)

	else if( BEADTYPE(i,SpID) == 'ION' ) then

		ionbk = ionbk + 1

		IONBEAD(i) = ionbk

		LJBEAD(i) = ljbk

		DAMPion_new( ionbk ) = sqrt( BEADDAMP( i, tag_val + updown, SpID ) )

	end if

end do

ljbk = 0
ionbk = 0

! Calculation of ULJBOX.

if( stlj == 1 )	then

	if( Nmol(0) == 1 ) then
	
		dULJBOX	= 0.0

	else

		call e6molecule2( lenlj, Xlj_old, Ylj_old, Zlj_old, TYPElj(1:lenlj), &
						  DAMPlj2_old, DAMPlj3_old, DAMPlj2_new, DAMPlj3_new, &
						  Nlj - lenlj, Xlj(endlj+1:Nlj), Ylj(endlj+1:Nlj), &						  
						  Zlj(endlj+1:Nlj), TYPElj(endlj+1:Nlj), &
						  DAMPlj2(endlj+1:Nlj),	DAMPlj3(endlj+1:Nlj), &
						  DAMPlj2(endlj+1:Nlj),	DAMPlj3(endlj+1:Nlj), &
						  Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
						  BoxSize, dULJBOX )

	end if