Program Gibbs		  

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.

! Nphases = The number of phases to simulate.
!	Do not set Nphases to anything but 2.
! Ntemp = Number of temperatures.
! Nham = Number of Hamiltonians.
!	NOTE: Significant modifications are required to use this
!			code to simulate multiple Hamiltonians and/or 
!			temperatures during the same run.  For information
!			on how to do so, please contact JRE or AZP.
!			Otherwise, set Ntemp and Nham equal to 1. 

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

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

integer, parameter								:: Nphases = 2
integer, parameter								:: Ntemp = 1
integer, parameter								:: Nham = 1
integer, parameter								:: MaxSp = 1
integer, parameter								:: MaxMol = 300
integer, parameter								:: MaxLJ = 900
integer, parameter								:: MaxIon = 900
integer, parameter								:: MaxSteps = 6
integer, parameter								:: MaxBeads = 6
integer, parameter								:: MaxInt = 7
integer, parameter								:: MaxReal = 8
integer, parameter								:: Kmax	= 5
integer, parameter								:: Maxkvec = 276 

! 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											:: Nljgrs
integer											:: Niongrs
integer											:: Nsp
integer											:: Seed
integer											:: StartConf
integer											:: Nequil
integer											:: Nprod
integer											:: Nblocks
integer											:: NinCycle
integer											:: Nstorage
integer											:: Nadjust
integer											:: Nkvec = Maxkvec
integer											:: cyc
integer											:: Istore
integer											:: Moves
integer											:: h
integer											:: i
integer											:: j
integer											:: k
integer											:: SpID
integer											:: DispOrRot
integer											:: Increase
integer											:: Don
integer											:: EquilOrProd
integer											:: stc

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,0:Nphases)			:: Nmol = 0
integer, dimension(Maxkvec)						:: KX = 0
integer, dimension(Maxkvec)						:: KY = 0
integer, dimension(Maxkvec)						:: KZ = 0
integer, dimension(Nphases)						:: Nlj = 0
integer, dimension(Nphases)						:: Nion = 0
integer, dimension(MaxLJ,Nphases)				:: TYPElj = 0
integer, dimension(MaxIon,Nphases)				:: TYPEion = 0
integer, dimension(MaxMol,Nphases)				:: SPECIES = 0
integer, dimension(MaxMol,Nphases)				:: LENGTHlj = 0
integer, dimension(MaxMol,Nphases)				:: LENGTHion = 0
integer, dimension(MaxMol,Nphases)				:: STARTlj = 0
integer, dimension(MaxMol,Nphases)				:: STARTion = 0
integer, dimension(5,MaxSp,Nphases)				:: NATT = 0
integer, dimension(5,MaxSp,Nphases)				:: NSUC = 0
integer, dimension(5,MaxSp,Nphases)				:: TOT_NATT = 0
integer, dimension(5,MaxSp,Nphases)				:: TOT_NSUC = 0
integer, dimension(Nphases)						:: Nneg = 0
integer, dimension(Nphases)						:: Npos = 0
integer, dimension(Nphases)						:: ZERO_MOL = 0

logical											:: FromDisk
logical											:: Success
logical											:: Constant = .True.
logical											:: New

character*30									:: InputFile
character*30									:: InitialConf
character*30									:: FinalConf
character*30									:: PlotsFile
character*30									:: ResultsFile
character*5										:: Ensemble

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

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

real, parameter									:: kB = 1.380658e-23
real, parameter									:: Nav = 6.02214e23

real											:: Alpha = 5.6
real											:: Pressure = 0.0
real											:: tt
real											:: tstart
real											:: tot_time
real											:: Prob_ph_dr
real											:: Prob_ph_rg
real											:: Ranq
real											:: Convert
real											:: box

