writeresults.f90

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

subroutine WriteResults( FileName, EquilOrProd, StartConf, EndConf, &
						 Nsp, MaxSp, nmtype, NATT, NSUC, NAMEsp, &
						 MOVETIME, tot_time, &
						 MaxEE, EESTEPS, ALP_PROB, nin, nout, &
						 Usim_min, Usim_max, &
						 Umin, Umax, Ubins, Nsim_min, Nsim_max)

implicit none

character*30, intent(in)							:: FileName

integer, intent(in)									:: EquilOrProd

integer, intent(in)									:: StartConf
integer, intent(in)									:: EndConf

integer, intent(in)									:: Nsp

integer, intent(in)									:: MaxSp	

integer, intent(in)									:: nmtype	

integer, dimension(nmtype,MaxSp), intent(in)		:: NATT
integer, dimension(nmtype,MaxSp), intent(in)		:: NSUC

character*15, dimension(Nsp), intent(in)			:: NAMEsp

real, dimension(nmtype)								:: MOVETIME

real												:: tot_time

integer												:: MaxEE
integer, dimension(MaxSp)							:: EESTEPS
integer, dimension(0:MaxEE,MaxSp)					:: ALP_PROB
integer, dimension(MaxSp)							:: nin, nout

real												:: Usim_min, Usim_max
real												:: Umin, Umax

integer												:: Ubins

integer, dimension(MaxSp)							:: Nsim_min, Nsim_max

! Local Stuff

integer									:: i, j

integer, parameter						:: unit = 60

real									:: ratio, total


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

write(unit,"(70A)") ('=', j=1, 70)

if( EquilOrProd == 1 ) then

	write(unit,*) ' Equilibration period averages:  Steps ', StartConf, &
					' to ', EndConf

else

	write(unit,*) ' Production period averages:  Steps ', StartConf, &
					' to ', EndConf

end if

write(unit,"(70A)") ('=', j=1, 70)

write(unit,*)

write(unit, "(A, T48, A, F12.4, A)") ' The minimum energy observed ', '= ', Usim_min, ' K'
write(unit, "(A, T48, A, F12.4, A)") ' The maximum energy observed ', '= ', Usim_max, ' K'

write(unit,*)

do j = 1, MaxSp

	write(unit, "(A, I2, A, T48, A, I8)") ' The minimum number of species', j, ' observed ', '= ', Nsim_min(j)
	write(unit, "(A, I2, A, T48, A, I8)") ' The maximum number of species', j, ' observed ', '= ', Nsim_max(j)
	write(unit,*)

end do

write(unit, "(A, T48, A, F12.4, A)") ' Umin was set to ', '= ', Umin, ' K'
write(unit, "(A, T48, A, F12.4, A)") ' Umax was set to ', '= ', Umax, ' K'
write(unit, "(A, T48, A, I8)") ' The number of bins ', '= ', Ubins
write(unit,*)

do i = 1, Nsp

	total = real( sum( ALP_PROB(0:MaxEE, i) ) )
	if( sum( ALP_PROB(0:MaxEE, i) ) == 0 ) total = 1.0

	write(unit, "(A,I2)") '  Expanded Ensemble Probabilities for Species ', i
	write(unit, "(A/A)")  &
					   '   EE Step        Occurrences      % of total ', &
					   ' -----------    ---------------   ------------ '

	do j = 0, EESTEPS(i) - 1

		write(unit, "(5x, I2, T20, I8, T37, F6.2)")  &
			j, ALP_PROB(j,i), real(ALP_PROB(j,i)) / total * 100.0

	end do

	write(unit,*) 
	write(unit, "(A,I8)") '  Number of Complete Creations = ', nin(i)
	write(unit, "(A,I8)") '  Number of Complete Deletions = ', nout(i)

	write(unit,*) 

end do

write(unit, "(A/A/A)") '  Monte Carlo Move Statistics ', &
  '  Move Type   Species        Attempted   Successful     % Acceptance  ', &
  ' ----------- ---------      ----------- ------------   ------------- '

do j = 1, Nsp

	ratio = 0.0

	if( NATT(1,j) /= 0 ) then
	
		ratio = real( NSUC(1,j) ) / real( NATT(1,j) ) * 100.0

		write(unit, "(2X, A, T14, A, T29, I8, T42, I8, T58, F6.2)") &
			'Displace', NAMEsp(j), NATT(1,j), NSUC(1,j), ratio

	end if

end do

do j = 1, Nsp

	ratio = 0.0

	if( NATT(2,j) /= 0 ) then
	
		ratio = real( NSUC(2,j) ) / real( NATT(2,j) ) * 100.0

		write(unit, "(2X, A, T14, A, T29, I8, T42, I8, T58, F6.2)") &
			'Rotate', NAMEsp(j), NATT(2,j), NSUC(2,j), ratio

	end if

end do

do j = 1, Nsp

	ratio = 0.0

	if( NATT(3,j) /= 0 ) then
	
		ratio = real( NSUC(3,j) ) / real( NATT(3,j) ) * 100.0

		write(unit, "(2X, A, T14, A, T29, I8, T42, I8, T58, F6.2)") &
			'Creation', NAMEsp(j), NATT(3,j), NSUC(3,j), ratio

	end if

