gcmc.f90

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

Program gcmc

Use Portlib	    ! # MS Fortran
implicit none

! Things to Check:
! 1. Rmax for water simulations (BigRead.f90).
! 2. Molecule flip for n-alkanes (Regrow.f90).
! 3. Parameters below.

! MaxSp = The maximum number of species.
! MaxMol = The maximum number of total molecules.
! MaxLJ = The maximum number of LJ/exp-6 beads.
! MaxIon = The maximum number of point charges.
! MaxSteps = The maximum number of configurational-bias steps.
! MaxBeads = The maximum number of beads per molecule.
!            LJ/exp-6 + Coulombic
! MaxInt = The maximum number of integer parameters required for 
!          any given growth method used.
! MaxReal = The maximum number of real parameters required for 
!           any given growth method used.
! MaxEE = The maximum number of expanded ensemble steps used. 
! MaxResmv = The maximum number of reservoir moves used. 
! MaxResp = The maximum number of reservoir parameters for 
!           any given reservoir move used.
! MaxVib = The maximum number of vibration potentials for a
!          given reservoir structure. 
! MaxBend = The maximum number of bending potentials for a
!           given reservoir structure. 
! MaxTor = The maximum number of torsion potentials for a
!          given reservoir structure. 

! Example: Harris and Yung model for carbon dioxide.
! Minimum parameters for a simulation with a maximum
! of 200 particles.
! MaxSp = 1, MaxMol = 200, MaxLJ = 600, MaxIon = 600
! MaxSteps = 3, MaxBeads = 6, MaxInt = 3, MaxReal = 3
! MaxEE = 1

integer, parameter								:: MaxSp = 1
integer, parameter								:: MaxMol = 200
integer, parameter								:: MaxLJ = 1600
integer, parameter								:: MaxIon = 1
integer, parameter								:: MaxSteps = 8
integer, parameter								:: MaxBeads = 8
integer, parameter								:: MaxInt = 7
integer, parameter								:: MaxReal = 8
integer, parameter								:: MaxEE = 1

integer, parameter								:: MaxResmv = 20
integer, parameter								:: MaxResp = 10
integer, parameter								:: MaxVib = 1
integer, parameter								:: MaxBend = 10
integer, parameter								:: MaxTor = 10

integer, parameter								:: Kmax	= 5
integer, parameter								:: Maxkvec = 276 
integer, parameter								:: nmoves = 3 
integer, parameter								:: nmtype = 7 

! for Kmax = 5, Maxkvec = 276
! for Kmax = 6, Maxkvec = 501
! for Kmax = 7, Maxkvec = 754

real											:: f14_lj = 0.0
real											:: f14_ion = 0.0

integer											:: Nham
integer											:: Nlj = 0
integer											:: Nion = 0
integer											:: Nljgrs
integer											:: Niongrs
integer											:: Nsp
integer											:: Seed
integer											:: StartConf
integer											:: Nequil
integer											:: Nprod
integer											:: NinCycle
integer											:: Nstorage
integer											:: Nadjust
integer											:: Ncollect
integer											:: Nstages
integer											:: Neew
integer											:: Nresmol
integer											:: Nrefresh
integer											:: Nrefrmoves
integer											:: Nkvec = Maxkvec
integer											:: Istore
integer											:: Moves
integer											:: h
integer											:: i
integer											:: j
integer											:: k
integer											:: m
integer											:: SpID
integer											:: DispOrRot
integer											:: EquilOrProd
integer											:: stc
integer											:: Ubins = 0
integer											:: ljb
integer											:: ionb
integer											:: Wcounts
integer											:: Stage = 1
integer											:: nentry = 0
integer											:: rmol
integer											:: Nres

