subroutine Create( InputConf, 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 )				  
				  
implicit none

! InputConf has the information for continuing a previous run.

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

! beta is the reciprocal temperature.

real, dimension(Nham), intent(in)						:: BETA

! zeta is a gcmc parameter = f(mu, T)
! zeta = q(T) * exp(beta * mu_B) from Frenkel + Smit

real, dimension(MaxSp,Nham), intent(in)					:: ZETA

! 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_ij parameters for each hamiltonian.

real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: EPS
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: SIG
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: CP
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: ALP
real, dimension(Nljgrs,Nljgrs,Nham), intent(in)			:: RMAX

! CHARGE contains the charge of a bead for each hamiltonian.

real, dimension(Niongrs,Nham), intent(in)				:: CHARGE

! 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.
! 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.
! BONDSAPART indicates how many bonds separate two beads.

integer, dimension(MaxSp), intent(in)					:: NSTEPS
integer, dimension(MaxSp), intent(in)					:: LENLJ 
integer, dimension(MaxSp), intent(in)					:: LENION
character*5, dimension(MaxBeads,MaxSp), intent(in)		:: BEADTYPE
integer, dimension(MaxBeads,MaxSp), intent(in)			:: GROUPTYPE 
real, dimension(MaxBeads,MaxEE,MaxSp), intent(in)		:: BEADDAMP
integer, dimension(MaxSteps,MaxSp), intent(in)			:: NTRIALS
integer, dimension(MaxSteps,MaxSp), intent(in)			:: STEPSTART 
integer, dimension(MaxSteps,MaxSp), intent(in)			:: STEPLENGTH
character*10, dimension(MaxBeads,MaxSp), intent(in)		:: METHOD
integer, dimension(MaxInt,MaxBeads,MaxSp), intent(in)	:: INTPARAM	
real, dimension(MaxReal,MaxBeads,MaxSp), intent(in)		:: REALPARAM
integer, dimension(MaxBeads,MaxBeads,MaxSp), intent(in)	:: BONDSAPART

integer, dimension(MaxSp), intent(in)					:: EESTEPS

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

logical, intent(in)										:: FromDisk

! BoxSize contains the length of the simulation boxes.

real, intent(in)										:: BoxSize

! Nlj is the number of LJ beads in the simulation boxes.
! Nion is the number of ionic beads in the simulation boxes.

integer, intent(out)									:: Nlj
integer, intent(out)									:: Nion

! MaxMol is the maximum number of molecules per box.
! MaxLJ is the maximum number of LJ beads per box.
! MaxIon is the maximum number of ionic beads per box.

integer, intent(in)										:: MaxMol
integer, intent(in)										:: MaxLJ
integer, intent(in)										:: MaxIon

! Nmol is the number of molecules in the simulation box.

integer, dimension(0:MaxSp), intent(in)					:: Nmol

! Xlj, Ylj, and Zlj are the coordinates of the LJ beads.
! Xion, Yion, and Zion are the coordinates of the ionic beads.

real, dimension( MaxLJ ), intent(out)					:: Xlj, Ylj, Zlj
real, dimension( MaxIon ), intent(out)					:: Xion, Yion, Zion

! TYPElj contains the group identity of each LJ bead.
! TYPEion contains the group identity of each ionic bead.

integer, dimension( MaxLJ ), intent(out)				:: TYPElj
integer, dimension( MaxIon ), intent(out)				:: TYPEion

real, dimension( MaxLJ ), intent(out)					:: DAMPlj2, DAMPlj3
real, dimension( MaxIon ), intent(out)					:: DAMPion

! SPECIES contains the species identity of each molecule.
! LENGTHlj contains the number of LJ beads in each molecule.
! LENGTHion	contains the number of ionic beads in each molecule.
! STARTlj contains the starting LJ bead number for each molecule.
! STARTion contains the starting ionic bead number for each molecule.

integer, dimension(MaxMol), intent(out)					:: SPECIES
integer, dimension(MaxMol), intent(out)					:: LENGTHlj, LENGTHion
integer, dimension(MaxMol), intent(out)					:: STARTlj, STARTion

! LNW contains the weight of each hamiltonian.

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

! LNPSI contains the log(psi) for each hamiltonian.

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

! LnPi contains the log(pi).

real, intent(out)										:: LnPi

! XLRCORR and ELRCORR contain parameters for the long range LJ energy.

real, dimension(605), intent(in)						:: XLRCORR, ELRCORR

! Alpha is an Ewald sum parameter, Alpha = kappa * L, for kappa in A + T.

real, intent(in)										:: Alpha

! Kmax is an Ewald sum parameter.
! Nkvec is the number of k-vectors used in the Fourier sum.
! KX, KY, KZ contain the vector identity of the Nkvec vectors.
! CONST contains the constant part of the Fourier summation for a given Nkvec.

integer, intent(in)										:: Kmax
integer, intent(in)										:: Nkvec
integer, dimension(Nkvec), intent(in)					:: KX, KY, KZ
real, dimension(Nkvec), intent(in)						:: CONST

! EXPX contains the value of exp( i*kx*x ) for a given kx and ion.

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

! SUMQEXPV contains the summation of qi*exp(i*(kx*x + ky*y + kz*z)) 
! for a given k-vector and hamiltonian.

complex, dimension(Nkvec, Nham), intent(inout)			:: SUMQEXPV

! SUMQX is the summation of qi * xi for all ions in the box.

real, dimension(Nham), intent(inout)					:: SUMQX, SUMQY, SUMQZ

! ULJ is the LJ energy of the system without the long range correction.
! UFOURIER is the coulombic fourier energy of the system.
! UREAL is the coulombic real energy of the system.
! USURF is the coulombic surface energy of the system.
! USELF_CH is the summation of the square of all the charges.
! USELF_MOL is the self energy of a given molecule.
! ULJ_MOL is the non bonded intramolecular LJ energy of a molecule.
! UION_MOL is the non bonded intramolecular ionic energy of a molecule.
! UINTRA contains the bending and torsional energy of each molecule.

real, dimension(Nham), intent(inout)					:: ULJBOX, ULJLR
real, dimension(Nham), intent(inout)					:: UREAL, UFOURIER
real, dimension(Nham), intent(inout)					:: USURF, USELF_CH
real, dimension(MaxMol,Nham), intent(inout)				:: ULJ_MOL
real, dimension(MaxMol,Nham), intent(inout)				:: USELF_MOL
real, dimension(MaxMol,Nham), intent(inout)				:: UION_MOL
real, dimension(MaxMol), intent(inout)					:: UINTRA

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

integer, intent(in)										:: Nres, Nresmol

integer, dimension(Nres), intent(in)					:: reslen

real, dimension(Nresmol, MaxBeads, Nres), intent(in)	:: Xresmols, Yresmols, Zresmols
real, dimension(Nresmol, Nres), intent(in)				:: Uint_resm
real, dimension(Nresmol, Nham, Nres), intent(in)		:: Ulj_resm, Uion_resm