!real, dimension(2)								:: xxx			! # U_ALPHA
real, dimension(4)								:: PROB_MOVE
real, dimension(2)								:: PROB_DR
real, dimension(MaxSp)							:: PROB_SP_DR
real, dimension(MaxSp)							:: PROB_SP_TR
real, dimension(MaxSp)							:: PROB_SP_RG
real, dimension(Nphases)						:: BoxSize
real, dimension(Nphases)						:: CDV
real, dimension(MaxSp)							:: MASSsp  = 0.0
real, dimension(MaxReal,MaxBeads,MaxSp)			:: REALPARAM = 0.0
real, dimension(MaxSp,Nphases)					:: DXYZ
real, dimension(MaxSp,Nphases)					:: DROT
real, dimension(605)							:: XLRCORR
real, dimension(605)							:: ELRCORR
real, dimension(Maxkvec)						:: CONST
real, dimension(MaxLJ,Nphases)					:: Xlj
real, dimension(MaxLJ,Nphases)					:: Ylj
real, dimension(MaxLJ,Nphases)					:: Zlj
real, dimension(MaxIon,Nphases)					:: Xion
real, dimension(MaxIon,Nphases)					:: Yion
real, dimension(MaxIon,Nphases)					:: Zion
real, dimension(MaxMol,Nphases)					:: UINTRA
real, dimension(5)								:: MOVETIME = 0.0
real, dimension(Nphases)						:: AVG_RHO = 0.0
real, dimension(Nphases)						:: AVG_UTOT = 0.0
real, dimension(Nphases)						:: AVG_P = 0.0
real, dimension(Nphases)						:: AVG_Nmol = 0.0
real, dimension(MaxSp,Nphases)					:: AVG_MOLE_FRAC = 0.0
real, dimension(Nphases)						:: STD_RHO = 0.0
real, dimension(Nphases)						:: STD_UTOT = 0.0
real, dimension(Nphases)						:: STD_P = 0.0
real, dimension(Nphases)						:: STD_Nmol = 0.0
real, dimension(MaxSp,Nphases)					:: STD_MOLE_FRAC = 0.0
real, dimension(Nphases)						:: ERR_RHO = 0.0
real, dimension(Nphases)						:: ERR_UTOT = 0.0
real, dimension(Nphases)						:: ERR_P = 0.0
real, dimension(Nphases)						:: ERR_Nmol = 0.0
real, dimension(MaxSp,Nphases)					:: ERR_MOLE_FRAC = 0.0

real, allocatable, dimension(:,:)				:: BETA
real, allocatable, dimension(:,:)				:: LNWSTATE
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(:,:)				:: UINT
real, allocatable, dimension(:,:)				:: ULJ
real, allocatable, dimension(:,:)				:: UTOT
real, allocatable, dimension(:,:)				:: PRESS_PLUS
real, allocatable, dimension(:,:)				:: PRESS_MINUS
real, allocatable, dimension(:,:,:,:)			:: AVG_PRESS

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(:,:)				:: CAVG_RHO
real, allocatable, dimension(:,:)				:: CAVG_VOL
real, allocatable, dimension(:,:,:)				:: CAVG_UINT
real, allocatable, dimension(:,:,:)				:: CAVG_ULJ
real, allocatable, dimension(:,:,:)				:: CAVG_UION
real, allocatable, dimension(:,:,:)				:: CAVG_UTOT
real, allocatable, dimension(:,:,:)				:: CAVG_Nmol
real, allocatable, dimension(:,:,:)				:: CAVG_MOLE_FRAC
real, allocatable, dimension(:,:)				:: CAVG_P
real, allocatable, dimension(:,:)				:: BAVG_RHO
real, allocatable, dimension(:,:)				:: BAVG_UTOT
real, allocatable, dimension(:,:)				:: BAVG_P
real, allocatable, dimension(:,:)				:: BAVG_Nmol
real, allocatable, dimension(:,:,:)				:: BAVG_MOLE_FRAC

real, external									:: ran2
!real, external									:: etime	! # U_ALPHA

complex, dimension(0:Kmax,MaxIon,Nphases)		:: EXPX
complex, dimension(-Kmax:Kmax,MaxIon,Nphases)	:: EXPY
complex, dimension(-Kmax:Kmax,MaxIon,Nphases)	:: EXPZ

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


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

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