integer, dimension(MaxSp)						:: nt1
integer, dimension(MaxSp)						:: tag_val
integer, dimension(MaxSp)						:: tag_mol
integer, dimension(MaxSp)						:: nin = 0
integer, dimension(MaxSp)						:: nout = 0
integer, dimension(MaxSp)						:: Nsim_max = 0
integer, dimension(MaxSp)						:: Nsim_min	= MaxMol
integer, dimension(MaxSp)						:: NSTEPS = 0
integer, dimension(MaxSp)						:: LENLJ = 0
integer, dimension(MaxSp)						:: LENION = 0
integer, dimension(MaxBeads,MaxSp)				:: GROUPTYPE = 0
integer, dimension(MaxSteps,MaxSp)				:: NTRIALS = 0
integer, dimension(MaxSteps,MaxSp)				:: STEPSTART = 0
integer, dimension(MaxSteps,MaxSp)				:: STEPLENGTH = 0
integer, dimension(MaxInt,MaxBeads,MaxSp)		:: INTPARAM = 0
integer, dimension(MaxBeads,MaxBeads,MaxSp)		:: BONDSAPART = 0
integer, dimension(0:MaxSp)						:: Nmol = 0
integer, dimension(Maxkvec)						:: KX = 0
integer, dimension(Maxkvec)						:: KY = 0
integer, dimension(Maxkvec)						:: KZ = 0
integer, dimension(MaxLJ)						:: TYPElj = 0
integer, dimension(MaxIon)						:: TYPEion = 0
integer, dimension(MaxMol)						:: SPECIES = 0
integer, dimension(MaxMol)						:: LENGTHlj = 0
integer, dimension(MaxMol)						:: LENGTHion = 0
integer, dimension(MaxMol)						:: STARTlj = 0
integer, dimension(MaxMol)						:: STARTion = 0
integer, dimension(nmtype,MaxSp)				:: NATT = 0
integer, dimension(nmtype,MaxSp)				:: NSUC = 0
integer, dimension(nmtype,MaxSp)				:: TOT_NATT = 0
integer, dimension(nmtype,MaxSp)				:: TOT_NSUC = 0
integer, dimension(nmtype)						:: RES_ATT = 0
integer, dimension(nmtype)						:: RES_SUC = 0
integer, dimension(MaxBeads,MaxSp)				:: GROUPTYPE_LJ
integer, dimension(MaxBeads,MaxSp)				:: GROUPTYPE_ION
integer, dimension(MaxSteps+1,MaxSp,2)			:: ER_HIST
integer, dimension(MaxSp)						:: ERSTEPS
integer, dimension(MaxSteps,MaxSp)				:: ERSTART
integer, dimension(MaxSteps,MaxSp)				:: EREND
integer, dimension(MaxSp)						:: EESTEPS
integer, dimension(0:MaxEE,MaxSp)				:: ALP_PROB = 0
integer, dimension(0:180)						:: th_h	= 0
integer, dimension(-180:180)					:: phi_h = 0

integer, allocatable, dimension(:)				:: NST_TOT
integer, allocatable, dimension(:)				:: NST_WT
integer, allocatable, dimension(:)				:: STOP_H
integer, allocatable, dimension(:,:)			:: Npart
integer, allocatable, dimension(:)				:: Nresmv
integer, allocatable, dimension(:)				:: reslen
integer, allocatable, dimension(:)				:: Nvib
integer, allocatable, dimension(:)				:: Nbend
integer, allocatable, dimension(:)				:: Ntor
integer, allocatable, dimension(:,:,:)			:: RESPARAM
integer, allocatable, dimension(:,:,:)			:: IVIB
integer, allocatable, dimension(:,:,:)			:: IBEND
integer, allocatable, dimension(:,:,:)			:: ITOR

logical											:: FromDisk
logical											:: Success
logical											:: New
logical											:: exists

character*30									:: InputFile
character*30									:: InitialConf
character*30									:: FinalConf
character*30									:: PlotsFile
character*30									:: ResultsFile
character*30									:: UN_HistFile
character*30									:: WFile
character*30									:: filename
character*30, dimension(0:MaxSp)				:: N_HistFile
character*30									:: InFile

character*80									:: Timestring
character*80									:: Datestring
character*80									:: Hostname
character*80									:: Infostring

character*10									:: histtype

character*5, dimension(MaxBeads, MaxSp)			:: BEADTYPE
character*10, dimension(MaxBeads, MaxSp)		:: METHOD
character*15, dimension(MaxSp)					:: NAMEsp

character*15, allocatable, dimension(:)			:: NAMElj
character*15, allocatable, dimension(:)			:: NAMEion
character*10, allocatable, dimension(:,:)		:: RESMOVE
					
real, parameter									:: kB = 1.380658e-23
real, parameter									:: Nav = 6.02214e23
real, parameter									:: SafetyFactor = 0.45

real											:: Alpha = 5.6
real											:: tt
real											:: tstart
real											:: tot_time
real											:: Ranq
real											:: Ranq2
real											:: BoxSize
real											:: Volume
real											:: Umin
real											:: Umax
real											:: Urange
real											:: Uwidth
real											:: Usim_min = 1.0e10
real											:: Usim_max	= -1.0e10
real											:: estart
real											:: LnPi
real											:: Largest

