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

integer, parameter										:: NPhases = 2

! InputConf has the information for continuing a previous run.

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.

real, dimension(Ntemp,1), intent(in)					:: BETA

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

! 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, dimension(NPhases), 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, dimension(NPhases), intent(out)				:: Nlj
integer, dimension(NPhases), 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 boxes.
! Nmol(0, iphase) ==> total number of molecules in that phase.
! Nmol(1:MaxSp, iphase) ==> number of molecules of that species in that phase.

integer, dimension(0:MaxSp, 0:NPhases), 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, NPhases ), intent(out)			:: Xlj, Ylj, Zlj
real, dimension( MaxIon, NPhases ), 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, NPhases ), intent(out)		:: TYPElj
integer, dimension( MaxIon, NPhases ), intent(out)		:: TYPEion

! 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, NPhases), intent(out)		:: SPECIES
integer, dimension(MaxMol, NPhases), intent(out)		:: LENGTHlj, LENGTHion
integer, dimension(MaxMol, NPhases), intent(out)		:: STARTlj, STARTion

! LNWSTATE contains the weight of each temperature and hamiltonian.

real, dimension(NTemp, Nham), intent(in)				:: LNWSTATE

! 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, NPhases), intent(out) &
														:: EXPX
complex, dimension(-Kmax:Kmax, MaxIon, NPhases), 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, NPhases), intent(inout) :: SUMQEXPV

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

real, dimension(Nham,NPhases), 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,NPhases), intent(inout)			:: ULJBOX, ULJLR
real, dimension(Nham,NPhases), intent(inout)			:: UREAL, UFOURIER
real, dimension(Nham,NPhases), intent(inout)			:: USURF, USELF_CH
real, dimension(MaxMol,Nham,NPhases), intent(inout)		:: USELF_MOL
real, dimension(MaxMol,Nham,NPhases), intent(inout)		:: ULJ_MOL
real, dimension(MaxMol,Nham,NPhases), intent(inout)		:: UION_MOL
real, dimension(MaxMol,NPhases), intent(inout)			:: UINTRA

! 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, phase, 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, 0:NPhases)			:: 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											:: Ubend, Utor, Uvib, Utemp
real											:: Random
real											:: BigW_new, BigW_old
real											:: CoulCombo
real											:: rtrial

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(15)								:: 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)							:: 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(Nsp)							:: PROB_SP
real, dimension(MaxLJ)							:: Xlj_tmp, Ylj_tmp, Zlj_tmp
real, dimension(MaxIon)							:: Xion_tmp, Yion_tmp, Zion_tmp
real, dimension(Nljgrs,Nljgrs,Nham)				:: EPS_TMP

real, allocatable, dimension(:)					:: Xlj_new, Ylj_new, Zlj_new
real, allocatable, dimension(:)					:: Xion_new, Yion_new, Zion_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