read(*,*) InputFile
read(*,*) InitialConf
read(*,*) FinalConf
read(*,*) PlotsFile
read(*,*) ResultsFile

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

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

call LittleRead( InputFile, Ensemble, Nljgrs, Niongrs, NSp )

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(Ntemp,1) )
allocate( LNWSTATE(Ntemp,Nham) )
allocate( SUMQX(Nham,Nphases) )
allocate( SUMQY(Nham,Nphases) )
allocate( SUMQZ(Nham,Nphases) )
allocate( ULJBOX(Nham,Nphases) )
allocate( ULJLR(Nham,Nphases) )
allocate( ULJ_MOL(MaxMol,Nham,Nphases) )
allocate( UREAL(Nham,Nphases) )
allocate( USURF(Nham,Nphases) )
allocate( UFOURIER(Nham,Nphases) )
allocate( USELF_CH(Nham,Nphases) )
allocate( USELF_MOL(MaxMol,Nham,Nphases) )
allocate( UION_MOL(MaxMol,Nham,Nphases) )
allocate( UION(Nham,Nphases) )
allocate( UINT(Nham,Nphases) )
allocate( ULJ(Nham,Nphases) )
allocate( UTOT(Nham,Nphases) )
allocate( PRESS_PLUS(Ntemp,Nham) )
allocate( PRESS_MINUS(Ntemp,Nham) )
allocate( AVG_PRESS(Ntemp,Nham,Nphases,2) )

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
UINT	= 0.0
ULJ = 0.0
UTOT	= 0.0
PRESS_PLUS = 0.0
PRESS_MINUS = 0.0
AVG_PRESS = 0.0

call BigRead( InputFile, InitialConf, Ntemp, Nham, Nljgrs, Niongrs, Nsp, &
			  MaxSp, MaxSteps, MaxBeads, MaxInt, MaxReal, &
			  BETA, Pressure, NAMElj, MASSlj, EPS, SIG, CP, ALP, RMAX, KAPPA, &
			  LAMDA, GAMMA, NAMEion, MASSion, CHARGE, NAMEsp, NSTEPS, &
			  LENLJ, LENION, BEADTYPE, GROUPTYPE, MASSsp, NTRIALS, STEPSTART, &
			  STEPLENGTH, METHOD, INTPARAM , REALPARAM , BONDSAPART, &
			  FromDisk, Nmol, StartConf, BoxSize, PROB_MOVE, &
			  PROB_DR, PROB_SP_TR, PROB_SP_RG, Nequil, Nprod, Nblocks, &
			  NinCycle, Nstorage, Nadjust, DXYZ, DROT, CDV, Seed )

open( 60, file = ResultsFile )

New = .True.

call banner( 60, ' ', Infostring, Hostname, Datestring, Timestring, New )

call WriteData( Ensemble, InitialConf, Ntemp, Nham, Nljgrs, Niongrs, &
				Nsp, MaxSp, MaxSteps, MaxBeads, MaxInt, MaxReal, &
				BETA, Pressure, NAMElj, MASSlj, EPS, SIG, CP, ALP, &
				KAPPA, LAMDA, GAMMA, NAMEion, MASSion, CHARGE, NAMEsp, &
				NSTEPS, LENLJ, LENION, BEADTYPE, GROUPTYPE, MASSsp, NTRIALS, &
				STEPSTART, STEPLENGTH, METHOD, INTPARAM, REALPARAM, &
				FromDisk, Nmol, BoxSize, PROB_MOVE, PROB_DR, &
				PROB_SP_TR, PROB_SP_RG, NinCycle, Nstorage, Nadjust, &
				Seed, 60 )

close(60)

cyc = max( int( Nequil / NinCycle ), int( Nprod / NinCycle ) ) + 1

allocate( CAVG_RHO( cyc, Nphases ) )
allocate( CAVG_VOL( cyc, Nphases ) )
allocate( CAVG_UINT( cyc, Nham, Nphases ) )
allocate( CAVG_ULJ( cyc, Nham, Nphases ) )
allocate( CAVG_UION( cyc, Nham, Nphases ) )
allocate( CAVG_UTOT( cyc, Nham, Nphases ) )
allocate( CAVG_Nmol( cyc, 0:Nsp, 0:Nphases ) )
allocate( CAVG_MOLE_FRAC( cyc, Nsp, Nphases ) )
allocate( CAVG_P( cyc, Nphases ) )