! f14_lj and f14_ion are damping factors for the 1-4 intramolecular interactions

real, intent(in)										:: f14_lj
real, intent(in)										:: f14_ion

! Seed is the current random number generator seed value.

integer, intent(inout)									:: Seed

real, external											:: ran2

! Local stuff

integer											:: h, i, j, k, m, n
integer											:: ljb, ionb
integer											:: mol2, Nbeads
integer											:: mol,  bead, sp
integer											:: stlj, stion
integer											:: llj, lion
integer											:: endlj, endion
integer											:: Count
integer											:: ljbk, ljbm
integer											:: ionbk, ionbm
integer											:: ntemp1, ntemp2
integer											:: nbendp

integer, dimension(MaxBeads, MaxSp)				:: LJBEAD, IONBEAD
integer, dimension(Nljgrs)						:: NGROUPS
integer, dimension(0:MaxSp)						:: Nmol_new
integer, dimension(MaxLJ)						:: TYPElj_tmp
integer, dimension(MaxIon)						:: TYPEion_tmp

integer, allocatable, dimension(:)				:: TYPElj_new
integer, allocatable, dimension(:)				:: TYPEion_new

logical											:: Ions	= .False.
logical											:: New

real											:: x, y, z
real											:: Uintra_new
real											:: Uvib, Ubend, Utor
real											:: Random
real											:: CoulCombo
real											:: Largest
real											:: rtrial
real											:: Utemp

real, parameter									:: Pi = 3.14159265359
real, parameter									:: ec = 1.60217733e-19
real, parameter									:: eps0 = 8.854187817e-12
real, parameter									:: kB = 1.380658e-23

real, dimension(20)								:: Xc, Yc, Zc
real, dimension(4)								:: ZERO2 = 0.0
real, dimension(Nham)							:: ULJBOX_new
real, dimension(Nham)							:: ULJLR_new
real, dimension(Nham)							:: ULJ_MOL_new
real, dimension(Nham)							:: UREAL_new
real, dimension(Nham)							:: USURF_new
real, dimension(Nham)							:: UFOURIER_new
real, dimension(Nham)							:: USELF_CH_new
real, dimension(Nham)							:: USELF_MOL_new
real, dimension(Nham)							:: UION_MOL_new
real, dimension(Nham)							:: ULJ
real, dimension(Nham)							:: UION
real, dimension(Nham)							:: UTOT
real, dimension(Nham)							:: dULJBOX
real, dimension(Nham)							:: dUREAL
real, dimension(Nham)							:: dUSURF
real, dimension(Nham)							:: dUFOURIER
real, dimension(Nham)							:: dUSELF_CH
real, dimension(Nham)							:: ULJ_MOL_part
real, dimension(Nham)							:: UION_MOL_part
real, dimension(Nham)							:: SUMQX_new 
real, dimension(Nham)							:: SUMQY_new
real, dimension(Nham)							:: SUMQZ_new
real, dimension(Nham)							:: BIGW_new, BIGW_old
real, dimension(Nsp)							:: PROB_SP
real, dimension(MaxLJ)							:: Xlj_tmp, Ylj_tmp, Zlj_tmp
real, dimension(MaxIon)							:: Xion_tmp, Yion_tmp, Zion_tmp
real, dimension(MaxLJ)							:: DAMPlj2_tmp, DAMPlj3_tmp 
real, dimension(MaxIon)							:: DAMPion_tmp
real, dimension(Nljgrs, Nljgrs, Nham)			:: EPS_CP

real, allocatable, dimension(:)					:: Xlj_new, Ylj_new, Zlj_new
real, allocatable, dimension(:)					:: Xion_new, Yion_new, Zion_new
real, allocatable, dimension(:)					:: DAMPlj2_new, DAMPlj3_new
real, allocatable, dimension(:)					:: DAMPion_new

complex, dimension(0:Kmax, MaxIon)				:: EXPX_tmp
complex, dimension(-Kmax:Kmax, MaxIon)			:: EXPY_tmp, EXPZ_tmp
complex, dimension(Nkvec, Nham)					:: SUMQEXPV_new

complex, allocatable, dimension(:,:)			:: EXPX_new, EXPY_new, EXPZ_new
complex, allocatable, dimension(:,:)			:: CZERO1, CZERO2


CoulCombo = ec * ec * 1.0e10 / ( 4.0 * Pi * eps0 * kB )

do i = 1, Nsp

	if( LENION(i) > 0 ) Ions = .True.

end do

LJBEAD = 0
IONBEAD = 0

ljb = 0
ionb = 0

do j = 1, Nsp

	ljb = 0
	ionb = 0

	do i = 1, LENLJ(j) + LENION(j)
	
		if(	BEADTYPE(i,j) == 'LJ' ) then

			ljb = ljb + 1
		
			LJBEAD(i,j) = ljb
		
			IONBEAD(i,j) = ionb
		
		else if( BEADTYPE(i,j) == 'ION' ) then

			ionb = ionb + 1

			IONBEAD(i,j) = ionb

			LJBEAD(i,j) = ljb

		end if

	end do

end do

DAMPlj2 = 1.0
DAMPlj3 = 1.0
DAMPion = 1.0

