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), &
					RES_ATT, RES_SUC, &
					BEADTYPE(:,i), GROUPTYPE(:,i), &
					BONDSAPART(:,:,i), Nljgrs, &
					EPS, SIG, CP, ALP, RMAX, &
					Niongrs, CHARGE, &
					f14_lj, f14_ion, &
					BETA, LNW, BoxSize, &
					UResInt(i), UResLj(:,i), UResIon(:,i), &
					Seed, th_h, phi_h )

end do

if( Nres > 0 ) then

	write(*,*) RES_ATT(1), RES_SUC(1)
	write(*,*) RES_ATT(2), RES_SUC(2)

	write(*,*) ' Uintra =   ', sum(Uint_resm) / real(Nresmol)
	write(*,*) ' Ulj =   ', sum(Ulj_resm) / real(Nresmol)
	write(*,*) ' Uion =   ', sum(Uion_resm) / real(Nresmol)
	write(*,*) ' Time =   ', timef() - tstart, '  sec.'			! # C_ALPHA
	!write(*,*) ' Time =   ', etime(xxx) - tstart, '  sec.'		! # U_ALPHA

end if

open( 50, file = PlotsFile )	
close(50, status = 'delete' )

open( 90, file = WFile )	
close(90, status = 'delete' )

if( histtype == 'list' ) then

	UN_HistFile = trim(UN_HistFile)//'.dat'

	do i = 1, MaxSp

		N_HistFile(i) = trim(N_HistFile(0))//'_'//char(48+i)//'.dat'

	end do

	allocate( Npart( Nstorage/Ncollect + 1, MaxSp ) )
	allocate( Energy( Nstorage/Ncollect + 1 ) )

else

	N_HistFile(1) = trim(N_HistFile(0))//'.dat'

end if

if( histtype == 'matrix' ) then

	inquire (file = N_HistFile(1), exist = exists)

	if( .NOT. exists ) then

		Nsim_min = MaxMol
		Nsim_max = 0

	end if

	if( FromDisk .AND. Nequil > 0 ) then

		Urange = Umax - Umin
		Urange = Urange / ( 1.0 + 2.0 * SafetyFactor )

		Usim_min = Umin	+ SafetyFactor * Urange
		Usim_max = Umax	- SafetyFactor * Urange

		Ubins = 0

	else if( FromDisk .AND. Nequil == 0 .AND. exists ) then

		Usim_min = 1.0e10
		Usim_max = -1.0e10

		allocate( UN_HIST( Ubins, 0:MaxMol, Nham ) )

		UN_HIST = 0.0

		open(80, file = N_HistFile(1))

		do i = Nsim_min(1), Nsim_max(1)

			read(80,*) j, N_HIST(i,1,1:Nham)

		end do

		close(80)

		do h = 1, Nham

			filename = trim(UN_HistFile)//'_h'//char( 48 + h/10 )//char(48 + mod(h,10))//'.dat'

			open( 70, file = filename )
		 
			read(70,*)
			read(70,*)

			do i = Nsim_min(1), Nsim_max(1)

				read(70,*) j, j, estart

				k = nint( ( estart - 0.5*Uwidth - Umin ) / Uwidth ) + 1

				read(70,*) UN_HIST( k:k+j-1, i, h )

			end do

			close(70)

		end do 

	else if( FromDisk .AND. Nequil == 0 .AND. .NOT. exists ) then

		allocate( UN_HIST( Ubins, 0:MaxMol, Nham ) )

		UN_HIST = 0.0

	end if

end if

New = .True.

call banner( ResultsFile, ' ', Infostring, Hostname, Datestring, &
			 Timestring, New )

call WriteData( InitialConf, 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 )

open( 10, file = 'lrcorr.dat' )

do i = 1, 605

	read(10,*) XLRCORR(i), ELRCORR(i)

end do

close(10)

if( histtype == 'list' ) then

	open(100, file=UN_HistFile )

	write(100, "(1x,3(G11.5), 3(F10.4,2x))")	& 
			1.0/BETA(1), log(ZETA(1:Nsp,1)) / BETA(1), &
			BoxSize, BoxSize, BoxSize

	close(100)

	Nsim_max = 0
	Nsim_min = MaxMol