if( FromDisk ) then

	open( 30, file = InputConf )

	do i = 1, 13 + Nsp
		
		read(30,*)

	end do

	Nlj = 0
	Nion = 0

	mol2 = 0

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

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

		if( mol /= mol2 ) then

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

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

			mol2 = mol

		end if

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

			Nlj(phase) = Nlj(phase) + 1

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

			TYPElj( Nlj(phase), phase ) = GROUPTYPE(bead, sp)

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

			Nion(phase) = Nion(phase) + 1

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

			TYPEion( Nion(phase), phase ) = GROUPTYPE(bead, sp)

		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 )


	do i = 1, Nphases

		if( Nmol(0,i) == 0 .AND. BoxSize(i) == 0.0 ) cycle
	
		call e6interact( Nlj(i), Xlj(:,i), Ylj(:,i), Zlj(:,i),  &
						 TYPElj(:,i), Nmol(0,i), LENGTHlj(:,i), &
						 Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
						 BoxSize(i), ULJBOX(:,i) )

		NGROUPS = 0
	
		do j = 1, Nlj(i)

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

		end do

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

			EPS_TMP = EPS * CP * ALP / ( ALP - 6.0 ) / 4.0

		else

			EPS_TMP = EPS

		end if

		call lrcorr( Nljgrs, Nham, BoxSize(i), NGROUPS, EPS_TMP, SIG, XLRCORR, &
					 ELRCORR, ULJLR(:,i) )

		do j = 1, Nmol(0,i)

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

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

					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,i), Ylj(ljbm,i), &
										 Zlj(ljbm,i), TYPElj(ljbm,i), &
										 1, Xlj(ljbk,i), Ylj(ljbk,i), Zlj(ljbk,i), &
										 TYPElj(ljbk,i), Nham, Nljgrs, EPS, &
										 SIG, CP, ALP, RMAX, &
										 BoxSize(i), ULJ_MOL_part )

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

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

					end if

				end do

			end do

		end do

		if( Ions ) then

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

			USURF(:,i) = USURF(:,i) * CoulCombo

			SUMQX(:,i) = SUMQX_new(:)
			SUMQY(:,i) = SUMQY_new(:)
			SUMQZ(:,i) = SUMQZ_new(:)

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

			UFOURIER(:,i) = UFOURIER(:,i) * CoulCombo

			EXPX(:,:,i) = EXPX_tmp(:,:)
			EXPY(:,:,i) = EXPY_tmp(:,:)
			EXPZ(:,:,i) = EXPZ_tmp(:,:)
			
			SUMQEXPV(:,:,i) = SUMQEXPV_new(:,:)

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

			do j = 1, Nmol(0,i)

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

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

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

						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,i), Yion(ionbm,i), &
										   	  Zion(ionbm,i), TYPEion(ionbm,i), &
										 	  1, Xion(ionbk,i), Yion(ionbk,i), Zion(ionbk,i), &
											  TYPEion(ionbk,i), Nham, Niongrs, &
											  CHARGE, BoxSize(i), UION_MOL_part )

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

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

						end if

					end do

				end do

			end do

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

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

			end do

		end if

		do j = 1, Nmol(0,i)

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

			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)

					Uvib = 0.0
					Ubend = 0.0
					Utor = 0.0

					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, i )
									Yc(m) = Ylj( stlj + LJBEAD( INTPARAM(m,k,sp), sp ) - 1, i )
									Zc(m) = Zlj( stlj + LJBEAD( INTPARAM(m,k,sp), sp ) - 1, i )

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

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

							end do

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

							end if

							call Outfold( 4, 0, BoxSize(i), 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, i )
								Yc(1) = Ylj( stlj + LJBEAD( INTPARAM(1,k,sp), sp ) - 1, i )
								Zc(1) = Zlj( stlj + LJBEAD( INTPARAM(1,k,sp), sp ) - 1, i )

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

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

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

							end if

							call Outfold( 2, 0, BoxSize(i), 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, i )
								Yc(2) = Ylj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1, i )
								Zc(2) = Zlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1, i )

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

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

							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, i )
									Yc(m+3) = Ylj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1, i )
									Zc(m+3) = Zlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1, i )

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

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

							end do

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

							call Outfold( ntemp1, 0, BoxSize(i), 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, i )
								Yc(3) = Ylj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1, i )
								Zc(3) = Zlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1, i )

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

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

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

							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, i )
									Yc(m+4) = Ylj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1, i )
									Zc(m+4) = Zlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1, i )

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

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

							end do

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

							call Outfold( ntemp1, 0, BoxSize(i), &
										  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, i )
								Yc(2) = Ylj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1, i )
								Zc(2) = Zlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1, i )

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

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

							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, i )
									Yc(m+3) = Ylj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1, i )
									Zc(m+3) = Zlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1, i )

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

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

							end do

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

							call Outfold( ntemp1, 0, BoxSize(i), 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, i )
								Yc(3) = Ylj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1, i )
								Zc(3) = Zlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1, i )

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

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

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

							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, i )
									Yc(m+4) = Ylj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1, i )
									Zc(m+4) = Zlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1, i )

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

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

							end do

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

							call Outfold( ntemp1, 0, BoxSize(i), &
										  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,i) = UINTRA(j,i) + Uvib + Ubend + Utor

				end do

			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 )

	do i = 1, Nphases

		if( Nmol(0,i) == 0 ) exit

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

		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,i) < Nmol(0,i) )
		
			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,i) == Nmol(sp,i) ) cycle

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

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

			allocate( TYPElj_new(llj) )

			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( EXPX_new( 0:Kmax, lion ) )
				allocate( EXPY_new( -Kmax:Kmax, lion ) )
				allocate( EXPZ_new( -Kmax:Kmax, lion ) )

				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,i) == 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,i) * ran2(Seed) ) + 1

				j = 0
				Count = 0

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

				end do


				Mol = j

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

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

				if( stlj == 1 ) then

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

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

				else

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

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

				end if

				if( lion > 0 ) then

					if( stion == 1 ) then

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

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

					else

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

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

					end if

				end if

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

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

				ULJBOX_new(:) = ULJBOX(:,i) - dULJBOX(:)

				NGROUPS = 0
	
				do k = 1, Nlj(i)-llj

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

				end do

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

					EPS_TMP = EPS * CP * ALP / ( ALP - 6.0 ) / 4.0
			
				else

					EPS_TMP = EPS

				end if

				call lrcorr( Nljgrs, Nham, BoxSize, NGROUPS, EPS_TMP, 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, i )
					Yion_new( 1:lion ) = Yion( stion:endion, i )
					Zion_new( 1:lion ) = Zion( stion:endion, i )

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

					UREAL_new(:) = UREAL(:,i) - dUREAL(:) * CoulCombo

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

					USURF_new(:) = USURF(:,i) + 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, i ), Nham, &
									   Niongrs, -CHARGE, BoxSize(i), Kmax, &
									   Nkvec, KX, KY, KZ, CONST, CZERO1, &
									   CZERO2, CZERO2, EXPX_new, EXPY_new, &
									   EXPZ_new, SUMQEXPV(:,:,i), &
									   SUMQEXPV_new, dUFOURIER )

					UFOURIER_new(:) = UFOURIER(:,i) + dUFOURIER(:) * CoulCombo

					deallocate( CZERO1 )
					deallocate( CZERO2 )

					dUSELF_CH = 0.0

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

						USELF_CH_new(h) = USELF_CH(h,i) - 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, i ), Xion_new, Yion_new, &
						   Zion_new, TYPEion( stion:endion, i ), 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, BigW_old, NTemp, BETA, LNWSTATE, &
						   Nlj(i)-llj, Xlj_tmp, Ylj_tmp, Zlj_tmp, TYPElj_tmp, &
						   Nion(i)-lion, Xion_tmp, Yion_tmp, Zion_tmp, TYPEion_tmp, &
						   BoxSize(i), Nljgrs, EPS, SIG, CP, ALP, RMAX, &
						   XLRCORR, ELRCORR, BONDSAPART(:,:,sp), &
						   f14_lj, f14_ion, Seed )

			end if

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

			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, &
					   Xion_new, Yion_new, Zion_new, TYPEion_new, &
					   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, &
					   BigW_new, NTemp, BETA, LNWSTATE, &
					   Nlj(i), Xlj(:,i), Ylj(:,i), Zlj(:,i), TYPElj(:,i), &
					   Nion(i), Xion(:,i), Yion(:,i), Zion(:,i), TYPEion(:,i), &
					   BoxSize(i), Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					   XLRCORR, ELRCORR, BONDSAPART(:,:,sp), &
					   f14_lj, f14_ion, Seed )

			
			if( ran2(Seed) < BigW_new / BigW_old .OR. Nmol_new(0,i) == 0 ) then

				write(*,*) Nmol_new(0,i) + 1

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

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

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

				if( lion > 0 ) then

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

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

					UREAL(:,i) = UREAL_new(:)
					USURF(:,i) = USURF_new(:)
					UFOURIER(:,i) = UFOURIER_new(:)
					USELF_CH(:,i) = USELF_CH_new(:)

					SUMQX(:,i) = SUMQX_new(:)
					SUMQY(:,i) = SUMQY_new(:)
					SUMQZ(:,i) = SUMQZ_new(:)

					SUMQEXPV(:,:,i) = SUMQEXPV_new(:,:)

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

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

				end if

				UINTRA( Nmol_new(0,i)+1,i ) = Uintra_new

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

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

				SPECIES( Nmol_new(0,i)+1,i ) = sp

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

				Nlj(i) = Nlj(i) + llj
				Nion(i) = Nion(i) + lion
				
			end if

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

			deallocate( TYPElj_new )

			if( lion > 0 ) then

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

				deallocate( TYPEion_new )

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

			end if

		end do

	end do

end if

return

end subroutine Create