if( FromDisk ) then

	open( 30, file = InputConf )

	do i = 1, 15 + Nsp + Nham + 2*Nsp + sum(EESTEPS(1:NSp))
		
		read(30,*)

	end do

	Nlj = 0
	Nion = 0

	mol2 = 0

	Nbeads = sum( Nmol(1:Nsp) * ( LENLJ(1:Nsp) + LENION(1:Nsp) ) )
	
	do i = 1, Nbeads

		read(30,*) x, y, z, mol, bead, sp

		if( mol /= mol2 ) then

			SPECIES(mol) = sp
			
			STARTlj(mol) = Nlj + 1
			STARTion(mol) = Nion + 1

			LENGTHlj(mol) = LENLJ(sp)
			LENGTHion(mol) = LENION(sp)

			mol2 = mol

		end if

		if( BEADTYPE(bead,sp) == 'LJ' ) then

			Nlj = Nlj + 1

			Xlj( Nlj ) = x
			Ylj( Nlj ) = y
			Zlj( Nlj ) = z

			TYPElj( Nlj ) = GROUPTYPE(bead, sp)

			if( mol == tag_mol(sp) ) then

				DAMPlj2( Nlj ) = sqrt( BEADDAMP(bead, tag_val(sp), sp) )
				DAMPlj3( Nlj ) = BEADDAMP(bead, tag_val(sp), sp) ** (1.0/3.0)

			end if

		else if( BEADTYPE(bead,sp) == 'ION' ) then

			Nion = Nion + 1

			Xion( Nion ) = x
			Yion( Nion ) = y
			Zion( Nion ) = z

			TYPEion( Nion ) = GROUPTYPE(bead, sp)

			if( mol == tag_mol(sp) ) then

				DAMPion( Nion ) = sqrt( BEADDAMP(bead, tag_val(sp), sp) )

			end if

		end if

	end do

	SUMQX = 0.0
	SUMQY = 0.0
	SUMQZ = 0.0

	ULJBOX = 0.0
	ULJLR = 0.0
	ULJ_MOL = 0.0
	UINTRA = 0.0

	UREAL = 0.0
	UFOURIER = 0.0
	USELF_CH = 0.0
	USURF = 0.0
	USELF_MOL = 0.0
	UION_MOL = 0.0

	EXPX = ( 0.0, 0.0 )
	EXPY = ( 0.0, 0.0 )
	EXPZ = ( 0.0, 0.0 )

	SUMQEXPV = ( 0.0, 0.0 )

	call e6interact( Nlj, Xlj, Ylj, Zlj,  &
					 TYPElj, DAMPlj2, DAMPlj3, &
					 Nmol(0), LENGTHlj, &
					 Nham, Nljgrs, EPS, SIG, CP, &
					 ALP, RMAX, BoxSize, ULJBOX )

	NGROUPS = 0
	
	do j = 1, Nlj

		NGROUPS( TYPElj(j) ) = NGROUPS( TYPElj(j) ) + 1

	end do

	do h = 1, Nham

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

			EPS_CP(:,:,h) = EPS(:,:,h) * CP(:,:,h) * ALP(:,:,h) / &
							( ALP(:,:,h) - 6.0 )	/ 4.0

		else

			EPS_CP(:,:,h) = EPS(:,:,h) * ALP(:,:,h) / ( ALP(:,:,h) - 6.0 ) * &
							( ALP(:,:,h) / 6.0 ) ** ( 6.0 / ( ALP(:,:,h) - 6.0 ) ) / 4.0

		end if

	end do

	call lrcorr( Nljgrs, Nham, BoxSize, NGROUPS, EPS_CP, SIG, XLRCORR, &
				 ELRCORR, ULJLR )

	do j = 1, Nmol(0)

		sp = SPECIES(j)
		stlj = STARTlj(j)

		do k = 1, LENGTHlj(j) + LENGTHion(j) - 1
				
			do m = k + 1, LENGTHlj(j) + LENGTHion(j)

				if( BONDSAPART(k,m,sp) >= 3 .AND. BEADTYPE(k,sp) == 'LJ' .AND. &
					BEADTYPE(m,sp) == 'LJ') then

					ljbk = stlj + LJBEAD(k,sp) - 1
					ljbm = stlj + LJBEAD(m,sp) - 1

					call e6molecule( 1, Xlj(ljbm), Ylj(ljbm), &
									 Zlj(ljbm), TYPElj(ljbm), &
									 DAMPlj2(ljbm), DAMPlj3(ljbm), &
									 1, Xlj(ljbk), Ylj(ljbk), &
									 Zlj(ljbk), TYPElj(ljbk), &
									 DAMPlj2(ljbk), DAMPlj3(ljbk), &
									 Nham, Nljgrs, EPS, SIG, CP, &
									 ALP, RMAX, BoxSize, &
									 ULJ_MOL_part )

					if( BONDSAPART(k,m,sp) == 3 ) ULJ_MOL_part = f14_lj * ULJ_MOL_part 

					ULJ_MOL(j,:) = ULJ_MOL(j,:) + ULJ_MOL_part(:)

				end if

			end do

		end do

	end do

	if( Ions ) then

		call Surf_Move( Nion, Xion, Yion, Zion, &
					    TYPEion, DAMPion, &
						Nham, Niongrs, CHARGE, BoxSize, &
					    SUMQX, SUMQY, SUMQZ, &
						SUMQX_new, SUMQY_new, SUMQZ_new, &
						USURF)

		USURF = USURF * CoulCombo

		SUMQX = SUMQX_new
		SUMQY = SUMQY_new
		SUMQZ = SUMQZ_new

		call Fourier_Move( Nion, Xion, Yion, Zion, &
						   TYPEion, DAMPion, &
						   Nham, Niongrs, CHARGE, BoxSize, &
						   Kmax, Nkvec, KX, KY, KZ, CONST, &
						   EXPX, EXPY, EXPZ, &
						   EXPX_tmp, EXPY_tmp, EXPZ_tmp, &
						   SUMQEXPV, SUMQEXPV_new, UFOURIER )

		UFOURIER = UFOURIER * CoulCombo

		EXPX = EXPX_tmp
		EXPY = EXPY_tmp
		EXPZ = EXPZ_tmp
		
		SUMQEXPV = SUMQEXPV_new

		call RealInteract( Nion, Xion, Yion, Zion, &
						   TYPEion, DAMPion, &
						   Nmol(0), LENGTHion, Nham, &
						   Niongrs, CHARGE, BoxSize, &
						   Alpha, UREAL)
						
		UREAL = UREAL * CoulCombo

		do j = 1, Nmol(0)

			sp = SPECIES(j)
			stion = STARTion(j)
			lion = LENGTHion(j)
			endion = stion + lion - 1
		
			call SelfMolecule( lion, Xion(stion:endion), &
							   Yion(stion:endion), &
							   Zion(stion:endion), &
							   TYPEion(stion:endion), &
							   DAMPion(stion:endion), &
							   Nham, Niongrs, CHARGE, BoxSize, &
							   Alpha, USELF_MOL(j,:) )

			USELF_MOL(j,:) = USELF_MOL(j,:) * CoulCombo

			do k = 1, LENGTHlj(j) + LENGTHion(j) - 1
				
				do m = k + 1, LENGTHlj(j) + LENGTHion(j)

					if( BONDSAPART(k,m,sp) >= 3 .AND. BEADTYPE(k,sp) == 'ION' .AND. &
						BEADTYPE(m,sp) == 'ION') then

						ionbk = stion + IONBEAD(k,sp) - 1
						ionbm = stion + IONBEAD(m,sp) - 1

						call IonMolecule( 1, Xion(ionbm), Yion(ionbm), &
									   	  Zion(ionbm), TYPEion(ionbm), &
										  DAMPion(ionbm), &
									 	  1, Xion(ionbk), Yion(ionbk), &
										  Zion(ionbk), TYPEion(ionbk), &
										  DAMPion(ionbk), &
										  Nham, Niongrs, CHARGE, &
										  BoxSize, UION_MOL_part )

						if( BONDSAPART(k,m,sp) == 3 ) UION_MOL_part = f14_ion * UION_MOL_part 

						UION_MOL(j,:) = UION_MOL(j,:) + UION_MOL_part(:) * CoulCombo

					end if

				end do

			end do

		end do

		do h = 1, Nham
				
			do j = 1, Nion
		
				USELF_CH(h) = USELF_CH(h) + DAMPion(j) * DAMPion(j) * &
											CHARGE( TYPEion(j), h ) * &
											CHARGE( TYPEion(j), h )
				
			end do

			USELF_CH(h) = USELF_CH(h) * &
							Alpha / sqrt(Pi) / BoxSize * CoulCombo

		end do

	end if

	do j = 1, Nmol(0)

		stlj = STARTlj(j)
		stion = STARTion(j)
		sp = SPECIES(j)

		do n = 1, NSTEPS(sp)

			do k = STEPSTART(n,sp), STEPSTART(n,sp) + STEPLENGTH(n,sp) - 1

				ljb = LJBEAD(k,sp)
				ionb = IONBEAD(k,sp)

				select case ( METHOD(k,sp) )

						case( 'ConeTor' )
					
							do m = 1, 3

								if( BEADTYPE( INTPARAM(m,k,sp), sp ) == 'LJ' )	then
						
									Xc(m) = Xlj( stlj + LJBEAD( INTPARAM(m,k,sp), sp ) - 1 )
									Yc(m) = Ylj( stlj + LJBEAD( INTPARAM(m,k,sp), sp ) - 1 )
									Zc(m) = Zlj( stlj + LJBEAD( INTPARAM(m,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m,k,sp), sp ) == 'ION' ) then

									Xc(m) = Xion( stion + IONBEAD( INTPARAM(m,k,sp), sp ) - 1 )
									Yc(m) = Yion( stion + IONBEAD( INTPARAM(m,k,sp), sp ) - 1 )
									Zc(m) = Zion( stion + IONBEAD( INTPARAM(m,k,sp), sp ) - 1 )
						
								end if

							end do

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(4) = Xlj( stlj + ljb - 1 )
								Yc(4) = Ylj( stlj + ljb - 1 )
								Zc(4) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(4) = Xion( stion + ionb - 1 )
								Yc(4) = Yion( stion + ionb - 1 )
								Zc(4) = Zion( stion + ionb - 1 )

							end if

							call Outfold( 4, 0, BoxSize, Xc, Yc, Zc, &
										  ZERO2, ZERO2, ZERO2 )

							call IntraTorsion( Xc, Yc, Zc, INTPARAM(4,k,sp), &
											   REALPARAM(3:8,k,sp), Utor )

						case( 'Stretch' )
					
							if( BEADTYPE( INTPARAM(1,k,sp), sp ) == 'LJ' )	then
						
								Xc(1) = Xlj( stlj + LJBEAD( INTPARAM(1,k,sp), sp ) - 1 )
								Yc(1) = Ylj( stlj + LJBEAD( INTPARAM(1,k,sp), sp ) - 1 )
								Zc(1) = Zlj( stlj + LJBEAD( INTPARAM(1,k,sp), sp ) - 1 )

							else if( BEADTYPE( INTPARAM(1,k,sp), sp ) == 'ION' ) then

								Xc(1) = Xion( stion + IONBEAD( INTPARAM(1,k,sp), sp ) - 1 )
								Yc(1) = Yion( stion + IONBEAD( INTPARAM(1,k,sp), sp ) - 1 )
								Zc(1) = Zion( stion + IONBEAD( INTPARAM(1,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(2) = Xlj( stlj + ljb - 1 )
								Yc(2) = Ylj( stlj + ljb - 1 )
								Zc(2) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(2) = Xion( stion + ionb - 1 )
								Yc(2) = Yion( stion + ionb - 1 )
								Zc(2) = Zion( stion + ionb - 1 )

							end if

							call Outfold( 2, 0, BoxSize, Xc, Yc, Zc, &
										  ZERO2, ZERO2, ZERO2 )

							rtrial = sqrt( ( Xc(2) - Xc(1) ) ** 2 + ( Yc(2) - Yc(1) ) ** 2 + &
										   ( Zc(2) - Zc(1) ) ** 2 )

							Uvib = 0.5 * REALPARAM(1,k,sp) * ( rtrial - REALPARAM(2,k,sp) ) ** 2.0


						case( 'FxBend' )

							! INTPARAM(1,k) = # of bending potentials evaluated
							! INTPARAM(2,k) = bead # from which the current bead is to be grown
						
							if( BEADTYPE( INTPARAM(2,k,sp), sp ) == 'LJ' )	then

								Xc(2) = Xlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Yc(2) = Ylj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Zc(2) = Zlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )

							else if( BEADTYPE( INTPARAM(2,k,sp), sp ) == 'ION' ) then
		
								Xc(2) = Xion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Yc(2) = Yion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Zc(2) = Zion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(3) = Xlj( stlj + ljb - 1 )
								Yc(3) = Ylj( stlj + ljb - 1 )
								Zc(3) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(3) = Xion( stion + ionb - 1 )
								Yc(3) = Yion( stion + ionb - 1 )
								Zc(3) = Zion( stion + ionb - 1 )

							end if

							do m = 1, INTPARAM(1,k,sp)
	
								if( BEADTYPE( INTPARAM(m+2,k,sp), sp ) == 'LJ' )	then
						
									Xc(m+3) = Xlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Yc(m+3) = Ylj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Zc(m+3) = Zlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m+2,k,sp), sp ) == 'ION' ) then

									Xc(m+3) = Xion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Yc(m+3) = Yion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Zc(m+3) = Zion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
						
								end if

							end do

							ntemp1 = 2+INTPARAM(1,k,sp)

							call Outfold( ntemp1, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), &
										  ZERO2, ZERO2, ZERO2 )

							do m = 1, INTPARAM(1,k,sp)

								Xc(1) = Xc(m+3)
								Yc(1) = Yc(m+3)
								Zc(1) = Zc(m+3)

								ntemp1 = 2*m
								ntemp2 = 1+2*m
		
								call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1,k,sp), &
												REALPARAM(ntemp2,k,sp), Utemp )

								Ubend = Ubend + Utemp

							end do

						case( 'FxBendTor' )

							if( BEADTYPE( INTPARAM(4,k,sp), sp ) == 'LJ' )	then
						
								Xc(3) = Xlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Yc(3) = Ylj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Zc(3) = Zlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )

							else if( BEADTYPE( INTPARAM(4,k,sp), sp ) == 'ION' ) then

								Xc(3) = Xion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Yc(3) = Yion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Zc(3) = Zion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(4) = Xlj( stlj + ljb - 1 )
								Yc(4) = Ylj( stlj + ljb - 1 )
								Zc(4) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(4) = Xion( stion + ionb - 1 )
								Yc(4) = Yion( stion + ionb - 1 )
								Zc(4) = Zion( stion + ionb - 1 )

							end if
					
							do m = 1, INTPARAM(1,k,sp)	+ 2 * INTPARAM(2,k,sp)

								if( BEADTYPE( INTPARAM(m+4,k,sp), sp ) == 'LJ' )	then
						
									Xc(m+4) = Xlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Yc(m+4) = Ylj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Zc(m+4) = Zlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m+4,k,sp), sp ) == 'ION' ) then

									Xc(m+4) = Xion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Yc(m+4) = Yion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Zc(m+4) = Zion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
						
								end if

							end do

							ntemp1 = 2 + INTPARAM(1,k,sp)	+ 2 * INTPARAM(2,k,sp)

							call Outfold( ntemp1, 0, BoxSize, &
										  Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 )

							do m = 1, INTPARAM(1,k,sp)

								Xc(2) = Xc(m+4)
								Yc(2) = Yc(m+4)
								Zc(2) = Zc(m+4)

								ntemp1 = 2*m
								ntemp2 = 1+2*m
		
								call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1,k,sp), &
												REALPARAM(ntemp2,k,sp), Utemp )

								Ubend = Ubend + Utemp

							end do

							nbendp = INTPARAM(1,k,sp)

							do m = 1, INTPARAM(2,k,sp)

								Xc(1) = Xc(2*m+nbendp+3)
								Yc(1) = Yc(2*m+nbendp+3)
								Zc(1) = Zc(2*m+nbendp+3)

								Xc(2) = Xc(2*m+nbendp+4)
								Yc(2) = Yc(2*m+nbendp+4)
								Zc(2) = Zc(2*m+nbendp+4)

								if( INTPARAM(3,k,sp) == 1 ) then

									ntemp1 = 2*nbendp-4+6*m
									ntemp2 = 2*nbendp+1+6*m
		
									call IntraTorsion( Xc, Yc, Zc, 1, &
													   REALPARAM(ntemp1:ntemp2,k,sp), &
													   Utemp )

								else if( INTPARAM(3,k,sp) == 2 ) then

									ntemp1 = 2*nbendp-2+4*m
									ntemp2 = 2*nbendp+1+4*m

									call IntraTorsion( Xc, Yc, Zc, 2, &
													   REALPARAM(ntemp1:ntemp2,k,sp), &
													   Utemp )

								end if

								Utor = Utor + Utemp

							end do


						case( 'StBend' )

							! INTPARAM(1,k) = # of bending potentials evaluated
							! INTPARAM(2,k) = bead # from which the current bead is to be grown
						
							if( BEADTYPE( INTPARAM(2,k,sp), sp ) == 'LJ' )	then

								Xc(2) = Xlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Yc(2) = Ylj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Zc(2) = Zlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )

							else if( BEADTYPE( INTPARAM(2,k,sp), sp ) == 'ION' ) then
		
								Xc(2) = Xion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Yc(2) = Yion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Zc(2) = Zion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(3) = Xlj( stlj + ljb - 1 )
								Yc(3) = Ylj( stlj + ljb - 1 )
								Zc(3) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(3) = Xion( stion + ionb - 1 )
								Yc(3) = Yion( stion + ionb - 1 )
								Zc(3) = Zion( stion + ionb - 1 )

							end if

							do m = 1, INTPARAM(1,k,sp)
	
								if( BEADTYPE( INTPARAM(m+2,k,sp), sp ) == 'LJ' )	then
						
									Xc(m+3) = Xlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Yc(m+3) = Ylj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Zc(m+3) = Zlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m+2,k,sp), sp ) == 'ION' ) then

									Xc(m+3) = Xion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Yc(m+3) = Yion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Zc(m+3) = Zion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
						
								end if

							end do

							ntemp1 = 2+INTPARAM(1,k,sp)

							call Outfold( ntemp1, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), &
										  ZERO2, ZERO2, ZERO2 )

							rtrial = sqrt( ( Xc(3) - Xc(2) ) ** 2 + ( Yc(3) - Yc(2) ) ** 2 + &
										   ( Zc(3) - Zc(2) ) ** 2 )

							Uvib = 0.5 * REALPARAM(1,k,sp) * ( rtrial - REALPARAM(2,k,sp) ) ** 2.0

							do m = 1, INTPARAM(1,k,sp)

								Xc(1) = Xc(m+3)
								Yc(1) = Yc(m+3)
								Zc(1) = Zc(m+3)

								ntemp1 = 1+2*m
								ntemp2 = 2+2*m
		
								call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1,k,sp), &
												REALPARAM(ntemp2,k,sp), Utemp )

								Ubend = Ubend + Utemp

							end do

						case( 'StBendTor' )

							if( BEADTYPE( INTPARAM(4,k,sp), sp ) == 'LJ' )	then
						
								Xc(3) = Xlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Yc(3) = Ylj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Zc(3) = Zlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )

							else if( BEADTYPE( INTPARAM(4,k,sp), sp ) == 'ION' ) then

								Xc(3) = Xion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Yc(3) = Yion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Zc(3) = Zion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(4) = Xlj( stlj + ljb - 1 )
								Yc(4) = Ylj( stlj + ljb - 1 )
								Zc(4) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(4) = Xion( stion + ionb - 1 )
								Yc(4) = Yion( stion + ionb - 1 )
								Zc(4) = Zion( stion + ionb - 1 )

							end if
					
							do m = 1, INTPARAM(1,k,sp)	+ 2 * INTPARAM(2,k,sp)

								if( BEADTYPE( INTPARAM(m+4,k,sp), sp ) == 'LJ' )	then
						
									Xc(m+4) = Xlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Yc(m+4) = Ylj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Zc(m+4) = Zlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m+4,k,sp), sp ) == 'ION' ) then

									Xc(m+4) = Xion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Yc(m+4) = Yion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Zc(m+4) = Zion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
						
								end if

							end do

							ntemp1 = 2 + INTPARAM(1,k,sp)	+ 2 * INTPARAM(2,k,sp)

							call Outfold( ntemp1, 0, BoxSize, &
										  Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 )

							rtrial = sqrt( ( Xc(4) - Xc(3) ) ** 2 + ( Yc(4) - Yc(3) ) ** 2 + &
										   ( Zc(4) - Zc(3) ) ** 2 )

							Uvib = 0.5 * REALPARAM(1,k,sp) * ( rtrial - REALPARAM(2,k,sp) ) ** 2.0

							do m = 1, INTPARAM(1,k,sp)

								Xc(2) = Xc(m+4)
								Yc(2) = Yc(m+4)
								Zc(2) = Zc(m+4)

								ntemp1 = 1+2*m
								ntemp2 = 2+2*m
		
								call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1,k,sp), &
												REALPARAM(ntemp2,k,sp), Utemp )

								Ubend = Ubend + Utemp

							end do

							nbendp = INTPARAM(1,k,sp)

							do m = 1, INTPARAM(2,k,sp)

								Xc(1) = Xc(2*m+nbendp+3)
								Yc(1) = Yc(2*m+nbendp+3)
								Zc(1) = Zc(2*m+nbendp+3)

								Xc(2) = Xc(2*m+nbendp+4)
								Yc(2) = Yc(2*m+nbendp+4)
								Zc(2) = Zc(2*m+nbendp+4)

								if( INTPARAM(3,k,sp) == 1 ) then

									ntemp1 = 2*nbendp-3+6*m
									ntemp2 = 2*nbendp+2+6*m
		
									call IntraTorsion( Xc, Yc, Zc, 1, &
													   REALPARAM(ntemp1:ntemp2,k,sp), &
													   Utemp )

								else if( INTPARAM(3,k,sp) == 2 ) then

									ntemp1 = 2*nbendp-1+4*m
									ntemp2 = 2*nbendp+2+4*m

									call IntraTorsion( Xc, Yc, Zc, 2, &
													   REALPARAM(ntemp1:ntemp2,k,sp), &
													   Utemp )

								end if

								Utor = Utor + Utemp

							end do

						case default

							Uvib = 0.0
							Ubend = 0.0
							Utor = 0.0

					end select

				UINTRA(j) = UINTRA(j) + Uvib + Ubend + Utor

			end do

		end do

	end do

	close(30)

