subroutine BigRead( InputFile, InputConf, 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 )

implicit none

! InputFile contains the parameters for the run.
! InputConf has the information for continuing a previous run.

character*30, intent(in)								:: InputFile
character*30, intent(in)								:: InputConf

! Ntemp is the number of temperatures.
! 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)										:: Ntemp
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.
! Pressure is the pressure for a constant pressure run.

real, dimension(Ntemp,1), intent(out)					:: BETA
real, intent(out)										:: Pressure

! NAMElj containes the names of the LJ groups.
! MASSlj contains the mass of the LJ groups.
! EPS is a rank 3 array containing the eps_ij parameters for each hamiltonian.
! SIG is a rank 3 array containing the sigma_ij parameters for each hamiltonian.
! CP is a rank 3 array containing the C_ij parameters for each hamiltonian.
! ALP is a rank 3 array containing the alpha_ij parameters for each hamiltonian.
! RMAX is a rank 3 array containing the Rmax_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(out)			:: NAMElj
real, dimension(Nljgrs), intent(out)					:: MASSlj
real, dimension(Nljgrs, Nljgrs, Nham), intent(out)		:: EPS
real, dimension(Nljgrs, Nljgrs, Nham), intent(out)		:: SIG
real, dimension(Nljgrs, Nljgrs, Nham), intent(out)		:: CP
real, dimension(Nljgrs, Nljgrs, Nham), intent(out)		:: ALP
real, dimension(Nljgrs, Nljgrs, Nham), intent(out)		:: RMAX
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: KAPPA
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: LAMDA
real, dimension(Nljgrs,Nljgrs,Nham), intent(out)		:: GAMMA

! 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(out)			:: NAMEion
real, dimension(Niongrs), intent(out)					:: MASSion
real, dimension(Niongrs,Nham), intent(out)				:: 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.
! MASSsp contains the mass of the species.
! 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.
! BONDAPART indicates how many bonds separate two beads.

character*15, dimension(Nsp), intent(out)				:: NAMEsp
integer, dimension(MaxSp), intent(out)					:: NSTEPS
integer, dimension(MaxSp), intent(out)					:: LENLJ 
integer, dimension(MaxSp), intent(out)					:: LENION
character*5, dimension(MaxBeads,MaxSp), intent(out)		:: BEADTYPE
integer, dimension(MaxBeads,MaxSp), intent(out)			:: GROUPTYPE 
real, dimension(MaxSp), intent(out)						:: MASSsp
integer, dimension(MaxSteps,MaxSp), intent(out)			:: NTRIALS
integer, dimension(MaxSteps,MaxSp), intent(out)			:: STEPSTART 
integer, dimension(MaxSteps,MaxSp), intent(out)			:: STEPLENGTH
character*10, dimension(MaxBeads,MaxSp), intent(out)	:: METHOD
integer, dimension(MaxInt,MaxBeads,MaxSp), intent(out)	:: INTPARAM	
real, dimension(MaxReal,MaxBeads,MaxSp), intent(out)	:: REALPARAM
integer, dimension(MaxBeads,MaxBeads,MaxSp), intent(out):: BONDSAPART

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

logical, intent(out)									:: FromDisk

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

! StartConf is the starting configuration number.

integer, intent(out)									:: StartConf

! BoxSize contains the length of the simulation boxes.

real, dimension(2), intent(out)							:: BoxSize

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

real, dimension(4), intent(out)							:: PROB_MOVE 
real, dimension(2), intent(out)							:: PROB_DR 
real, dimension(MaxSp), intent(out)						:: PROB_SP_TR
real, dimension(MaxSp), intent(out)						:: PROB_SP_RG

! Nequil is the number of equilibration steps.
! Nprod is the number of production steps.
! Nblocks is the number of blocks used for the error analysis.
! 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.

integer, intent(out)									:: Nequil
integer, intent(out)									:: Nprod
integer, intent(out)									:: Nblocks
integer, intent(out)									:: NinCycle
integer, intent(out)									:: Nstorage
integer, intent(out)									:: Nadjust

