subroutine BigRead( InputFile, InputConf, 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, Nmin, Nmax, LNW, &
					Nstages, NST_TOT, NST_WT, STOP_H, 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

! Nham is the number of hamiltonians.
! Nljgrs is the number of Lennard-Jones groups in the simulation.
! Niongrs is the number of ionic groups in the simulation.
! Nsp is the number of species in the simulation.

integer, intent(in)										:: Nham
integer, intent(in)										:: Nljgrs
integer, intent(in)										:: Niongrs
integer, intent(in)										:: Nsp

! MaxSp is the maximum number of species in the simulation.
! MaxSteps is the maximum number of CB steps for a molecule.
! MaxBeads is the maximum number of LJ and ionic beads in a molecule.
! MaxInt is the maximum number of integer parameters for a CB method.
! MaxReal is the maximum number of real parameters for a CB method.

integer, intent(in)										:: MaxSp
integer, intent(in)										:: MaxSteps
integer, intent(in)										:: MaxBeads
integer, intent(in)										:: MaxInt
integer, intent(in)										:: MaxReal
integer, intent(in)										:: MaxEE

integer, intent(in)										:: nmoves

! Uwidth is the energy width of the histogram bins.
! beta is the reciprical temperature.
! mu is the chemical potential term

real, intent(out)										:: Uwidth
real, dimension(Nham), intent(out)						:: BETA
real, dimension(MaxSp,Nham), intent(out)				:: ZETA

! NAMElj containes the names of the LJ groups.
! MASSlj contains the mass of the LJ groups.
! EPS contains the eps_ij parameters for each hamiltonian.
! SIG contains the sigma_ij parameters for each hamiltonian.
! CP contains the C_ij parameters for each hamiltonian.
! ALP contains the alpha_ij parameters for each hamiltonian.
! RMAX contains the Rmax 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(MaxBeads,MaxEE,MaxSp), intent(out)		:: BEADDAMP
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

integer, dimension(MaxSp), intent(out)					:: ERSTEPS
integer, dimension(MaxSteps,MaxSp), intent(out)			:: ERSTART
integer, dimension(MaxSteps,MaxSp), intent(out)			:: EREND
integer, dimension(MaxSp), intent(out)					:: EESTEPS
real, dimension(MaxEE,MaxSp), intent(out)				:: EEW

real, dimension(MaxSp), intent(out)						:: MASSsp

! 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), intent(out)				:: Nmol

! StartConf is the starting configuration number.

integer, intent(out)									:: StartConf

! BoxSize contains the length of the simulation boxes.

real, intent(out)										:: Volume

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

real, dimension(nmoves), intent(out)					:: PROB_MOVE 
real, dimension(2), intent(out)							:: PROB_DR 
real, dimension(MaxSp), intent(out)						:: PROB_SP_CD
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 MC steps before histogram data collection.
! 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)									:: Ncollect
integer, intent(out)									:: Neew
integer, intent(out)									:: NinCycle
integer, intent(out)									:: Nstorage
integer, intent(out)									:: Nadjust
integer, intent(out)									:: Nresmol
integer, intent(out)									:: Nrefresh
integer, intent(out)									:: Nrefrmoves

character*10, intent(out)								:: histtype

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

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

integer, dimension(MaxSp), intent(inout)			:: tag_val
integer, dimension(MaxSp), intent(inout)			:: tag_mol

! Ubins is the number of energy histogram bins.
! Umin is the minimum energy in the energy histogram.
! Umax is the maximum energy in the energy histogram.

integer, intent(inout)								:: Ubins
real, intent(out)									:: Umin
real, intent(out)									:: Umax

! Nmin is the minimum number of particles in the previous simulation.
! Nmax is the maximum number of particles in the previous simulation.

integer, dimension(MaxSp), intent(out)				:: Nmin
integer, dimension(MaxSp), intent(out)				:: Nmax

real, dimension(Nham), intent(out)					:: LNW

! Nstages is the number of stages for determining the weights
! NST_TOT is the number of MC steps in each stage.
! NST_WT is the number of MC steps before averaging in each stage.
! STOP_H is the decision hamiltonian number for which the weights will 
!	be changed after that stage.

integer													:: Nstages
integer, dimension(Nstages), intent(out)				:: NST_TOT
integer, dimension(Nstages), intent(out)				:: NST_WT
integer, dimension(Nstages), intent(out)				:: STOP_H

! 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, repeat
integer										:: ns, nb, nh
integer										:: Nmol1
integer										:: ref1, ref2
integer										:: len
integer, dimension(MaxBeads)				:: REFBEAD
integer, dimension(MaxBeads)				:: BASE
integer, dimension(MaxBeads)				:: BONDS

character*5									:: btype
character*15								:: Species
character*30								:: InFile

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:Nham) 

BETA = 1.0 / BETA