else

	Nlj = 0
	Nion = 0

	Nmol_new = 0

	SUMQX = 0.0
	SUMQY = 0.0
	SUMQZ = 0.0

	ULJBOX = 0.0
	ULJLR = 0.0
	ULJ_MOL = 0.0
	UINTRA = 0.0

	UREAL = 0.0
	UFOURIER = 0.0
	USELF_CH = 0.0
	USURF = 0.0
	USELF_MOL = 0.0
	UION_MOL = 0.0

	EXPX = ( 0.0, 0.0 )
	EXPY = ( 0.0, 0.0 )
	EXPZ = ( 0.0, 0.0 )

	SUMQEXPV = ( 0.0, 0.0 )

	PROB_SP(1:Nsp) = real( Nmol(1:Nsp) )

	if( Nmol(0) == 0 ) PROB_SP(Nsp)	= 1.0

	do j = 2, Nsp

		PROB_SP(j) = PROB_SP(j-1) + PROB_SP(j)

	end do

	PROB_SP = PROB_SP / PROB_SP(Nsp)

	do while ( Nmol_new(0) < Nmol(0) )
		
		Random = ran2(Seed)

		sp = 0
		j = 1

		do while ( sp == 0 )

			if( Random < PROB_SP(j) ) sp = j

			j = j + 1

		end do

		if( Nmol_new(sp) == Nmol(sp) ) cycle

		llj = LENLJ(sp)
		lion = LENION(sp)

		allocate( Xlj_new(llj) )
		allocate( Ylj_new(llj) )
		allocate( Zlj_new(llj) )

		allocate( TYPElj_new(llj) )

		allocate( DAMPlj2_new(llj) )
		allocate( DAMPlj3_new(llj) )

		DAMPlj2_new	= 1.0
		DAMPlj3_new	= 1.0
		
		Count = 0

		do j = 1, llj + lion

			if( BEADTYPE(j,sp) == 'LJ' ) then

				Count = Count + 1

				TYPElj_new(Count) = GROUPTYPE(j,sp)

			end if

		end do
		
		if( lion > 0 ) then

			allocate( Xion_new( lion ) )
			allocate( Yion_new( lion ) )
			allocate( Zion_new( lion ) )

			allocate( TYPEion_new( lion ) )

			allocate( DAMPion_new( lion ) )

			allocate( EXPX_new( 0:Kmax, lion ) )
			allocate( EXPY_new( -Kmax:Kmax, lion ) )
			allocate( EXPZ_new( -Kmax:Kmax, lion ) )

			DAMPion_new = 1.0

			Count = 0

			do j = 1, llj + lion

				if( BEADTYPE(j,sp) == 'ION' ) then

					Count = Count + 1

					TYPEion_new(Count) = GROUPTYPE(j,sp)

				end if

			end do
			
		end if
		
		if(	Nmol_new(sp) == 0 ) then

			BigW_old = 1.0
				
			do j = 1, NSTEPS(sp)
				
				BigW_old = BigW_old * NTRIALS(j,sp)

			end do

		else
			
			Mol = int( Nmol_new(sp) * ran2(Seed) ) + 1

			j = 0
			Count = 0

			do while ( Count < Mol )
	
				j = j + 1
	
				if( SPECIES(j) == sp ) Count = Count + 1

			end do


			Mol = j

			stlj = STARTlj(Mol)
			endlj = stlj + llj - 1

			stion = STARTion(Mol)
			endion = stion + lion - 1

			if( stlj == 1 ) then

				Xlj_tmp( 1:Nlj-llj ) = Xlj( endlj+1:Nlj )
				Ylj_tmp( 1:Nlj-llj ) = Ylj( endlj+1:Nlj )
				Zlj_tmp( 1:Nlj-llj ) = Zlj( endlj+1:Nlj )

				TYPElj_tmp( 1:Nlj-llj ) = TYPElj( endlj+1:Nlj )

				DAMPlj2_tmp( 1:Nlj-llj ) = DAMPlj2( endlj+1:Nlj )
				DAMPlj3_tmp( 1:Nlj-llj ) = DAMPlj3( endlj+1:Nlj )

			else

				Xlj_tmp( 1:stlj-1 ) = Xlj( 1:stlj-1 )
				Xlj_tmp( stlj:Nlj-llj ) = Xlj( endlj+1:Nlj )
				Ylj_tmp( 1:stlj-1 ) = Ylj( 1:stlj-1 )
				Ylj_tmp( stlj:Nlj-llj ) = Ylj( endlj+1:Nlj )
				Zlj_tmp( 1:stlj-1 ) = Zlj( 1:stlj-1 )
				Zlj_tmp( stlj:Nlj-llj ) = Zlj( endlj+1:Nlj )

				TYPElj_tmp( 1:stlj-1 ) = TYPElj( 1:stlj-1 )
				TYPElj_tmp( stlj:Nlj-llj ) = TYPElj( endlj+1:Nlj )

				DAMPlj2_tmp( 1:stlj-1 ) = DAMPlj2( 1:stlj-1 )
				DAMPlj2_tmp( stlj:Nlj-llj ) = DAMPlj2( endlj+1:Nlj )

				DAMPlj3_tmp( 1:stlj-1 ) = DAMPlj3( 1:stlj-1 )
				DAMPlj3_tmp( stlj:Nlj-llj ) = DAMPlj3( endlj+1:Nlj )

			end if

			if( lion > 0 ) then

				if( stion == 1 ) then

					Xion_tmp( 1:Nion-lion ) = Xion( endion+1:Nion )
					Yion_tmp( 1:Nion-lion ) = Yion( endion+1:Nion )
					Zion_tmp( 1:Nion-lion ) = Zion( endion+1:Nion )

					TYPEion_tmp( 1:Nion-lion ) = TYPEion( endion+1:Nion )

					DAMPion_tmp( 1:Nion-lion ) = DAMPion( endion+1:Nion )

				else

					Xion_tmp( 1:stion-1 ) = Xion( 1:stion-1 )
					Xion_tmp( stion:Nion-lion ) = Xion( endion+1:Nion )
					Yion_tmp( 1:stion-1 ) = Yion( 1:stion-1 )
					Yion_tmp( stion:Nion-lion ) = Yion( endion+1:Nion )
					Zion_tmp( 1:stion-1 ) = Zion( 1:stion-1 )
					Zion_tmp( stion:Nion-lion ) = Zion( endion+1:Nion )

					TYPEion_tmp( 1:stion-1 ) = TYPEion( 1:stion-1 )
					TYPEion_tmp( stion:Nion-lion ) = TYPEion( endion+1:Nion )

					DAMPion_tmp( 1:stion-1 ) = DAMPion( 1:stion-1 )
					DAMPion_tmp( stion:Nion-lion ) = DAMPion( endion+1:Nion )

				end if

			end if

			Xlj_new( 1:llj ) = Xlj( stlj:endlj )
			Ylj_new( 1:llj ) = Ylj( stlj:endlj )
			Zlj_new( 1:llj ) = Zlj( stlj:endlj )

			call e6molecule( llj, Xlj_new, Ylj_new, Zlj_new, &
							 TYPElj(stlj:endlj), DAMPlj2(stlj:endlj), &
							 DAMPlj3(stlj:endlj), &
							 Nlj-llj, Xlj_tmp, Ylj_tmp, Zlj_tmp, &
							 TYPElj_tmp, DAMPlj2_tmp, DAMPlj3_tmp, &
							 Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
							 BoxSize, dULJBOX )

			ULJBOX_new = ULJBOX - dULJBOX

			NGROUPS = 0
	
			do k = 1, Nlj-llj

				NGROUPS( TYPElj_tmp(k) ) = NGROUPS( TYPElj_tmp(k) ) + 1

			end do

			do h = 1, Nham

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

					EPS_CP(:,:,h) = EPS(:,:,h) * CP(:,:,h) * ALP(:,:,h) / &
									( ALP(:,:,h) - 6.0 )	/ 4.0

				else

					EPS_CP(:,:,h) = EPS(:,:,h) * ALP(:,:,h) / ( ALP(:,:,h) - 6.0 ) * &
									( ALP(:,:,h) / 6.0 ) ** ( 6.0 / ( ALP(:,:,h) - 6.0 ) ) / 4.0

				end if

			end do

			call lrcorr( Nljgrs, Nham, BoxSize, NGROUPS, EPS_CP, SIG, &
						 XLRCORR, ELRCORR, ULJLR_new )

			ULJ_MOL_new = 0.0

			Uintra_new = 0.0

			if( lion > 0 ) then

				Xion_new( 1:lion ) = Xion( stion:endion )
				Yion_new( 1:lion ) = Yion( stion:endion )
				Zion_new( 1:lion ) = Zion( stion:endion )

				call RealMolecule( lion, Xion_new, Yion_new, Zion_new, &
								   TYPEion(stion:endion), &
								   DAMPion(stion:endion), &
								   Nion-lion, Xion_tmp, Yion_tmp, &
								   Zion_tmp, TYPEion_tmp, DAMPion_tmp, &
								   Nham, Niongrs, CHARGE, &
								   BoxSize, Alpha, dUREAL )

				UREAL_new = UREAL - dUREAL * CoulCombo

				call Surf_Move( lion, -Xion_new, -Yion_new, -Zion_new, &
								TYPEion(stion:endion), &
								DAMPion(stion:endion), &
								Nham, Niongrs, CHARGE, BoxSize, &
								SUMQX, SUMQY, SUMQZ, &
								SUMQX_new, SUMQY_new, SUMQZ_new, dUSURF )

				USURF_new = USURF + dUSURF * CoulCombo

				allocate( CZERO1(0:Kmax,lion) )
				allocate( CZERO2(-Kmax:Kmax,lion) )

				CZERO1 = ( 0.0, 0.0 )
				CZERO2 = ( 0.0, 0.0 )

				call Fourier_Move( lion, Xion_new, Yion_new, Zion_new, &
								   TYPEion( stion:endion ), &
								   DAMPion(stion:endion), &
								   Nham, Niongrs, -CHARGE, BoxSize, &
								   Kmax, Nkvec, KX, KY, KZ, CONST, &
								   CZERO1, CZERO2, CZERO2, &
								   EXPX_new, EXPY_new, EXPZ_new, &
								   SUMQEXPV, SUMQEXPV_new, dUFOURIER )

				UFOURIER_new = UFOURIER + dUFOURIER * CoulCombo

				deallocate( CZERO1 )
				deallocate( CZERO2 )

				dUSELF_CH = 0.0

				do h = 1, Nham
					
					do m = stion, endion
	
						dUSELF_CH(h) = dUSELF_CH(h) + DAMPion(m) * DAMPion(m) * &
													CHARGE( TYPEion(m), h ) * &
													CHARGE( TYPEion(m), h )
				
					end do
					
					dUSELF_CH(h) = dUSELF_CH(h) * Alpha * CoulCombo / sqrt(Pi) / BoxSize 

					USELF_CH_new(h) = USELF_CH(h) - dUSELF_CH(h)

				end do

				USELF_MOL_new = 0.0
				UION_MOL_new = 0.0
	
			end if

			New = .False.

			call Grow( NSTEPS(sp), STEPSTART(:,sp), STEPLENGTH(:,sp), &
					   NTRIALS(:,sp), llj, lion, MaxBeads, &
					   METHOD(:,sp), MaxInt, MaxReal, &
					   INTPARAM(:,:,sp), REALPARAM(:,:,sp), &
					   BEADTYPE(:,sp), New, 1, &
					   Xlj_new, Ylj_new, Zlj_new, &
					   TYPElj(stlj:endlj), DAMPlj2(stlj:endlj), &
					   DAMPlj3(stlj:endlj), &
					   Xion_new, Yion_new, Zion_new, &
					   TYPEion(stion:endion), DAMPion(stion:endion ), &
					   Nham, ULJBOX_new, ULJLR_new, ULJ_MOL_new, &
					   UREAL_new, USURF_new, UFOURIER_new, USELF_CH_new, &
					   USELF_MOL_new, UION_MOL_new, Uintra_new, &
					   Niongrs, CHARGE, Alpha, &
					   Kmax, Nkvec, KX, KY, KZ, CONST, &
					   SUMQX_new, SUMQY_new, SUMQZ_new, &
					   EXPX_new, EXPY_new, EXPZ_new, &
					   SUMQEXPV_new, BETA, BIGW_old, LNW, &
					   Nlj-llj, Xlj_tmp, Ylj_tmp, Zlj_tmp, &
					   TYPElj_tmp, DAMPlj2_tmp, DAMPlj3_tmp, &
					   Nion-lion, Xion_tmp, Yion_tmp, Zion_tmp, &
					   TYPEion_tmp, DAMPion_tmp, &
					   BoxSize, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					   XLRCORR, ELRCORR, &
					   Nres, Nresmol, reslen, Xresmols, &
					   Yresmols, Zresmols, Uint_resm, &
					   Ulj_resm, Uion_resm, &
					   BONDSAPART(:,:,sp), f14_lj, f14_ion, Seed )

		end if

		ULJBOX_new = ULJBOX				  
		ULJLR_new = ULJLR				  
		ULJ_MOL_new = 0.0				   
		UREAL_new = UREAL				 
		USURF_new = USURF
		UFOURIER_new = UFOURIER
		USELF_CH_new = USELF_CH				  
		USELF_MOL_new = 0.0
		UION_MOL_new = 0.0
		Uintra_new = 0.0
	
		SUMQX_new = SUMQX
		SUMQY_new = SUMQY
		SUMQZ_new = SUMQZ
	
		SUMQEXPV_new = SUMQEXPV

		New = .True.

		call Grow( NSTEPS(sp), STEPSTART(:,sp), STEPLENGTH(:,sp), &
				   NTRIALS(:,sp), llj, lion, MaxBeads, &  
				   METHOD(:,sp), MaxInt, MaxReal, &
				   INTPARAM(:,:,sp), REALPARAM(:,:,sp), &
				   BEADTYPE(:,sp), New, 1, &
				   Xlj_new, Ylj_new, Zlj_new, &
				   TYPElj_new, DAMPlj2_new, DAMPlj3_new, &
				   Xion_new, Yion_new, Zion_new, &
				   TYPEion_new, DAMPion_new, &
				   Nham, ULJBOX_new, ULJLR_new, ULJ_MOL_new, &
				   UREAL_new, USURF_new, UFOURIER_new, USELF_CH_new, &
				   USELF_MOL_new, USELF_MOL_new, Uintra_new, &
				   Niongrs, CHARGE, Alpha, &
				   Kmax, Nkvec, KX, KY, KZ, CONST, &
				   SUMQX_new, SUMQY_new, SUMQZ_new, &
				   EXPX_new, EXPY_new, EXPZ_new, &
				   SUMQEXPV_new, BETA, BIGW_new, LNW, &
				   Nlj, Xlj, Ylj, Zlj, &
				   TYPElj, DAMPlj2, DAMPlj3, &
				   Nion, Xion, Yion, Zion, &
				   TYPEion, DAMPion, &
				   BoxSize, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
				   XLRCORR, ELRCORR, &
				   Nres, Nresmol, reslen, Xresmols, &
				   Yresmols, Zresmols, Uint_resm, &
				   Ulj_resm, Uion_resm, &
				   BONDSAPART(:,:,sp), f14_lj, f14_ion, Seed )

			
		if( ran2(Seed) < BigW_new(1) / BigW_old(1) ) then	   ! Approximation

			Xlj( Nlj+1:Nlj+llj ) = Xlj_new( 1:llj )
			Ylj( Nlj+1:Nlj+llj ) = Ylj_new( 1:llj )
			Zlj( Nlj+1:Nlj+llj ) = Zlj_new( 1:llj )

			TYPElj( Nlj+1:Nlj+llj ) = TYPElj_new( 1:llj )

			DAMPlj2( Nlj+1:Nlj+llj ) = DAMPlj2_new( 1:llj )
			DAMPlj3( Nlj+1:Nlj+llj ) = DAMPlj3_new( 1:llj )

			ULJBOX = ULJBOX_new
			ULJLR = ULJLR_new
			ULJ_MOL( Nmol_new(0)+1,: ) = ULJ_MOL_new(:)

			if( lion > 0 ) then

				Xion( Nion+1:Nion+lion ) = Xion_new( 1:lion )
				Yion( Nion+1:Nion+lion ) = Yion_new( 1:lion )
				Zion( Nion+1:Nion+lion ) = Zion_new( 1:lion )

				TYPEion( Nion+1:Nion+lion ) = TYPEion_new( 1:lion )

				DAMPion( Nion+1:Nion+lion ) = DAMPion_new( 1:lion )

				UREAL = UREAL_new
				USURF = USURF_new
				UFOURIER = UFOURIER_new
				USELF_CH = USELF_CH_new

				SUMQX = SUMQX_new
				SUMQY = SUMQY_new
				SUMQZ = SUMQZ_new

				SUMQEXPV = SUMQEXPV_new

				EXPX( :, Nion+1:Nion+lion ) = EXPX_new( :, 1:lion )
				EXPY( :, Nion+1:Nion+lion ) = EXPY_new( :, 1:lion )
				EXPZ( :, Nion+1:Nion+lion ) = EXPZ_new( :, 1:lion )

				USELF_MOL( Nmol_new(0)+1,: ) = USELF_MOL_new(:)
				UION_MOL( Nmol_new(0)+1,: ) = UION_MOL_new(:)

			end if

			UINTRA( Nmol_new(0)+1 ) = Uintra_new

			LENGTHlj( Nmol_new(0)+1 ) = llj
			LENGTHion( Nmol_new(0)+1 ) = lion

			STARTlj( Nmol_new(0)+1 ) = Nlj + 1 
			STARTion( Nmol_new(0)+1 ) = Nion + 1

			SPECIES( Nmol_new(0)+1 ) = sp

			Nmol_new(0) = Nmol_new(0) + 1
			Nmol_new(sp) = Nmol_new(sp) + 1

			Nlj = Nlj + llj
			Nion = Nion + lion
			
		end if

		deallocate( Xlj_new )
		deallocate( Ylj_new )
		deallocate( Zlj_new )

		deallocate( TYPElj_new )

		deallocate( DAMPlj2_new )
		deallocate( DAMPlj3_new )

		if( lion > 0 ) then

			deallocate( Xion_new )
			deallocate( Yion_new )
			deallocate( Zion_new )

			deallocate( TYPEion_new )

			deallocate( DAMPion_new )

			deallocate( EXPX_new )
			deallocate( EXPY_new )
			deallocate( EXPZ_new )

		end if

	end do

end if

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

	ULJ(h) = ULJBOX(h) + ULJLR(h)
			   	
		
	UTOT(h) = UION(h) + ULJ(h) + sum( ULJ_MOL( 1:Nmol(0), h) ) + &
			  sum( UION_MOL( 1:Nmol(0), h) )
						
end do

LNPSI = 3.0*Nmol(0)*log(BoxSize) - BETA*UTOT

do i = 1, Nsp

	LNPSI(:) = LNPSI(:) + Nmol(i)*log(ZETA(i,:))

	do k = 1, Nmol(i)

		LNPSI = LNPSI - log(real(k))

	end do

end do

Largest = maxval( LNW + LNPSI )

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


return

end subroutine Create