!real, dimension(2)								:: xxx			! # U_ALPHA
real, dimension(nmoves)							:: PROB_MOVE
real, dimension(2)								:: PROB_DR
real, dimension(MaxSp)							:: PROB_SP_DR
real, dimension(MaxSp)							:: PROB_SP_CD
real, dimension(MaxSp)							:: PROB_SP_RG
real, dimension(MaxSp)							:: MASSsp  = 0.0
real, dimension(MaxReal,MaxBeads,MaxSp)			:: REALPARAM = 0.0
real, dimension(MaxBeads,MaxEE,MaxSp)			:: BEADDAMP
real, dimension(MaxSp)							:: DXYZ
real, dimension(MaxSp)							:: DROT
real, dimension(605)							:: XLRCORR
real, dimension(605)							:: ELRCORR
real, dimension(Maxkvec)						:: CONST
real, dimension(MaxLJ)							:: Xlj
real, dimension(MaxLJ)							:: Ylj
real, dimension(MaxLJ)							:: Zlj
real, dimension(MaxIon)							:: Xion
real, dimension(MaxIon)							:: Yion
real, dimension(MaxIon)							:: Zion
real, dimension(MaxLJ)							:: DAMPlj2
real, dimension(MaxLJ)							:: DAMPlj3
real, dimension(MaxIon)							:: DAMPion
real, dimension(MaxMol)							:: UINTRA = 0.0
real, dimension(nmtype)							:: MOVETIME = 0.0
real, dimension(MaxEE,MaxSp)					:: EEW

real, allocatable, dimension(:)					:: BETA
real, allocatable, dimension(:,:)				:: ZETA
real, allocatable, dimension(:)					:: SUMQX 
real, allocatable, dimension(:)					:: SUMQY 
real, allocatable, dimension(:)					:: SUMQZ 
real, allocatable, dimension(:)					:: ULJBOX 
real, allocatable, dimension(:)					:: ULJLR 
real, allocatable, dimension(:,:)				:: ULJ_MOL 
real, allocatable, dimension(:)					:: UREAL 
real, allocatable, dimension(:)					:: USURF 
real, allocatable, dimension(:)					:: UFOURIER	
real, allocatable, dimension(:)					:: USELF_CH
real, allocatable, dimension(:,:)				:: USELF_MOL 
real, allocatable, dimension(:,:)				:: UION_MOL 
real, allocatable, dimension(:)					:: UION	
real, allocatable, dimension(:)					:: ULJ 
real, allocatable, dimension(:)					:: UTOT	
real, allocatable, dimension(:)					:: LNW
real, allocatable, dimension(:)					:: LNW_old
real, allocatable, dimension(:)					:: LNPSI
real, allocatable, dimension(:)					:: PSI_PI
real, allocatable, dimension(:)					:: AVG_PSI_PI
real, allocatable, dimension(:)					:: PSI_PI_AVG
real, allocatable, dimension(:)					:: U_AVG
real, allocatable, dimension(:,:,:)				:: N_HIST
real, allocatable, dimension(:,:,:)				:: UN_HIST

real, allocatable, dimension(:,:,:)				:: EPS
real, allocatable, dimension(:,:,:)				:: SIG
real, allocatable, dimension(:,:,:)				:: CP
real, allocatable, dimension(:,:,:)				:: ALP
real, allocatable, dimension(:,:,:)				:: RMAX
real, allocatable, dimension(:,:,:)				:: KAPPA
real, allocatable, dimension(:,:,:)				:: LAMDA
real, allocatable, dimension(:,:,:)				:: GAMMA
real, allocatable, dimension(:)					:: MASSlj
real, allocatable, dimension(:,:)				:: CHARGE
real, allocatable, dimension(:)					:: MASSion

real, allocatable, dimension(:)					:: Energy
real, allocatable, dimension(:,:)				:: PROB_RMV
real, allocatable, dimension(:,:)				:: Xres
real, allocatable, dimension(:,:)				:: Yres
real, allocatable, dimension(:,:)				:: Zres
real, allocatable, dimension(:)					:: UResInt
real, allocatable, dimension(:,:)				:: UResLj
real, allocatable, dimension(:,:)				:: UResIon
real, allocatable, dimension(:,:,:)				:: RVIB
real, allocatable, dimension(:,:,:)				:: RBEND
real, allocatable, dimension(:,:,:)				:: RTOR
real, allocatable, dimension(:,:,:)				:: Xresmols
real, allocatable, dimension(:,:,:)				:: Yresmols
real, allocatable, dimension(:,:,:)				:: Zresmols
real, allocatable, dimension(:,:)				:: Uint_resm
real, allocatable, dimension(:,:,:)				:: Ulj_resm
real, allocatable, dimension(:,:,:)				:: Uion_resm