! DXYZ is the maximum displacement per species.
! DROT is the maximum rotation per species.
! CDV is the constant volume change.

real, dimension(MaxSp,2), intent(out)					:: DXYZ
real, dimension(MaxSp,2), intent(out)					:: DROT
real, dimension(2), intent(out)							:: CDV

! Seed is the current random number generator seed value.

integer, intent(out)									:: Seed


! Local Stuff
integer										:: h, i, j, k, m
integer										:: bst, bend
integer										:: group
integer										:: step
integer										:: Nmol1, Nmol2
integer										:: ref1, ref2
integer										:: len
integer, dimension(MaxBeads)				:: REFBEAD
integer, dimension(MaxBeads)				:: BASE
integer, dimension(MaxBeads)				:: BONDS


character*5									:: btype
character*15								:: Species

real										:: rho1, rho2
real, parameter								:: Pi = 3.14159265359
real, parameter								:: Nav = 6.02214e23

real, external								:: zbrent2, funk2, zbrent3, funk3


open( 20, file = InputFile )

read(20,*)
read(20,*)
read(20,*)
read(20,*) BETA( 1:Ntemp, 1 ) 

BETA = 1.0 / BETA

read(20,*)
read(20,*) Pressure
read(20,*)
read(20,*)
read(20,*)

! Sigma is read in from input file.

do i = 1, Nljgrs

	read(20,*) j, NAMElj(j), MASSlj(j), EPS(j,j,1), SIG(j,j,1), CP(j,j,1), ALP(j,j,1)

end do

read(20,*)
read(20,*)
read(20,*) ( ( KAPPA( i, j, 1 ), j = 1, Nljgrs ), i = 1, Nljgrs )
read(20,*)
read(20,*) ( ( LAMDA( i, j, 1 ), j = 1, Nljgrs ), i = 1, Nljgrs )
read(20,*)
read(20,*) ( ( GAMMA( i, j, 1 ), j = 1, Nljgrs ), i = 1, Nljgrs )

! The following combining rules are used
! eps_ij = sqrt( eps_ii * eps_jj )
! sig_ij = 0.5*( sig_ii + sig_jj )
! alp_ij = sqrt( alp_ii * alp_jj )

do h = 1, Nham

	do i = 1, Nljgrs - 1

		do j = i + 1, Nljgrs

			EPS(i,j,h) = sqrt( EPS(i,i,h) * EPS(j,j,h) ) 
			EPS(j,i,h) = EPS(i,j,h)	* ( 1 - KAPPA(j,i,h) )
			EPS(i,j,h) = EPS(i,j,h)	* ( 1 - KAPPA(i,j,h) )

			ALP(i,j,h) = sqrt( ALP(i,i,h) * ALP(j,j,h) ) 
			ALP(j,i,h) = ALP(i,j,h)	* ( 1 - GAMMA(j,i,h) )
			ALP(i,j,h) = ALP(i,j,h)	* ( 1 - GAMMA(i,j,h) )

			SIG(i,j,h) = 0.5 * ( SIG(i,i,h) + SIG(j,j,h) ) 
			SIG(j,i,h) = SIG(i,j,h)	* ( 1 - LAMDA(j,i,h) )
			SIG(i,j,h) = SIG(i,j,h)	* ( 1 - LAMDA(i,j,h) )

			if( CP(1,1,1) > 0.0 ) then

				CP(i,j,h) = 1.0
				CP(j,i,h) = 1.0
			
				ALP(i,j,h) = sqrt( ALP(i,i,h) * ALP(j,j,h) ) 
				ALP(j,i,h) = ALP(i,j,h)	* ( 1 - KAPPA(j,i,h) )
				ALP(i,j,h) = ALP(i,j,h)	* ( 1 - KAPPA(i,j,h) )

				! Convert from sig_ij to rm_ij

				SIG(i,j,h) = SIG(i,j,h) / zbrent3( funk3, 0.5, 1.0, 1.0e-7, &
													ALP(i,j,h), CP(i,j,h) )

				SIG(j,i,h) = SIG(j,i,h) / zbrent3( funk3, 0.5, 1.0, 1.0e-7, &
													ALP(j,i,h), CP(j,i,h) )

				! find rmax_ij from rm_ij

				RMAX(i,j,h) = SIG(i,j,h) * zbrent2( funk2, 0.0, 0.8, 1.0e-7, &
													ALP(i,j,h), CP(i,j,h) )

				RMAX(j,i,h) = SIG(j,i,h) * zbrent2( funk2, 0.0, 0.8, 1.0e-7, &
													ALP(j,i,h), CP(j,i,h) )

			else

				CP(i,j,h) = -1.0
				CP(j,i,h) = -1.0

			end if

		end do

	end do