allocate( BAVG_RHO( Nblocks, Nphases ) )
allocate( BAVG_UTOT( Nblocks, Nphases ) )
allocate( BAVG_P( Nblocks, Nphases ) )
allocate( BAVG_Nmol( Nblocks, Nphases ) )
allocate( BAVG_MOLE_FRAC( Nblocks, Nsp, Nphases ) )

CAVG_RHO = 0.0
CAVG_VOL = 0.0
CAVG_UINT = 0.0
CAVG_ULJ = 0.0
CAVG_UION = 0.0
CAVG_UTOT = 0.0
CAVG_Nmol = 0.0
CAVG_MOLE_FRAC = 0.0
CAVG_P = 0.0

BAVG_RHO = 0.0
BAVG_UTOT = 0.0
BAVG_Nmol = 0.0
BAVG_MOLE_FRAC = 0.0
BAVG_P = 0.0

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

do i = 1, 605

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

end do

close(10)

LNWSTATE = -log( real( Ntemp * Nham ) ) 

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

allocate( SUMQEXPV( Nkvec, Nham, Nphases ) )

call Create( InitialConf, Ntemp, Nham, Nljgrs, Niongrs, Nsp, &
			 MaxSp, MaxSteps, MaxBeads, MaxInt, MaxReal, &
			 BETA, EPS, SIG, CP, ALP, RMAX, CHARGE, NSTEPS, &
			 LENLJ, LENION, BEADTYPE, GROUPTYPE, NTRIALS, &
			 STEPSTART, STEPLENGTH, METHOD, INTPARAM, &
			 REALPARAM, BONDSAPART, FromDisk, BoxSize,  &
			 Nlj, Nion, MaxMol, MaxLJ, MaxIon, Nmol, &
			 Xlj, Ylj, Zlj, TYPElj, Xion, Yion, Zion, TYPEion, &
			 SPECIES, LENGTHlj, LENGTHion, STARTlj, STARTion,	&
			 LNWSTATE, 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, &
			 f14_lj, f14_ion, Seed )

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

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

do j = 1, Nphases

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

		ULJ(h,j) = ULJBOX(h,j) + ULJLR(h,j)
		
		UINT(h,j) = sum( UINTRA( 1:Nmol(0,j), j) ) + &
					sum( ULJ_MOL( 1:Nmol(0,j), h, j) ) + &
					sum( UION_MOL( 1:Nmol(0,j), h, j) )
					
		UTOT(h,j) = UION(h,j) + ULJ(h,j)
						
		if( Ensemble == 'NPT' ) UTOT(h,j) = UTOT(h,j) + &
			Pressure * ( BoxSize(j) ** 3 ) * 1.0e-25 / kB
	end do

end do

write(*,*)
write(*,"(A,T15,I8)") ' MC Step = ', StartConf

if( BoxSize(1) == 0.0 ) then

	box = 1.0
	
else 

	box = BoxSize(1)
	
end if
 
write(*,"(A,T15,2F15.4)") ' Rho1 = ', dot_product( Nmol(1:Nsp,1) , MASSsp(1:Nsp) ) * &
						  1.0e27 / Nav / ( box ** 3 )

if( BoxSize(2) == 0.0 ) then

	box = 1.0
	
else 

	box = BoxSize(2)
	
end if
 
write(*,"(A,T15,2F15.4)") ' Rho2 = ', dot_product( Nmol(1:Nsp,2) , MASSsp(1:Nsp) ) * &
						  1.0e27 / Nav / ( box ** 3 )

write(*,"(A,T15,2F15.4)") ' U_Inter = ', UTOT(1:Nham,1:Nphases)
write(*,"(A,T15,2F15.4)") ' U_Intra = ', UINT(1:Nham,1:Nphases)