real, external									:: ran2
!real, external									:: etime	! # U_ALPHA
															
complex, dimension(0:Kmax,MaxIon)				:: EXPX
complex, dimension(-Kmax:Kmax,MaxIon)			:: EXPY
complex, dimension(-Kmax:Kmax,MaxIon)			:: EXPZ

complex, allocatable, dimension(:,:)			:: SUMQEXPV

tstart = timef()				! # C_ALPHA
!tstart = etime(xxx)				! # U_ALPHA

write(*,"(A)") ' Enter Input, Initial Conf., FinalConf., Plots, &
				&Results, Energy-Particles Histogram, and &
				&Particle Histogram file names. '

read(*,*) InputFile
read(*,*) InitialConf
read(*,*) FinalConf
read(*,*) PlotsFile
read(*,*) ResultsFile
read(*,*) UN_HistFile
read(*,*) N_HistFile(0)
read(*,*) WFile

InputFile = trim(InputFile)//'.dat'
InitialConf = trim(InitialConf)//'.dat'
FinalConf = trim(FinalConf)//'.dat'
PlotsFile = trim(PlotsFile)//'.dat'
ResultsFile = trim(ResultsFile)//'.dat'
WFile = trim(WFile)//'.dat'

call LittleRead( InputFile, Nljgrs, Niongrs, NSp, Nham, Nstages, Nres )

allocate( EPS( Nljgrs, Nljgrs, Nham ) )
allocate( SIG( Nljgrs, Nljgrs, Nham ) )
allocate( CP( Nljgrs, Nljgrs, Nham ) )
allocate( ALP( Nljgrs, Nljgrs, Nham ) )
allocate( RMAX( Nljgrs, Nljgrs, Nham ) )
allocate( KAPPA( Nljgrs, Nljgrs, Nham ) )
allocate( LAMDA( Nljgrs, Nljgrs, Nham ) )
allocate( GAMMA( Nljgrs, Nljgrs, Nham ) )
allocate( MASSlj( Nljgrs ) )
allocate( NAMElj( Nljgrs ) )

allocate( CHARGE( Niongrs, Nham ) )
allocate( MASSion( Niongrs ) )
allocate( NAMEion( Niongrs ) )

allocate( BETA( Nham ) )
allocate( ZETA( MaxSp, Nham ) )
allocate( SUMQX( Nham ) )
allocate( SUMQY( Nham ) )
allocate( SUMQZ( Nham ) )
allocate( ULJBOX( Nham ) )
allocate( ULJLR( Nham ) )
allocate( ULJ_MOL( MaxMol,Nham ) )
allocate( UREAL( Nham ) )
allocate( USURF( Nham ) )
allocate( UFOURIER( Nham ) )
allocate( USELF_CH( Nham ) )
allocate( USELF_MOL( MaxMol,Nham ) )
allocate( UION_MOL( MaxMol,Nham ) )
allocate( UION( Nham ) )
allocate( ULJ( Nham ) )
allocate( UTOT( Nham ) )
allocate( LNW( Nham ) )
allocate( LNW_old( Nham ) )
allocate( LNPSI( Nham ) )
allocate( PSI_PI( Nham ) )
allocate( AVG_PSI_PI( Nham ) )
allocate( U_AVG( Nham ) )
allocate( N_HIST( 0:MaxMol,MaxSp,Nham ) )

allocate( PSI_PI_AVG( Nham ) )
allocate( NST_TOT( Nstages ) )
allocate( NST_WT( Nstages ) )
allocate( STOP_H( Nstages ) )

allocate( reslen(Nres) )
allocate( Xres(MaxBeads,Nres) )
allocate( Yres(MaxBeads,Nres) )
allocate( Zres(MaxBeads,Nres) )
allocate( Nresmv(Nres) )
allocate( PROB_RMV(MaxResmv,Nres) )
allocate( RESMOVE(MaxResmv,Nres) )
allocate( RESPARAM(MaxResp,MaxResmv,Nres) )
allocate( UResInt(Nres) )
allocate( UResLj(Nham,Nres) )
allocate( UResIon(Nham,Nres) )
allocate( Nvib(Nres) )
allocate( Nbend(Nres) )
allocate( Ntor(Nres) )