end do

if( CP(1,1,1) > 0.0 ) then

	do i = 1, Nljgrs 

		do h = 1, Nham

			! Convert from sig_ii to rm_ii

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

			! find rmax_ii from rm_ii

			RMAX(i,i,h) = SIG(i,i,h) * zbrent2( funk2, 0.0, 0.8, 1.0e-7, &
												ALP(i,i,h), CP(i,i,h) )

			! Water, set O to LJ bead 1.
			!if( i == 1 ) RMAX(i,i,h) = 0.45 * SIG(i,i,h)
			
		end do

	end do

end if

read(20,*)
read(20,*)

if( Niongrs > 0 ) then

	read(20,*)

	do i = 1, Niongrs

		read(20,*) j, NAMEion(j), MASSion(j), CHARGE(j,1)

	end do

end if

read(20,*)
read(20,*)

do i = 1, Nsp

	read(20,*)
	read(20,*) NAMEsp(i), NSTEPS(i), LENLJ(i), LENION(i)
	read(20,*)

	bend = 0

	do while ( bend < LENLJ(i) + LENION(i) )

		read(20,*) bst, bend, btype, group

		do j = bst, bend

			BEADTYPE(j,i) = btype
			GROUPTYPE(j,i) = group

			if( btype == 'LJ' ) then

				MASSsp(i) = MASSsp(i) + MASSlj(group)

			else if( btype == 'ION' ) then

				MASSsp(i) = MASSsp(i) + MASSion(group)

			end if

		end do

	end do

	read(20,*)

	step = 0

	do while ( step < NSTEPS(i) )

		read(20,*) step, NTRIALS(step,i), STEPSTART(step,i), bend

		STEPLENGTH(step,i) = bend - STEPSTART(step,i) + 1

		do k = STEPSTART(step,i), bend

			REFBEAD(k) = k

		end do

	end do

	read(20,*)

	bend = 0

	do while ( bend < LENLJ(i) + LENION(i) )

		read(20,*) bst, METHOD(bst,i)

		bend = bst

		select case ( METHOD(bst,i) )

			case( 'Random' )

			case( 'Sphere' )

				read(20,*) INTPARAM(1,bst,i), REALPARAM(1,bst,i)

			case( 'Cone' )

				read(20,*) INTPARAM(1:3,bst,i), &
						   REALPARAM( 1:3*INTPARAM(3,bst,i), bst, i )

				REALPARAM( 1:2*INTPARAM(3,bst,i), bst, i ) = &
				REALPARAM( 1:2*INTPARAM(3,bst,i), bst, i ) * Pi / 180.0

			case( 'FxBend' ) 

				read(20,*) INTPARAM(1:2,bst,i)
				read(20,*) INTPARAM(3:2+INTPARAM(1,bst,i),bst,i), &
						   REALPARAM(1:1+2*INTPARAM(1,bst,i),bst,i)

				do j = 1, INTPARAM(1,bst,i)

					REALPARAM( 1+2*j, bst, i ) = REALPARAM( 1+2*j, bst, i ) * Pi / 180.0

				end do

			case( 'FxBendTor' )

				read(20,*) INTPARAM(1:4,bst,i)

				if( INTPARAM(3,bst,i) == 1 ) then

					read(20,*) INTPARAM(5:4+INTPARAM(1,bst,i)+2*INTPARAM(2,bst,i),bst,i), &
							   REALPARAM(1:1+2*INTPARAM(1,bst,i)+6*INTPARAM(2,bst,i),bst,i)

				else if( INTPARAM(3,bst,i) == 2 ) then

					read(20,*) INTPARAM(5:4+INTPARAM(1,bst,i)+2*INTPARAM(2,bst,i),bst,i), &
							   REALPARAM(1:1+2*INTPARAM(1,bst,i)+4*INTPARAM(2,bst,i),bst,i)

				end if
				
				do j = 1, INTPARAM(1,bst,i)

					REALPARAM( 1+2*j, bst, i ) = REALPARAM( 1+2*j, bst, i ) * Pi / 180.0

				end do


			case( 'Stretch' )

				read(20,*) INTPARAM(1,bst,i), REALPARAM(1:2,bst,i)

			case( 'StBend' ) 

				read(20,*) INTPARAM(1:2,bst,i)
				read(20,*) INTPARAM(3:2+INTPARAM(1,bst,i),bst,i), &
						   REALPARAM(1:2+2*INTPARAM(1,bst,i),bst,i)

				do j = 1, INTPARAM(1,bst,i)

					REALPARAM( 2+2*j, bst, i ) = REALPARAM( 2+2*j, bst, i ) * Pi / 180.0

				end do

			case( 'StBendTor' )

				read(20,*) INTPARAM(1:4,bst,i)

				if( INTPARAM(3,bst,i) == 1 ) then

					read(20,*) INTPARAM(5:4+INTPARAM(1,bst,i)+2*INTPARAM(2,bst,i),bst,i), &
							   REALPARAM(1:2+2*INTPARAM(1,bst,i)+6*INTPARAM(2,bst,i),bst,i)

				else if( INTPARAM(3,bst,i) == 2 ) then

					read(20,*) INTPARAM(5:4+INTPARAM(1,bst,i)+2*INTPARAM(2,bst,i),bst,i), &
							   REALPARAM(1:2+2*INTPARAM(1,bst,i)+4*INTPARAM(2,bst,i),bst,i)

				end if
				
				do j = 1, INTPARAM(1,bst,i)

					REALPARAM( 2+2*j, bst, i ) = REALPARAM( 2+2*j, bst, i ) * Pi / 180.0

				end do

			case( 'Match' )

				read(20,*) INTPARAM(1,bst,i)

			case( 'Complete' )

		end select

	end do