write(*,"(A,T15,2F15.4)") ' U_lj_mol = ', ULJ_MOL(1,1,:)
write(*,"(A,T15,2F15.4)") ' U_Intra = ', UINTRA(1,:)

cyc = 1

Istore = 1

Prob_ph_dr = real( Nmol(0,1) ) / real( Nmol(0,0) )
Prob_ph_rg = real( Nmol(0,1) ) / real( Nmol(0,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( ran2(Seed) < Prob_ph_dr ) then

			j = 1

		else

			j = 2

		end if

		PROB_SP_DR(1) = Nmol(1,j)

		do k = 2, Nsp

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

		end do

		PROB_SP_DR = PROB_SP_DR / PROB_SP_DR(Nsp)

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

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

		if( Success ) NSUC( DispOrRot, SpID, j ) = NSUC( DispOrRot, SpID, j ) + 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( Ensemble == 'NVT' ) then

			call Volume_NVT( Nmol(0,1:2), Nlj, Nion, Maxmol, Maxlj, Maxion, &
							 Xlj, Ylj, Zlj, TYPElj, &
							 Xion, Yion, Zion, TYPEion, BoxSize, &
							 Constant, CDV(1), &
							 LENGTHlj, LENGTHion, STARTlj, STARTion, &
							 NTemp, BETA, Nham, LNWSTATE, XLRCORR, ELRCORR, &
							 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, ULJLR, &
							 UREAL, UFOURIER, USURF, USELF_CH, USELF_MOL, &
							 Increase, Success, PRESS_PLUS, PRESS_MINUS, Seed )

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

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

			if( i > StartConf + Nequil ) then

				if( Increase == 1 ) then

					AVG_PRESS(:,:,1,2) = AVG_PRESS(:,:,1,2) + PRESS_PLUS(:,:)
					AVG_PRESS(:,:,2,1) = AVG_PRESS(:,:,2,1) + PRESS_MINUS(:,:)

					Npos(1) = Npos(1) + 1
					Nneg(2) = Nneg(2) + 1

				else

					AVG_PRESS(:,:,2,2) = AVG_PRESS(:,:,2,2) + PRESS_PLUS(:,:)
					AVG_PRESS(:,:,1,1) = AVG_PRESS(:,:,1,1) + PRESS_MINUS(:,:)

					Npos(2) = Npos(2) + 1
					Nneg(1) = Nneg(1) + 1

				end if

			end if

		else

			if( BoxSize(2) == 0.0 ) then

				j = 1

			else
		
				j = int( 2.0 * ran2(Seed) ) + 1

			end if

			call Volume_NPT( Nmol(0,j), Nlj(j), Nion(j),  MaxMol, &
					  		 Xlj(:,j), Ylj(:,j), Zlj(:,j), TYPElj(:,j), &
						     Xion(:,j), Yion(:,j), Zion(:,j), TYPEion(:,j), & 
							 BoxSize(j), Constant, CDV(j), Pressure, &
							 LENGTHlj(:,j), LENGTHion(:,j), &
							 STARTlj(:,j), STARTion(:,j), &
						     NTemp, BETA, Nham, LNWSTATE, XLRCORR, ELRCORR,&
						     Nljgrs, EPS, SIG, CP, ALP, RMAX, MASSlj, Niongrs, &
							 CHARGE, MASSion, Alpha, Kmax, Nkvec, KX, KY, KZ, &
							 CONST, EXPX(:,:,j), EXPY(:,:,j), EXPZ(:,:,j), &
							 SUMQEXPV(:,:,j), SUMQX(:,j), SUMQY(:,j), SUMQZ(:,j), &
						     ULJBOX(:,j), ULJLR(:,j), UREAL(:,j), UFOURIER(:,j), &
							 USURF(:,j), USELF_CH(:,j), USELF_MOL(:,:,j), &
							 Increase, Success, PRESS_PLUS, Seed )

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

			if( Success ) NSUC( 3, 1, j ) = NSUC( 3, 1, j ) + 1
			
		end if

		MOVETIME(3) = MOVETIME(3) + timef() - tt		! # C_ALPHA
!		MOVETIME(3) = MOVETIME(3) + etime(xxx) - tt		! # U_ALPHA
	
	else if( Ranq < PROB_MOVE(3) ) then

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

		call Trans( Nlj, Nion, MaxMol, MaxSp, MaxLJ, MaxIon, Nmol, &
			  	    Xlj, Ylj, Zlj, TYPElj, Xion, Yion, Zion, TYPEion, &
					NSTEPS, MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
					MaxBeads, METHOD, MaxInt, MaxReal, INTPARAM, &
					REALPARAM, BEADTYPE, BONDSAPART, BoxSize, SPECIES, &
					LENGTHlj, LENGTHion, STARTlj, STARTion, PROB_SP_TR, &
					NTemp, BETA, Nham, LNWSTATE, 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, Don, &
					Success, f14_lj, f14_ion, Seed )

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

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

		Prob_ph_dr = real( Nmol(0,1) ) / real( Nmol(0,0) )
		Prob_ph_rg = real( Nmol(0,1) ) / real( Nmol(0,0) )

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

	else 

		tt = timef()						! # C_ALPHA
!		tt = etime(xxx)						! # U_ALPHA
		
		if( ran2(Seed) < Prob_ph_rg ) then

			j = 1

		else

			j = 2

		end if

		call ReGrow( Nlj(j), Nion(j), MaxMol, MaxSp, Nmol(:,j), &
				     Xlj(:,j), Ylj(:,j), Zlj(:,j), TYPElj(:,j), &
					 Xion(:,j), Yion(:,j), Zion(:,j), TYPEion(:,j), &
				     NSTEPS, MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
				     MaxBeads, METHOD, MaxInt, MaxReal, INTPARAM, &
				     REALPARAM, BEADTYPE, BONDSAPART, BoxSize(j), &
					 SPECIES(:,j), LENGTHlj(:,j), LENGTHion(:,j), &
					 STARTlj(:,j), STARTion(:,j), PROB_SP_RG, &
				     NTemp, BETA, Nham, LNWSTATE, XLRCORR, ELRCORR, &
				     Nljgrs, EPS, SIG, CP, ALP, RMAX, Niongrs, CHARGE, &
				     Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
				     EXPX(:,:,j), EXPY(:,:,j), EXPZ(:,:,j), &
					 SUMQEXPV(:,:,j), SUMQX(:,j), SUMQY(:,j), SUMQZ(:,j), &
				     ULJBOX(:,j), ULJLR(:,j), UREAL(:,j), UFOURIER(:,j), &
					 USURF(:,j), USELF_CH(:,j), USELF_MOL(:,:,j), &
					 ULJ_MOL(:,:,j), UION_MOL(:,:,j), UINTRA(:,j), &
					 SpID, Success, f14_lj, f14_ion, Seed )


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

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

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

	end if

	do j = 1, Nphases

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

			ULJ(h,j) = ULJBOX(h,j) + ULJLR(h,j) 

			UINT(h,j) = sum( UINTRA( 1:Nmol(0,j), j) ) + &
						sum( ULJ_MOL( 1:Nmol(0,j), h, j) ) + &
						sum( UION_MOL( 1:Nmol(0,j), h, j) )
			
			UTOT(h,j) = UION(h,j) + ULJ(h,j)
						
			if( Ensemble == 'NPT' ) UTOT(h,j) = UTOT(h,j) + &
				Pressure * ( BoxSize(j) ** 3 ) * 1.0e-25 / kB
				 
		end do

	end do

	do j = 1, Nphases

		box = BoxSize(j)

		if( box > -1.0e-7 .AND. box < 1.0e-7 ) box = 1.0

		CAVG_RHO(cyc,j) = CAVG_RHO(cyc,j) + &
						  dot_product( Nmol(1:Nsp,j) , MASSsp(1:Nsp) ) * &
						  1.0e27 / Nav / ( box ** 3 )

		CAVG_VOL(cyc,j) = CAVG_VOL(cyc,j) + BoxSize(j) ** 3

		if( Nmol(0,j) /= 0 ) then

			Convert = kB * Nav / dot_product( Nmol(1:Nsp,j) , MASSsp(1:Nsp) )

		else

			Convert = 0.0

		end if

		CAVG_UINT(cyc,:,j) = CAVG_UINT(cyc,:,j) + UINT(:,j) * Convert

		CAVG_ULJ(cyc,:,j) = CAVG_ULJ(cyc,:,j) + ULJ(:,j) * Convert

		CAVG_UION(cyc,:,j) = CAVG_UION(cyc,:,j) + UION(:,j)	* Convert

		CAVG_UTOT(cyc,:,j) = CAVG_UTOT(cyc,:,j) + UTOT(:,j)	* Convert

		CAVG_Nmol(cyc,:,j) = CAVG_Nmol(cyc,:,j) + Nmol(:,j)

		if( Nmol(0,j) /= 0 ) then
		
			CAVG_MOLE_FRAC(cyc,1:Nsp,j) = CAVG_MOLE_FRAC(cyc,1:Nsp,j) + &
										 real( Nmol(1:Nsp,j) ) / real( Nmol(0,j) )

		else

			ZERO_MOL(j) = ZERO_MOL(j) + 1

		end if

	end do

	CAVG_Nmol(cyc,:,0) = CAVG_Nmol(cyc,:,0) + Nmol(:,0)

	if( mod( i, NinCycle ) == 0 .AND. i - StartConf >= NinCycle ) then

		do j = 1, Nphases

			CAVG_RHO(cyc,j) = CAVG_RHO(cyc,j) / real( Moves )

			CAVG_VOL(cyc,j) = CAVG_VOL(cyc,j) / real( Moves )

			CAVG_UINT(cyc,:,j) = CAVG_UINT(cyc,:,j) / real( Moves )

			CAVG_ULJ(cyc,:,j) = CAVG_ULJ(cyc,:,j) / real( Moves )

			CAVG_UION(cyc,:,j) = CAVG_UION(cyc,:,j) / real( Moves )

			CAVG_UTOT(cyc,:,j) = CAVG_UTOT(cyc,:,j) / real( Moves )
		
			CAVG_Nmol(cyc,:,j) = CAVG_Nmol(cyc,:,j) / real( Moves )

			if( ZERO_MOL(j) == Moves ) then

				CAVG_MOLE_FRAC(cyc,:,j) = 0.0

			else
			
				CAVG_MOLE_FRAC(cyc,:,j) = CAVG_MOLE_FRAC(cyc,:,j) / &
										 real( Moves - ZERO_MOL(j) )

			end if
			
			if( Ensemble == 'NVT' .AND. i > StartConf + Nequil ) then

				if( AVG_PRESS(1,1,j,1) == 0.0 .OR. Nneg(j) == 0 .OR. Npos(j) == 0 ) then

					CAVG_P(cyc,j) = 0.0

				else

					AVG_PRESS(:,:,j,1) = AVG_PRESS(:,:,j,1) / real(Nneg(j))
					AVG_PRESS(:,:,j,2) = AVG_PRESS(:,:,j,2) / real(Npos(j))

					CAVG_P(cyc,j) = kB * 1.0e25 / BETA(1,1) / CDV(1) * &
						log( 0.5 * ( 1.0 / AVG_PRESS(1,1,j,1) + AVG_PRESS(1,1,j,2) ) )

				end if

			end if

		end do

		CAVG_Nmol(cyc,:,0) = CAVG_Nmol(cyc,:,0) / real( Moves )

		write(*,*)
		write(*,"(A,T15,I8)") ' MC Step = ', i 
		write(*,"(A,T15,2F15.4)") ' Rho = ', CAVG_RHO(cyc,1:2)
		write(*,"(A,T15,2F15.4)") ' U_Inter = ', CAVG_UTOT(cyc,1:Nham,1:2)
		write(*,"(A,T15,2F15.4)") ' U_Intra = ', CAVG_UINT(cyc,1:Nham,1:2)

		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' )

		if( Ensemble == 'NVT' ) then

			write(50,"(I8, 30(1X,G15.7))") i, CAVG_RHO(cyc,1:2), CAVG_VOL(cyc,1:2), &
				CAVG_UTOT(cyc,1:Nham,1:2), CAVG_UINT(cyc,1:Nham,1:2), &
				CAVG_Nmol(cyc,1:Nsp,1:2), CAVG_P(cyc,1:2), tot_time/60.0

		else if( Ensemble == 'NPT' ) then

			write(50,"(I8, 30(1X,G15.7))") i, CAVG_RHO(cyc,1:2), CAVG_VOL(cyc,1:2), &
				CAVG_UTOT(cyc,1:Nham,1:2), CAVG_UINT(cyc,1:Nham,1:2), &
				CAVG_Nmol(cyc,1:Nsp,1:2), tot_time/60.0

		end if
			
		close(50)	

		ZERO_MOL = 0
		
		Moves = 0

		cyc = cyc + 1

		AVG_PRESS = 0.0

		Nneg = 0
		Npos = 0

	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, Ensemble, DXYZ, DROT, CDV, &
					    NATT, NSUC, BoxSize )

	end if

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

		call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
					    i, BoxSize, DXYZ, DROT, CDV, MaxMol, &
					    LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
					    MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, Seed )

	end if

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

		do j = 1, Nphases
			
			call statistics( cyc-1, Nblocks, CAVG_RHO(:,j), BAVG_RHO(:,j), &
							  AVG_RHO(j), STD_RHO(j), ERR_RHO(j) ) 

			call statistics( cyc-1, Nblocks, CAVG_UTOT(:,1,j), BAVG_UTOT(:,j), &
							  AVG_UTOT(j), STD_UTOT(j), ERR_UTOT(j) )
							  
			call statistics( cyc-1, Nblocks, CAVG_P(:,j), BAVG_P(:,j), &
							  AVG_P(j), STD_P(j), ERR_P(j) ) 
					  
			call statistics( cyc-1, Nblocks, CAVG_Nmol(:,0,j), BAVG_Nmol(:,j), &
							  AVG_Nmol(j), STD_Nmol(j), ERR_Nmol(j) )
			
			
			do k = 1, Nsp
	
				call statistics( cyc-1, Nblocks, CAVG_MOLE_FRAC(:,k,j), BAVG_MOLE_FRAC(:,k,j), &
								  AVG_MOLE_FRAC(k,j), STD_MOLE_FRAC(k,j), ERR_MOLE_FRAC(k,j) )

			end do

		end do

		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

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

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

		call WriteResults( 60, EquilOrProd, Nblocks, stc, i, &
						   Nsp, BAVG_RHO, BAVG_UTOT, BAVG_NMOL, BAVG_MOLE_FRAC, &
						   BAVG_P, AVG_RHO, AVG_UTOT, AVG_NMOL, &
						   AVG_MOLE_FRAC, AVG_P, ERR_RHO, ERR_UTOT, &
						   ERR_NMOL, ERR_MOLE_FRAC, ERR_P,  MaxSp, &
						   TOT_NATT, TOT_NSUC, NAMEsp, CDV, MOVETIME, &
						   tot_time, Ensemble )

		close(60)

		cyc = 1

		moves = 0

		ZERO_MOL = 0

		CAVG_RHO = 0.0
		CAVG_VOL = 0.0
		CAVG_UINT = 0.0
		CAVG_ULJ = 0.0
		CAVG_UION = 0.0
		CAVG_UTOT = 0.0
		CAVG_Nmol = 0.0
		CAVG_MOLE_FRAC = 0.0

		NATT = 0
		NSUC = 0

		TOT_NATT = 0
		TOT_NSUC = 0

		MOVETIME = 0.0

	end if

end do

Istore = 0

i = StartConf + Nequil + Nprod

call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
			    i, BoxSize, DXYZ, DROT, CDV, MaxMol, &
			    LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
			    MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, Seed )

write(*,*) 
write(*,*) UTOT(1,1:2)
write(*,*) UINT(1,1:2)

write(*,*) ULJ_MOL(1,1,:)
write(*,*) UINTRA(1,:)

end program Gibbs