writedata.f90

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

subroutine WriteData( InputConf, Nham, Nljgrs, Niongrs, &
					  Nsp, MaxSp, MaxSteps, MaxBeads, MaxInt, MaxReal, &
					  BETA, ZETA, NAMElj, MASSlj, EPS, SIG, CP, ALP, KAPPA, &
					  LAMDA, NAMEion, MASSion, CHARGE, NAMEsp, NSTEPS, &
					  LENLJ, LENION, BEADTYPE, GROUPTYPE, NTRIALS, &
					  STEPSTART, STEPLENGTH, METHOD, INTPARAM, REALPARAM, &
					  ERSTEPS, ERSTART, EREND, EESTEPS, MaxEE, BEADDAMP, &
					  FromDisk, Nmol, Volume, nmoves, PROB_MOVE, PROB_DR, &
					  PROB_SP_CD, PROB_SP_RG, NinCycle, Nstorage, Nadjust, &
					  Ncollect, Uwidth, Seed, ResultsFile )

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

! beta is the reciprical temperature.
! mu is the chemical potential term

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

! NAMElj containes the names of the LJ groups.
! MASSlj contains the mass of the LJ groups.
! 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.
! KAPPA contains the correction to the epsilom_ij cross parameters.
! LAMDA contains the correction to the sigma_ij cross parameters.

character*15, dimension(Nljgrs), intent(in)				:: NAMElj
real, dimension(Nljgrs), intent(in)						:: MASSlj
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)			:: KAPPA
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: LAMDA

! NAMEion containes the names of the ionic groups.
! MASSion contains the mass of the ionic groups.
! CHARGE contains the charge of a bead for each hamiltonian.

character*15, dimension(Niongrs), intent(in)			:: NAMEion
real, dimension(Niongrs), intent(in)					:: MASSion
real, dimension(Niongrs,Nham), intent(in)				:: CHARGE

! NAMEsp contains the name of each species.
! 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.

character*15, dimension(Nsp), intent(in)				:: NAMEsp
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 
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)					:: REALPARAM

integer, dimension(MaxSp)								:: ERSTEPS
integer, dimension(MaxSteps,MaxSp)						:: ERSTART
integer, dimension(MaxSteps,MaxSp)						:: EREND
integer, dimension(MaxSp)								:: EESTEPS
integer													:: MaxEE
real, dimension(MaxBeads,MaxEE,MaxSp)					:: BEADDAMP

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

logical, intent(in)										:: FromDisk

! Nmol is the total and per species number of molecules in the simulation.
 
integer, dimension(0:MaxSp), intent(in)					:: Nmol

! BoxSize contains the length of the simulation boxes.

real, intent(in)										:: Volume

! PROB_MOVE is the accumulative probability of performing a 
! Displacement / Rotation, Volume change, transfer, or regrowth.
! PROB_DR is the accumulative probability of performing a Displacement or Rotation.
! PROB_SP_CD is the accumulative probability of selecting a given species for transfer.
! PROB_SP_RG is the accumulative probability of selecting a given species for regrowth.

integer, intent(in)										:: nmoves
real, dimension(nmoves), intent(in)						:: PROB_MOVE 
real, dimension(2), intent(in)							:: PROB_DR 
real, dimension(MaxSp), intent(in)						:: PROB_SP_CD
real, dimension(MaxSp), intent(in)						:: PROB_SP_RG

! NinCycle is the number of MC steps per cycle.
! Nstorage is the number of MC steps before storing the current configuration.
! Nadjust is the number of MC steps before adjusting the maximum displacement, etc.
! Ncollect is the number of MC steps before collecting histogram data.

integer, intent(in)										:: NinCycle
integer, intent(in)										:: Nstorage
integer, intent(in)										:: Nadjust
integer, intent(in)										:: Ncollect

! Uwidth is the width of an energy bin.

real, intent(in)										:: Uwidth

! Seed is the current random number generator seed value.

integer, intent(in)										:: Seed

! ResultsFile is the file name to write the data to

character*30, intent(in)								:: ResultsFile

! Local stuff

integer								:: h, i, j, k, m, n
integer, parameter					:: unit = 60
real								:: tmp
real, parameter						:: Pi = 3.14159265359
real, parameter						:: Nav = 6.02214e23