end do

read(20,*)
read(20,*) FromDisk
read(20,*)

if( FromDisk ) open( 30, file = InputConf )

do i = 1, Nsp

	if( FromDisk ) then

		read(30,*) Species, Nmol1, Nmol2
		read(20,*)

	else 

		read(20,*) Species, Nmol1, Nmol2

	end if
	
	do j = 1, Nsp

		if( Species == NAMEsp(j) ) then

			Nmol(j,1) = Nmol1
			Nmol(j,2) = Nmol2
			Nmol(j,0) = Nmol1 + Nmol2

		end if

	end do

end do

Nmol(0,1) = sum( Nmol(1:Nsp,1) )
Nmol(0,2) = sum( Nmol(1:Nsp,2) )
Nmol(0,0) = Nmol(0,1) + Nmol(0,2)

if( FromDisk ) then

	read(30,*)
	read(30,*) StartConf, Seed

	read(20,*)
	read(20,*)
	read(20,*)

else

	StartConf = 0
	
	read(20,*)
	read(20,*) Seed
	read(20,*)
 
end if


if( FromDisk ) then

	read(30,*) BoxSize(1:2)
	read(20,*) 

else

	read(20,*) rho1, rho2

	BoxSize(1) = dot_product( Nmol(1:Nsp,1) , MASSsp(1:Nsp) ) * &
					1.0e27 / Nav / rho1
	BoxSize(2) = dot_product( Nmol(1:Nsp,2) , MASSsp(1:Nsp) ) * &
					1.0e27 / Nav / rho2

	BoxSize = BoxSize ** (1.0 / 3.0)

end if

read(20,*)
read(20,*) PROB_MOVE(1:4)

do i = 2, 4

	PROB_MOVE(i) = PROB_MOVE(i) + PROB_MOVE(i-1)

end do

PROB_MOVE = PROB_MOVE / PROB_MOVE(4)