end if

call Fourier_Setup( Kmax, Alpha, CONST, KX, KY, KZ, Nkvec )

allocate( SUMQEXPV( Nkvec, Nham ) )

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

write(*,*)
write(*,"(A)") ' The system has been created. '

NATT = 0
NSUC = 0
TOT_NATT = 0
TOT_NSUC = 0

do h = 1, Nham
		
	UION(h) = USURF(h) + UFOURIER(h) + UREAL(h) - USELF_CH(h) - &
	             sum( USELF_MOL( 1:Nmol(0), h) ) + &
				 sum( UION_MOL( 1:Nmol(0), h) )

	ULJ(h) = ULJBOX(h) + ULJLR(h) + sum( ULJ_MOL( 1:Nmol(0), h) )

	UTOT(h) = UION(h) + ULJ(h)
						
end do

write(*,*)
write(*,"(A,T15,I8)") ' MC Step = ', StartConf 
write(*,"(A,T15,F15.4)") ' Rho = ', &
	dot_product( Nmol(1:Nsp) , MASSsp(1:Nsp) ) * 1.0e27 / Nav / Volume
write(*,11) ' UTotal = ', UTOT(1:Nham)
write(*,11) ' LNPSI = ', LNPSI(1:Nham)
write(*,"(A,T15,4I8)") ' Nmol = ', Nmol(1:Nsp)
write(*,"(A,T15,4I8)") ' Tag Value = ', tag_val(1:Nsp) 
write(*,"(A,T15,4I8)") ' Tag Mol. = ', tag_mol(1:Nsp) 

11 format( A, T15, <Nham>F15.4 )

Istore = 1

nt1 = Nmol(0)

do i = StartConf + 1, StartConf + Nequil + Nprod

	Moves = Moves + 1

	Ranq = ran2(Seed)

	if( Ranq < PROB_MOVE(1) ) then

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

		if( Nmol(0) > 0 ) then
		
			PROB_SP_DR(1) = Nmol(1)

			do k = 2, Nsp

				PROB_SP_DR(k) = PROB_SP_DR(k-1) + Nmol(k)

			end do

			PROB_SP_DR = PROB_SP_DR / PROB_SP_DR(Nsp)

		end if

		call Disp_Rot( MaxSp, Nmol, Nlj, Nion, &
					   Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
					   Xion, Yion, Zion, TYPEion, DAMPion, &
					   BoxSize, DXYZ, DROT, LENGTHlj, LENGTHion, SPECIES, &
					   STARTlj, STARTion, PROB_DR, PROB_SP_DR, &
					   Nham, BETA, LNW, LNPSI, LnPi, &
					   Nljgrs, EPS, SIG, CP, ALP, RMAX, MASSlj, &
					   Niongrs, CHARGE, MASSion, &
					   Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
					   EXPX, EXPY, EXPZ, SUMQEXPV, SUMQX, SUMQY, SUMQZ, &
					   ULJBOX, UFOURIER, UREAL, USURF, &
					   DispOrRot, SpID, Success, Seed )

		NATT( DispOrRot, SpID ) = NATT( DispOrRot, SpID ) + 1

		if( Success ) NSUC( DispOrRot, SpID ) = NSUC( DispOrRot, SpID ) + 1

		MOVETIME( DispOrRot ) = MOVETIME( DispOrRot ) + timef() - tt		! # C_ALPHA
!		MOVETIME( DispOrRot ) = MOVETIME( DispOrRot ) + etime(xxx) - tt		! # U_ALPHA

	else if( Ranq < PROB_MOVE(2) ) then

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

		if( ran2(Seed) < 0.5 ) then

			Ranq2 = ran2(Seed)
			SpID = 0
			j = 1

			do while ( SpID == 0 )

				if( Ranq2 < PROB_SP_CD(j) ) SpID = j
				
				j = j + 1

			end do

			if( tag_val(SpID) == 0 ) then

				call 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(SpID), tag_mol(SpID), &
						  Nres, Nresmol, reslen, Xresmols, &
						  Yresmols, Zresmols, Uint_resm, &
						  Ulj_resm, Uion_resm, &
						  f14_lj, f14_ion, ER_HIST, Seed )

				NATT( 3, SpID ) = NATT( 3, SpID ) + 1

				if( Success ) NSUC( 3, SpID ) = NSUC( 3, SpID ) + 1

				MOVETIME(3) = MOVETIME(3) + timef() - tt		! # C_ALPHA