real, external						:: zbrent3, funk3

open(unit, file = ResultsFile, position = 'append' )

write(unit,"(1x,A)")  'EHSGCMC Simulation'

write(unit,12 ) 'Temperature = ', 1.0/BETA, ' K '
12 format( 1x, A, <Nham>F10.3, A )

do i = 1, Nsp

	write(unit,13) 'ln(Zeta)', i, ' = ', log(ZETA(i,1:Nham))
	13 format( 1x, A, I2, A, <Nham>F10.3 )

end do

write(unit,"(1x,A, F10.3, A)" ) 'Volume = ', Volume, ' A^3 '

write(unit,*)

write(unit,"(1x,A)") 'Exp-6 / LJ Groups: '
write(unit,"(1x,A)") 'Name          Mass   Hamil.  Eps/k_B (K)  Sigma (A)     C      Alpha '

do h = 1, Nham

	do i = 1, Nljgrs

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

			tmp = SIG(i,i,h) * zbrent3( funk3, 0.5, 1.0, 1.0e-7, ALP(i,i,h), CP(i,i,h) )

		else

			tmp = SIG(i,i,h)

		end if
		
		write(unit,"(1x,A,T13,f8.4,T24,I2,T31,f7.3,T44,f7.4,T54,f7.4,T65,f5.2)") NAMElj(i), MASSlj(i), h, &
							    EPS(i,i,h), tmp, CP(i,i,h), ALP(i,i,h)

	end do

end do

do h = 1, Nham

	write(unit, "(1x,A,I2)") 'Cross Parameters for hamiltonian ', h
	
	write(unit, "(T21,A,T51,A)") 'Kappa_ij', 'Lamda_ij'
	write(unit, 10) '  Group i\j', (i, i=1,Nljgrs), (i, i=1,Nljgrs)

10 format( A,T14,<Nljgrs>I8,T44,<Nljgrs>I8 )

	do i = 1, Nljgrs

		write(unit, 20) i, (KAPPA(i,j,h),j=1,Nljgrs) , (LAMDA(i,j,h),j=1,Nljgrs)

20 format( I7,T16,<Nljgrs>F8.4,T46,<Nljgrs>F8.4 )

	end do

end do

if( Niongrs > 0 ) then

	write(unit,*)
	write(unit,"(1x,A)") 'Ionic Groups: '
	write(unit,"(1x,A)") 'Name          Mass      Hamil.       Charge (e) '

	do h = 1, Nham

		do i = 1, Niongrs

			write(unit,"(1x,A,T13,G12.5,T27,I2,T38,SP,F8.5,SS)") NAMEion(i), MASSion(i), h, &
																CHARGE(i,h)

		end do

	end do

end if

do i = 1, Nsp

	write(unit,*)
	write(unit,"(1x,A,A)") trim( NAMEsp(i) ), ' Data'
	write(unit,"(1x,A,I6,I6)") 'Starting number of molecules = ', Nmol(i)
	write(unit,"(1x,A,A)")	'Growth of ', NAMEsp(i)
	write(unit,"(1x,A)") 'Bead  Type  Group  Method   Int. Params    Real Params'

	do j = 1, LENLJ(i) + LENION(i)

		select case ( METHOD(j,i) )

			case( 'Random' )

				write(unit,"(1x,I3,T9,A,T13,I3,T21,A)") &
					j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i)
					
			case( 'Sphere' )
			
				write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,I3,T45,G11.5)") &
					j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), &
					INTPARAM(1,j,i), REALPARAM(1,j,i)

			case( 'Cone' )
			
				write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,3I3,T45,20(G11.5,1x))") &
					j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), &
					INTPARAM(1:3,j,i), REALPARAM(1:2*INTPARAM(3,j,i),j,i) * 180.0/Pi, &
					REALPARAM(2*INTPARAM(3,j,i)+1:3*INTPARAM(3,j,i),j,i)

			case( 'ConeTor' )

				if( INTPARAM(4,j,i) == 1 ) then

					write(unit,"(1x,I3,T9,A,T13,I3,T21,A,T30,4I3,T45,9G11.5)") &
						j, BEADTYPE(j,i), GROUPTYPE(j,i), METHOD(j,i), &