read(20,*)
read(20,*) PROB_DR(1:2)

PROB_DR(2) = PROB_DR(1) + PROB_DR(2)

PROB_DR = PROB_DR / PROB_DR(2)

read(20,*)
read(20,*) PROB_SP_TR(1:Nsp)

do i = 2, Nsp

	PROB_SP_TR(i) = PROB_SP_TR(i) + PROB_SP_TR(i-1)

end do

PROB_SP_TR = PROB_SP_TR / PROB_SP_TR(Nsp)

read(20,*)
read(20,*) PROB_SP_RG(1:Nsp)

do i = 2, Nsp

	PROB_SP_RG(i) = PROB_SP_RG(i) + PROB_SP_RG(i-1)

end do

PROB_SP_RG = PROB_SP_RG / PROB_SP_RG(Nsp)

read(20,*)
read(20,*) Nequil, Nprod
read(20,*)
read(20,*) Nblocks
read(20,*) 
read(20,*) NinCycle
read(20,*)
read(20,*) Nstorage
read(20,*)
read(20,*) Nadjust

if( FromDisk ) then

	read(30,*)
	read(30,*)
	read(30,*) DXYZ(1:Nsp,1)
	read(30,*) DROT(1:Nsp,1)
	read(30,*) CDV(1)
	read(30,*)
	read(30,*) DXYZ(1:Nsp,2)
	read(30,*) DROT(1:Nsp,2)
	read(30,*) CDV(2)

else

	DXYZ = 0.1
	DROT = 0.2

	CDV = 0.02 * ( minval( BoxSize )**3 )

	if( BoxSize(2) == 0.0 ) CDV = 0.02 * ( BoxSize(1)**3 ) 

end if

BONDSAPART = 0

do i = 1, Nsp

	BASE(1) = 1
	BONDS(1) = 0

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

		select case ( METHOD(j,i) )

			case( 'Sphere' )

				BASE(j) = INTPARAM(1,j,i)
				BONDS(j) = 1

			case( 'Cone' )

				do k = j, j + INTPARAM(3,j,i) - 1

					BASE(k) = INTPARAM(2,j,i)
					BONDS(k) = 1
				
				end do

			case( 'ConeTor' )

				BASE(j) = INTPARAM(3,j,i)
				BONDS(j) = 1

			case( 'FxBend' )

				BASE(j) = INTPARAM(2,j,i)
				BONDS(j) = 1

			case( 'FxBendTor' )

				BASE(j) = INTPARAM(4,j,i)
				BONDS(j) = 1

			case( 'Stretch' )

				BASE(j) = INTPARAM(1,j,i)
				BONDS(j) = 1

			case( 'StBend' )

				BASE(j) = INTPARAM(2,j,i)
				BONDS(j) = 1

			case( 'StBendTor' )

				BASE(j) = INTPARAM(4,j,i)
				BONDS(j) = 1

			case( 'Match' )

				BASE(j) = INTPARAM(1,j,i)
				BONDS(j) = 0

			case( 'Complete' )

		end select

	end do

	do j = 1, LENLJ(i) + LENION(i) - 1
	
		do k = j + 1, LENLJ(i) + LENION(i)
		
			ref1 = BASE(j)
			ref2 = BASE(k)

			len = BONDS(k)

			m = 0
			
			do while ( m == 0 )

				if( ref1 == ref2 ) then

					len = len + BONDS(j)

					m = 1
			
				else if( ref1 > ref2 ) then

					if(	ref2 /= j .AND. ref1 /= k ) then

						len = len + BONDS(ref1)
						ref1 = BASE(ref1)
						
					else

						m = 1	

					end if
				
				else if( ref2 > ref1 ) then

					if(	ref2 /= j .AND. ref1 /= k ) then

						len = len + BONDS(ref2)
						ref2 = BASE(ref2)

					else

						m = 1

					end if

				end if

			end do
				
			BONDSAPART(j,k,i) = len
			BONDSAPART(k,j,i) = len
			
		end do
		
	end do
	
end do						 

close(20)
close(30)

return

end subroutine BigRead