allocate( IVIB(2,MaxVib,Nres) )
allocate( IBEND(3,MaxBend,Nres) )
allocate( ITOR(4,MaxTor,Nres) )
allocate( RVIB(2,MaxVib,Nres) )
allocate( RBEND(2,MaxBend,Nres) )
allocate( RTOR(4,MaxTor,Nres) )

SUMQX = 0.0
SUMQY = 0.0
SUMQZ = 0.0
ULJBOX = 0.0
ULJLR = 0.0
ULJ_MOL = 0.0
UREAL = 0.0
USURF = 0.0
UFOURIER = 0.0
USELF_CH = 0.0
USELF_MOL = 0.0
UION_MOL = 0.0
UION = 0.0
ULJ = 0.0
UTOT = 0.0
AVG_PSI_PI = 0.0
PSI_PI_AVG = 0.0
U_AVG = 0.0
N_HIST = 0.0

ER_HIST = 0

call BigRead( InputFile, InitialConf, 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, Nsim_min, Nsim_max, LNW, &
			  Nstages, NST_TOT, NST_WT, STOP_H, Seed )

BoxSize = Volume ** ( 1.0 / 3.0 )

ZETA = exp(ZETA)

do j = 1, Nsp

	ljb = 0
	ionb = 0

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

			ljb = ljb + 1
		
			GROUPTYPE_LJ(ljb,j) = GROUPTYPE(i,j)
		
		else if( BEADTYPE(i,j) == 'ION' ) then

			ionb = ionb + 1

			GROUPTYPE_ION(ionb,j) = GROUPTYPE(i,j)

		end if

	end do

end do

do i = 1, Nres

	InFile = 'reserv_'//char(48+i)//'.dat'

	call ResRead( InFile, MaxBeads, reslen(i), &
					Xres(:,i), Yres(:,i), Zres(:,i), &
					Nresmv(i), MaxResmv, MaxResp, &
					PROB_RMV(:,i), RESMOVE(:,i), &
					RESPARAM(:,:,i), &
					Nvib(i), Nbend(i), Ntor(i), &
					MaxVib, MaxBend, MaxTor, &
					IVIB(:,:,i), IBEND(:,:,i), ITOR(:,:,i), &
					RVIB(:,:,i), RBEND(:,:,i), RTOR(:,:,i), &
					BEADTYPE(:,i), GROUPTYPE(:,i), &
					BONDSAPART(:,:,i), Nham, Nljgrs, &
					EPS, SIG, CP, ALP, RMAX, &
					Niongrs, CHARGE, &
					f14_lj, f14_ion, &
					BoxSize, UResInt(i), &
					UResLj(:,i), UResIon(:,i) )

end do

if( Nres == 0 ) then

	Nrefrmoves = Nequil + Nprod + 1
	Nresmol = 1

end if

allocate( Xresmols( Nresmol, MaxBeads, Nres ) )
allocate( Yresmols( Nresmol, MaxBeads, Nres ) )
allocate( Zresmols( Nresmol, MaxBeads, Nres ) )
allocate( Uint_resm( Nresmol, Nres ) )
allocate( Ulj_resm( Nresmol, Nham, Nres ) )
allocate( Uion_resm( Nresmol, Nham, Nres ) )

do i = 1, Nres

	call ResSetup( reslen(i), &
					Xres(:,i), Yres(:,i), Zres(:,i), &
					Nresmol, MaxBeads, Nrefrmoves, Nham, &
					Xresmols(:,:,i), Yresmols(:,:,i), &
					Zresmols(:,:,i), Uint_resm(:,i), &
					Ulj_resm(:,:,i), Uion_resm(:,:,i), &
					Nresmv(i), MaxResmv, MaxResp, &
					PROB_RMV(:,i), RESMOVE(:,i), &
					RESPARAM(:,:,i), &
					Nvib(i), Nbend(i), Ntor(i), &
					MaxVib, MaxBend, MaxTor, &
					IVIB(:,:,i), IBEND(:,:,i), ITOR(:,:,i), &
					RVIB(:,:,i), RBEND(:,:,i), RTOR(:,:,i), &