read(20,*)
read(20,*) ZETA(1:Nsp, 1:Nham)
read(20,*)
read(20,*) Volume
read(20,*)
read(20,*) Uwidth
read(20,*)
read(20,*)
read(20,*)
read(20,*)
read(20,*)

do i = 1, Nljgrs

	read(20,*) j, NAMElj(j), MASSlj(j)

end do

read(20,*)

! Sigma is read in from input file.

do i = 1, Nljgrs * Nham

	read(20,*) j, k, EPS(k,k,j), SIG(k,k,j), CP(k,k,j), ALP(k,k,j)

end do

do h = 1, Nham

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

end do

! 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(i,i,h) > 0.0 ) then

				CP(i,j,h) = sqrt( CP(i,i,h) * CP(j,j,h) ) 
				CP(j,i,h) = CP(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 ratio of rmax_ij / rm_ij

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

				RMAX(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

do i = 1, Nljgrs 

	do h = 1, Nham

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

			! 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 ratio of rmax_ii / rm_ii

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

			! Water simulations

!			if( i == 1 ) then

!				RMAX(i,i,h) = 0.45

!			end if

		end if

	end do

end do

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

if( Niongrs > 0 ) then

	read(20,*)

	do i = 1, Niongrs

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

	end do

	read(20,*)

	do i = 1, Niongrs * Nham

		read(20,*) j, k, CHARGE(k,j)

	end do

end if

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

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

do i = 1, Nsp

	read(20,*)
	read(20,*) NAMEsp(i), ERSTEPS(i), EESTEPS(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,*)

	do j = 1, ERSTEPS(i)

		read(20,*) ERSTART(j,i), ERSTART(j,i), EREND(j,i)

	end do

	! read in alpha values for expanded ensemble

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

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

		read(20,*) step, BEADDAMP(step,1:EESTEPS(i),i)

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

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

			case( 'ResRot' )

			case( 'ResCir' )

			case( 'ResCR4' )

			case( 'ResMatch' )

			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
		read(20,*)

	else 

		read(20,*) Species, Nmol1

	end if
	
	do j = 1, Nsp

		if( Species == NAMEsp(j) ) then

			Nmol(j) = Nmol1

		end if

	end do

end do

Nmol(0) = sum( Nmol(1:Nsp) )

if( FromDisk ) then

	read(30,*)
	read(30,*) StartConf, Seed
	read(30,*)
	read(30,*) tag_val(1:Nsp), tag_mol(1:Nsp)

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

else

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

	tag_val = 0
	tag_mol	= 0
 
end if

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

do i = 2, nmoves

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

end do

PROB_MOVE = PROB_MOVE / PROB_MOVE(nmoves)

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_CD(1:Nsp)

do i = 2, Nsp

	PROB_SP_CD(i) = PROB_SP_CD(i) + PROB_SP_CD(i-1)

end do

PROB_SP_CD = PROB_SP_CD / PROB_SP_CD(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,*) NinCycle
read(20,*)
read(20,*) Nstorage
read(20,*)
read(20,*) Nadjust
read(20,*)
read(20,*) Neew
read(20,*) 
read(20,*) Nrefresh
read(20,*) 
read(20,*) Nrefrmoves
read(20,*) 
read(20,*) Ncollect
read(20,*) 
read(20,*) histtype
read(20,*)
read(20,*) 

i = 1

do while( i <= Nstages )

	read(20,*) ns, nb, nh, repeat

	do j = i, i + repeat

		NST_TOT(j) = ns
		NST_WT(j) = nb
		STOP_H(j) = nh

	end do								  

	i = i + 1 + repeat

end do

if( FromDisk ) then

	read(30,*)
	read(30,*)
	read(30,*) DXYZ(1:Nsp)
	read(30,*) DROT(1:Nsp)
	read(30,*)
	read(30,*)
	read(30,*) Ubins, Umin, Umax
	read(30,*)
	read(30,*) Nmin(1), Nmax(1)

else

	DXYZ = 0.1
	DROT = 0.2

end if

if( FromDisk ) then

	read(30,*)
	read(30,*) 

	do i = 1, Nham

		read(30,*) LNW(i)

	end do

	LNW = log(LNW)

	do j = 1, Nsp

		read(30,*)
		read(30,*) 

		do i = 1, EESTEPS(j)

			read(30,*) EEW(i,j)

		end do

	end do

else

!	LNW = -log( real( Nham ) )

	LNW = -1.0e50
	LNW(1) = 0.0

	EEW = 0.0

end if

close(20)
close(30)


BONDSAPART = 0

do i = 1, Nsp

	InFile = 'ba_'//char(48+i)//'.dat'

	open( 33, file = InFile )

	k = LENLJ(i) + LENION(i)

	do j = 1, k

		read(33,*) BONDSAPART(j,1:k,i)

	end do

	close(33)

end do

return

! Code below is no longer used.

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						 

return

end subroutine BigRead