end do

do j = 1, Nsp

	ratio = 0.0

	if( NATT(4,j) /= 0 ) then
	
		ratio = real( NSUC(4,j) ) / real( NATT(4,j) ) * 100.0

		write(unit, "(2X, A, T14, A, T29, I8, T42, I8, T58, F6.2)") &
			'Deletion', NAMEsp(j), NATT(4,j), NSUC(4,j), ratio

	end if

end do

do j = 1, Nsp

	ratio = 0.0

	if( NATT(6,j) /= 0 ) then
	
		ratio = real( NSUC(6,j) ) / real( NATT(6,j) ) * 100.0

		write(unit, "(2X, A, T14, A, T29, I8, T42, I8, T58, F6.2)") &
			'Alpha Up', NAMEsp(j), NATT(6,j), NSUC(6,j), ratio

	end if

end do

do j = 1, Nsp

	ratio = 0.0

	if( NATT(7,j) /= 0 ) then
	
		ratio = real( NSUC(7,j) ) / real( NATT(7,j) ) * 100.0

		write(unit, "(2X, A, T14, A, T29, I8, T42, I8, T58, F6.2)") &
			'Alpha Dn', NAMEsp(j), NATT(7,j), NSUC(7,j), ratio

	end if

end do

do j = 1, Nsp

	ratio = 0.0

	if( NATT(5,j) /= 0 ) then
	
		ratio = real( NSUC(5,j) ) / real( NATT(5,j) ) * 100.0

		write(unit, "(2X, A, T14, A, T29, I8, T42, I8, T58, F6.2)") &
			'Regrowth', NAMEsp(j), NATT(5,j), NSUC(5,j), ratio

	end if

end do

total = sum( MOVETIME )

write(unit,*)

write(unit, "(A/A/A)") '  Time Distribution', &
					   '  Move Type        Move Time       % of total ', &
					   ' -----------    ---------------   ------------ '

if ( MOVETIME(1) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Displace', MOVETIME(1), ' sec.', MOVETIME(1) / total	* 100.0

else if ( MOVETIME(1) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Displace', MOVETIME(1)/60.0, ' min.', MOVETIME(1) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Displace', MOVETIME(1)/3600.0, ' hrs.', MOVETIME(1) / total	* 100.0

endif

if ( MOVETIME(2) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Rotate', MOVETIME(2), ' sec.', MOVETIME(2) / total	* 100.0

else if ( MOVETIME(2) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Rotate', MOVETIME(2)/60.0, ' min.', MOVETIME(2) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Rotate', MOVETIME(2)/3600.0, ' hrs.', MOVETIME(2) / total	* 100.0

endif

if ( MOVETIME(3) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Creation', MOVETIME(3), ' sec.', MOVETIME(3) / total	* 100.0

else if ( MOVETIME(3) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Creation', MOVETIME(3)/60.0, ' min.', MOVETIME(3) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Creation', MOVETIME(3)/3600.0, ' hrs.', MOVETIME(3) / total	* 100.0

endif

if ( MOVETIME(4) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Deletion', MOVETIME(4), ' sec.', MOVETIME(4) / total	* 100.0

else if ( MOVETIME(4) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Deletion', MOVETIME(4)/60.0, ' min.', MOVETIME(4) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Deletion', MOVETIME(4)/3600.0, ' hrs.', MOVETIME(4) / total	* 100.0

endif

if ( MOVETIME(6) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Alpha Up', MOVETIME(6), ' sec.', MOVETIME(6) / total	* 100.0

else if ( MOVETIME(6) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Alpha Up', MOVETIME(6)/60.0, ' min.', MOVETIME(6) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Alpha Up', MOVETIME(6)/3600.0, ' hrs.', MOVETIME(6) / total	* 100.0

endif

if ( MOVETIME(7) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Alpha Dn', MOVETIME(7), ' sec.', MOVETIME(7) / total	* 100.0

else if ( MOVETIME(6) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Alpha Dn', MOVETIME(7)/60.0, ' min.', MOVETIME(7) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Alpha Dn', MOVETIME(7)/3600.0, ' hrs.', MOVETIME(7) / total	* 100.0

endif

if ( MOVETIME(5) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Regrow', MOVETIME(5), ' sec.', MOVETIME(5) / total	* 100.0

else if ( MOVETIME(5) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Regrow', MOVETIME(5)/60.0, ' min.', MOVETIME(5) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Regrow', MOVETIME(5)/3600.0, ' hrs.', MOVETIME(5) / total	* 100.0

endif

write(unit,*)

if ( tot_time < 60.0 ) then
	
	write(unit, "(1X, A, F8.2, A)")  &
		' Total Time = ', tot_time, ' sec.'

else if ( tot_time < 3600.0 ) then

	write(unit, "(1X, A, F8.2, A)")  &
		' Total Time = ', tot_time/60.0, ' min.'

else

	write(unit, "(1X, A, F8.2, A)")  &
		' Total Time = ', tot_time/3600.0, ' hrs.'

endif

write(unit,*)

close(unit)

return

end subroutine WriteResults