!				MOVETIME(3) = MOVETIME(3) + etime(xxx) - tt		! # U_ALPHA

			else 

				call 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(SpID), tag_mol(SpID), +1, &
							   f14_lj, f14_ion, &
							   SpID, Success, Seed )

				NATT( 6, SpID ) = NATT( 6, SpID ) + 1

				if( Success ) NSUC( 6, SpID ) = NSUC( 6, SpID ) + 1

				MOVETIME(6) = MOVETIME(6) + timef() - tt		! # C_ALPHA
!				MOVETIME(6) = MOVETIME(6) + etime(xxx) - tt		! # U_ALPHA

			end if
		
		else

			Ranq2 = ran2(Seed)
			SpID = 0
			j = 1

			do while ( SpID == 0 )

				if( Ranq2 < PROB_SP_CD(j) ) SpID = j
				
				j = j + 1

			end do

			if( tag_val(SpID) == 0 ) tag_val(SpID) = EESTEPS(SpID)

			if( tag_val(SpID) == 1 ) then

				call Remove( Nlj, Nion, MaxMol, MaxSp, MaxEE, Nmol, &
							 Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
							 Xion, Yion, Zion, TYPEion, DAMPion, &
							 MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
							 MaxBeads, METHOD, MaxInt, MaxReal, INTPARAM, &
							 REALPARAM, BEADTYPE, BONDSAPART, &
							 ERSTEPS, ERSTART, EREND, EESTEPS, EEW, &
							 BoxSize, SPECIES, LENGTHlj, LENGTHion, STARTlj, &
							 STARTion, 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(SpID), tag_mol(SpID), rmol, &
							 Nres, Nresmol, reslen, Xresmols, &
							 Yresmols, Zresmols, Uint_resm, &
							 Ulj_resm, Uion_resm, &
							 f14_lj, f14_ion, ER_HIST, Seed )

				if( success ) then

					do j = 1, Nsp

						if( tag_mol(j) > rmol ) then
						
							tag_mol(j) = tag_mol(j) - 1

						end if

					end do

				end if

				NATT( 4, SpID ) = NATT( 4, SpID ) + 1

				if( Success ) NSUC( 4, SpID ) = NSUC( 4, SpID ) + 1

				MOVETIME(4) = MOVETIME(4) + timef() - tt			! # C_ALPHA
!				MOVETIME(4) = MOVETIME(4) + etime(xxx) - tt			! # U_ALPHA
			
			else

				call 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(SpID), tag_mol(SpID), -1, &
							   f14_lj, f14_ion, &
							   SpID, Success, Seed )

				NATT( 7, SpID ) = NATT( 7, SpID ) + 1

				if( Success ) NSUC( 7, SpID ) = NSUC( 7, SpID ) + 1

				MOVETIME(7) = MOVETIME(7) + timef() - tt		! # C_ALPHA
!				MOVETIME(7) = MOVETIME(7) + etime(xxx) - tt		! # U_ALPHA

			end if
			
		end if
	
	else

		tt = timef()							! # C_ALPHA
!		tt = etime(xxx)							! # U_ALPHA
		
		call ReGrow( Nlj, Nion, MaxMol, MaxSp, Nmol, &
					 Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
					 Xion, Yion, Zion, TYPEion, DAMPion, &
				     NSTEPS, MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
				     MaxBeads, METHOD, MaxInt, MaxReal, INTPARAM, &
				     REALPARAM, BEADTYPE, BONDSAPART, &
					 BoxSize, SPECIES, LENGTHlj, LENGTHion, &
					 STARTlj, STARTion, PROB_SP_RG, &
				     Nham, BETA, 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, &
					 Nres, Nresmol, reslen, Xresmols, &
					 Yresmols, Zresmols, Uint_resm, &
					 Ulj_resm, Uion_resm, &
					 f14_lj, f14_ion, Seed )

		NATT( 5, SpID ) = NATT( 5, SpID ) + 1

		if( Success ) NSUC( 5, SpID ) = NSUC( 5, SpID ) + 1

		MOVETIME(5) = MOVETIME(5) + timef() - tt		! # C_ALPHA
!		MOVETIME(5) = MOVETIME(5) + etime(xxx) - tt		! # U_ALPHA

	end if

	do h = 1, Nham
		
		UION(h) = USURF(h) + UFOURIER(h) + UREAL(h) - USELF_CH(h) - &
		             sum( USELF_MOL( 1:Nmol(0), h) ) + &
					 sum( UION_MOL( 1:Nmol(0), h) )

		ULJ(h) = ULJBOX(h) + ULJLR(h) + &
					sum( ULJ_MOL( 1:Nmol(0), h) )

		UTOT(h) = UION(h) + ULJ(h)
						
	end do

	do j = 1, Nsp
	
		if( tag_val(j) == 0 ) then

			if( Nmol(j) > nt1(j) ) nin(j) = nin(j) + 1
			if( Nmol(j) < nt1(j) ) nout(j) = nout(j) + 1

			nt1(j) = Nmol(j)

		end if

	end do

	! for EE weights

	do k = 1, Nsp

		ALP_PROB(tag_val(k),k) = ALP_PROB(tag_val(k),k) + 1

	end do

	! for HS stuff

	PSI_PI = exp( LNPSI - LnPi )

	AVG_PSI_PI = AVG_PSI_PI + PSI_PI

	if( i > StartConf + sum( NST_TOT(1:stage) ) - NST_TOT(stage) + NST_WT(stage) ) then	  ! For calculation of W

		U_AVG = U_AVG + UTOT
		
		PSI_PI_AVG = PSI_PI_AVG + PSI_PI

		WCounts = WCounts + 1

	end if
	
	if( .NOT. FromDisk .AND. i > int( 0.25 * Nequil ) &
		.AND. i <= Nequil ) then

		do j = 1, Nsp

			Nsim_min(j) = min( Nmol(j), Nsim_min(j) )
			Nsim_max(j) = max( Nmol(j), Nsim_max(j) )

		end do

		Usim_min = min( minval(UTOT), Usim_min )
		Usim_max = max( maxval(UTOT), Usim_max )

	else if( FromDisk .AND. i <= StartConf + Nequil ) then

		do j = 1, Nsp

			Nsim_min(j) = min( Nmol(j), Nsim_min(j) )
			Nsim_max(j) = max( Nmol(j), Nsim_max(j) )

		end do

		Usim_min = min( minval(UTOT), Usim_min )
		Usim_max = max( maxval(UTOT), Usim_max )

	end if

	if( mod( i, NinCycle ) == 0 ) then

		write(*,*)
		write(*,"(A,T15,I8)") ' MC Step = ', i 
		write(*,"(A,T15,F15.4)") ' Rho = ', &
			dot_product( Nmol(1:Nsp) , MASSsp(1:Nsp) ) * 1.0e27 / Nav / Volume

		write(*,11) ' UTotal = ', UTOT(1:Nham)
		write(*,11) ' LNPSI = ', LNPSI(1:Nham)

		! format 11 is format( A, T15, <Nham>F15.4 )

		write(*,"(A,T15,4I8)") ' Nmol = ', Nmol(1:Nsp)
		write(*,"(A,T15,4I8)") ' Tag Value = ', tag_val(1:Nsp) 
		write(*,"(A,T15,4I8)") ' Tag Mol. = ', tag_mol(1:Nsp) 
		write(*,"(A,T15,4I8)") ' n in = ', nin(1:Nsp) 
		write(*,"(A,T15,4I8)") ' n out = ', nout(1:Nsp)
		write(*,"(A,T15,10I8)") ' probs = ', ALP_PROB(0:EESTEPS(1)-1,1:Nsp)

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

		if ( tot_time < 60.0 ) then
	
			write(*, "(1X, A, T15, F15.4, A)")  &
				'Total Time = ', tot_time, ' sec.'

		else if ( tot_time < 3600.0 ) then

			write(*, "(1X, A, T15, F15.4, A)")  &
				'Total Time = ', tot_time/60.0, ' min.'

		else

			write(*, "(1X, A, T15, F15.4, A)")  &
				'Total Time = ', tot_time/3600.0, ' hrs.'

		endif
		
		open( 50, file = PlotsFile, position = 'append' )

		write(50,12) i,  &
			dot_product( Nmol(1:Nsp) , MASSsp(1:Nsp) ) * 1.0e27 / Nav / ( Volume ), &
			UTOT(:), LNPSI(:), Nmol(1:Nsp)

		12 format( I8, <2*Nham+1>(1X,G15.7), <Nsp>I8	)
			
		close(50)	

		Moves = 0

	end if

	if( mod( i, Nrefresh ) == 0 ) then

		do m = 1, Nres

			do j = 1, Nresmol

				Xres(1:reslen(m),m) = Xresmols(j,1:reslen(m),m)
				Yres(1:reslen(m),m) = Yresmols(j,1:reslen(m),m)
				Zres(1:reslen(m),m) = Zresmols(j,1:reslen(m),m)

				do k = 1, Nrefrmoves

					call ResUpdate( reslen(m), Xres(:,m), &
									Yres(:,m), Zres(:,m), &
									Nresmv(m), MaxResmv, MaxResp, &
									PROB_RMV(:,m), RESMOVE(:,m), &
									RESPARAM(:,:,m), &
									Nvib(m), Nbend(m), Ntor(m), &
									MaxVib, MaxBend, MaxTor, &
									IVIB(:,:,m), IBEND(:,:,m), ITOR(:,:,m), &
									RVIB(:,:,m), RBEND(:,:,m), RTOR(:,:,m), &
									RES_ATT, RES_SUC, &
									MaxBeads, &
									BEADTYPE(:,m), GROUPTYPE(:,m), &
									BONDSAPART(:,:,m), Nham, Nljgrs, &
									EPS, SIG, CP, ALP, RMAX, &
									Niongrs, CHARGE, &
									f14_lj, f14_ion, &
									BETA, LNW, BoxSize, &
									UResInt(m), UResLj(:,m), UResIon(:,m), &
									Seed, th_h, phi_h )

				end do

				Xresmols(j,1:reslen(m),m) = Xres(1:reslen(m),m)
				Yresmols(j,1:reslen(m),m) = Yres(1:reslen(m),m)
				Zresmols(j,1:reslen(m),m) = Zres(1:reslen(m),m)

				Uint_resm(j,m) = UResInt(m)
				Ulj_resm(j,:,m) = UResLj(:,m)
				Uion_resm(j,:,m) = UResIon(:,m)

			end do

		end do

	end if

	if( mod( i, Ncollect ) == 0 .AND. i > StartConf + Nequil &
		.AND. sum(tag_val) == 0 ) then

		PSI_PI = PSI_PI	* exp( LNW ) * real( Nham )
		
		do h = 1, Nham

			do j = 1, Nsp

				N_HIST( Nmol(j), j, h ) = N_HIST( Nmol(j), j, h ) + PSI_PI(h)

				Nsim_min(j) = min( Nmol(j), Nsim_min(j) )
				Nsim_max(j) = max( Nmol(j), Nsim_max(j) )

			end do

		end do

		if( histtype == 'list' ) then

			nentry = nentry + 1

			do j = 1, Nsp

				Npart( nentry, j ) = Nmol(j)

			end do

			Energy( nentry ) = UTOT(1)

			Usim_min = min( UTOT(1), Usim_min )
			Usim_max = max( UTOT(1), Usim_max )

		else if( histtype == 'matrix' ) then
		
			do h = 1, Nham

				if( UTOT(h) > Umin .AND. UTOT(h) < Umax ) then
			
					j = int( ( UTOT(h) - Umin ) / Uwidth ) + 1

					UN_HIST( j, Nmol(0), h ) = UN_HIST( j, Nmol(0), h ) + PSI_PI(h)

				else

					if( UTOT(h) < Umin ) then

						write(*,*) ' Umin was exceeded! '
						write(*,*) ' Umin = ', Umin
						write(*,*) ' Energy = ', UTOT(h)

					else if( UTOT(h) > Umax	) then

						write(*,*) ' Umax was exceeded! '
						write(*,*) ' Umax = ', Umax
						write(*,*) ' Energy = ', UTOT(h)

					end if

					write(*,*) ' The program will stop. '
					stop

				end if

				Usim_min = min( UTOT(h), Usim_min )
				Usim_max = max( UTOT(h), Usim_max )

			end do
		
		end if

	end if

	if( mod( i, Nadjust ) == 0 .AND. i <= StartConf + Nequil ) then

		TOT_NATT = TOT_NATT + NATT
		TOT_NSUC = TOT_NSUC + NSUC

		call NewMaxima( Nsp, MaxSp, nmtype, DXYZ, DROT, NATT, NSUC )

	end if

	if( mod( i, Neew ) == 0 .AND. i <= StartConf + Nequil ) then

		do j = 1, Nsp

			call EEWeights( EESTEPS(j), EEW(:,j), ALP_PROB(0:EESTEPS(j),j) )

			write(*,*)
			write(*,"(A,I2)") ' New EE Weights for Species ', j
			write(*,"(10f10.6)") EEW(1:EESTEPS(j),j)
			write(*,*)

		end do

		ALP_PROB = 0

	end if

	if( mod( i, Nstorage ) == 0 ) then

		if( i <= StartConf + Nequil ) then

			call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
						    i, DXYZ, DROT, MaxMol, Ubins, Usim_min, Usim_max, &
							tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
						    LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
						    MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
							Nsim_min(1), Nsim_max(1), Nham, LNW, Seed )

		else

			if( histtype == 'matrix' ) then

				call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
								i, DXYZ, DROT, MaxMol, Ubins, Umin, Umax, &
								tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
								LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
								MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
								Nsim_min(1), Nsim_max(1), Nham, LNW, Seed )

				call WriteHist( MaxMol, MaxSp, Ubins, Nham, BETA, ZETA, Uwidth, Umin, &	
				  				UN_HIST, N_HIST, BoxSize, UN_HistFile, N_HistFile(1), &	
								Nsim_min, Nsim_max )								  
							
			else if( histtype == 'list' ) then
			
				call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
								i, DXYZ, DROT, MaxMol, Ubins, Usim_min, Usim_max, &
								tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
								LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
								MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
								Nsim_min(1), Nsim_max(1), Nham, LNW, Seed )

				open(90, file=UN_HistFile, position='append')
				
				do j = 1, nentry
					
					write(90,14) Npart(j,1:Nsp), Energy(j)

					14 format(1x,<Nsp>I4,G16.7)

				end do

				close(90)

				nentry = 0	

				do j = 1, MaxSp
				
					open(100, file=N_HistFile(j))
					
					do k = Nsim_min(j), Nsim_max(j)

						write(100,"(I8,G16.7E3)") k, N_HIST(k,j,1)

					end do

					close(100)
				
				end do		 
							
			end if												

		end if

		open(19, file = 'er.dat')
		
		do k = 1, Nsp

			do j = 1, ERSTEPS(1) + 1

				write(19,*) j, ER_HIST(j, k, 1:2)

			end do

			write(19,*) 

		end do

		close(19)
	
	end if

	if( i == StartConf + sum( NST_TOT(1:stage) ) ) then

		open(90, file = WFile, position = 'append')

		write(90,*)  
		write(90,"(A,I4)") ' Results after stage ', stage 
		write(90,*)

		if( STOP_H(stage) < 0 ) then

			STOP_H(stage) = - STOP_H(stage)

!			U_AVG = U_AVG / real( WCounts )
			
!			do h = 2, Nham

!				BETA(h) = BETA(1) * U_AVG(1) / U_AVG(h)

!			end do 

			call predict( Nsim_min(1), Nsim_max(1), Ubins, Umin, Uwidth, Nham, &
						  MaxMol, UN_HIST, BETA, ZETA )

			if( STOP_H(stage) /= 100 ) then
			
				PSI_PI_AVG = PSI_PI_AVG / real( WCounts )

				call Weights( Nham, STOP_H(stage), LNW, PSI_PI_AVG )

			end if

			write(90,"(A)") ' Hamiltonian    Weight         Temperature    Activity ' 

			do h = 1, Nham

				write(90,"(5x,I3,T15,G14.7E3,T32,F9.4,T46,F8.4)") h, exp( LNW(h) ), 1.0/BETA(h), log(ZETA(1,h))

			end do

		else if( STOP_H(stage) > 0 ) then

			LNW_old = LNW

			PSI_PI_AVG = PSI_PI_AVG / real( WCounts )

			call Weights( Nham, STOP_H(stage), LNW, PSI_PI_AVG )

			write(90,"(A)") ' Hamiltonian    Weight         PSI_PI_AVG      Coverage ' 

			do h = 1, Nham

				write(90,"(5x,I3,T15,G14.7E3,T32,G14.7E3,T47,F8.5)") &
					h, exp( LNW(h) ), PSI_PI_AVG(h), PSI_PI_AVG(h) * exp(LNW_old(h)) * real(Nham)

			end do

		end if

		close(90)

		PSI_PI_AVG = 0.0

		WCounts = 0

		if( i > StartConf + Nequil ) then
		
			UN_HIST = 0.0
			N_HIST = 0.0

		end if

		Nsim_min = MaxMol
		Nsim_max = -1

		LNPSI(:) = Nmol(0)*log(ZETA(1,:)) + 3.0*Nmol(0)*log(BoxSize) - BETA(:)*UTOT(:)

		do j = 1, Nsp

			do k = 1, Nmol(j)

			LNPSI = LNPSI - log(real(k))

			end do

		end do

		Largest = maxval( LNW + LNPSI )

		LnPi = log( sum( exp( LNW + LNPSI - Largest ) ) ) + Largest

		stage = stage + 1

		if( Stage > Nstages ) stage = Nstages

	end if

	if( i == startConf + Nequil .OR. i == StartConf + Nequil + Nprod ) then

		TOT_NATT = TOT_NATT + NATT
		TOT_NSUC = TOT_NSUC + NSUC

		if( i == StartConf + Nequil ) then

			EquilOrProd = 1
			stc = StartConf + 1

		else

			EquilOrProd = 2
			stc = StartConf + Nequil + 1

		end if

		if( EquilOrProd == 1 ) then

			if( histtype == 'matrix' ) then

				Urange = Usim_max - Usim_min

				Umax = Usim_max + SafetyFactor * Urange

				Umin = Usim_min - SafetyFactor * Urange

				if( Umax > 0 ) Umax = ( real( int( Umax / Uwidth ) ) + 1.5 ) * Uwidth

				if( Umax < 0 ) Umax = ( real( int( Umax / Uwidth ) ) + 0.5 ) * Uwidth

				if( Umin > 0 ) Umin = ( real( int( Umin / Uwidth ) ) - 0.5 ) * Uwidth

				if( Umin < 0 ) Umin = ( real( int( Umin / Uwidth ) ) - 1.5 ) * Uwidth

				Ubins = nint( ( Umax - Umin ) / Uwidth )

				allocate( UN_HIST( Ubins, 0:MaxMol, Nham ) )

				UN_HIST = 0.0

			end if

		end if

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

		call WriteResults( ResultsFile, EquilOrProd, stc, i, &
						   Nsp, MaxSp, nmtype, TOT_NATT, TOT_NSUC, &
						   NAMEsp, MOVETIME, tot_time, &
						   MaxEE, EESTEPS, ALP_PROB, nin, nout, &
						   Usim_min, Usim_max, &
						   Umin, Umax, Ubins, Nsim_min, Nsim_max )

		if( EquilOrProd == 1 .AND. Nprod > 0 ) then

			Nsim_min = MaxMol
			Nsim_max = 0

			Usim_min = Umax
			Usim_max = Umin

		end if
					
		moves = 0

		NATT = 0
		NSUC = 0

		TOT_NATT = 0
		TOT_NSUC = 0

		ALP_PROB = 0

		nin = 0
		nout = 0

		MOVETIME = 0.0

	end if

end do

Istore = 0

i = StartConf + Nequil + Nprod

call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
			    i, DXYZ, DROT, MaxMol, Ubins, Umin, Umax, &
				tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
			    LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
				MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
				Nsim_min, Nsim_max, Nham, LNW, Seed )

write(*,*) 
write(*,11) ' UTotal = ', UTOT(1:Nham)
write(*,11) ' LNPSI = ', LNPSI(1:Nham)

write(*,*) 

close(50)

end program gcmc