Subroutine Add( Nlj, Nion, MaxMol, MaxSp, MaxBeads, MaxEE, Nmol, &
				Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
				Xion, Yion, Zion, TYPEion, DAMPion, &
				MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
				METHOD, MaxInt, MaxReal, INTPARAM, &
				REALPARAM, BEADTYPE, BEADDAMP, BONDSAPART, &
				ERSTEPS, ERSTART, EREND, EESTEPS, EEW, &
				BoxSize, SPECIES, LENGTHlj, LENGTHion, &
				STARTlj, STARTion, &
				LENLJ, LENION, GROUPTYPE_LJ, GROUPTYPE_ION, &
				Nham, BETA, ZETA, LNW, LNPSI, LnPi, XLRCORR, ELRCORR, &
				Nljgrs, EPS, SIG, CP, ALP, RMAX, Niongrs, CHARGE, &
				Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
				EXPX, EXPY, EXPZ, SUMQEXPV, SUMQX, SUMQY, SUMQZ, &
				ULJBOX, ULJLR, UREAL, UFOURIER, USURF, USELF_CH, &
				USELF_MOL, ULJ_MOL, UION_MOL, UINTRA, SpID, &
				Success, tag_val, tag_mol, &
				Nres, Nresmol, reslen, Xresmols, &
				Yresmols, Zresmols, Uint_resm, &
				Ulj_resm, Uion_resm, &
				f14_lj, f14_ion, ER_HIST, Seed )

implicit none

! This subroutine regrows a flexible molecule from StartGrowthStep to the 
! last step and then accepts or rejects the move.

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

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

! MaxMol is the maximum number of molecules per box.
! MaxSp is the maximum number of Species in the simulation.
! MaxBeads in the maximum number of beads per molecule.

integer, intent(in)										:: MaxMol
integer, intent(in)										:: MaxSp
integer, intent(in)										:: MaxBeads
integer, intent(in)										:: MaxEE

! 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), intent(inout)				:: 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(Nlj+MaxBeads), intent(inout)			:: Xlj, Ylj, Zlj
real, dimension(Nion+MaxBeads), intent(inout)			:: Xion, Yion, Zion

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

integer, dimension(Nlj+MaxBeads), intent(inout)			:: TYPElj
integer, dimension(Nion+MaxBeads), intent(inout)		:: TYPEion

real, dimension(Nlj+MaxBeads), intent(inout)			:: DAMPlj2, DAMPlj3
real, dimension(Nion+MaxBeads), intent(inout)			:: DAMPion

! MaxSteps is the maximum number of CB steps.

integer, intent(in)										:: MaxSteps

! STEPSTART is the start bead of each CB step.
! STEPLENGTH is the number of beads in each CB step.
! NTRIALS is the number of trial orientations/locations for each CB step.

integer, dimension(MaxSteps, MaxSp), intent(in)			:: STEPSTART
integer, dimension(MaxSteps, MaxSp), intent(in)			:: STEPLENGTH
integer, dimension(MaxSteps, MaxSp), intent(in)			:: NTRIALS

! METHOD is the method used to grow each bead.
! MaxInt is the maximum number of integer parameters.
! MaxReal is the maximum number of real	parameters.
! INTPARAM are the integer parameters needed by the method to grow each bead.
! REALPARAM are the real parameters needed by the method to grow each bead.
! BEADTYPE indicates whether the bead is 'LJ' or 'ION'.
! BONDSAPART gives the bondlengths any two LJ beads are away from each other.

character*10, dimension(MaxBeads, MaxSp), intent(in)	:: METHOD
integer, intent(in)										:: MaxInt
integer, intent(in)										:: MaxReal
integer, dimension(MaxInt, MaxBeads, MaxSp), intent(in)	:: INTPARAM
real, dimension(MaxReal, MaxBeads, MaxSp), intent(in)	:: REALPARAM
character*5, dimension(MaxBeads, MaxSp), intent(in)		:: BEADTYPE
real, dimension(MaxBeads, MaxEE, MaxSp), intent(in)		:: BEADDAMP
integer, dimension(MaxBeads,MaxBeads,MaxSp), intent(in)	:: BONDSAPART

integer, dimension(MaxSp), intent(in)					:: ERSTEPS
integer, dimension(MaxSteps,MaxSp), intent(in)			:: ERSTART
integer, dimension(MaxSteps,MaxSp), intent(in)			:: EREND

integer, dimension(MaxSp), intent(in)					:: EESTEPS
real, dimension(MaxEE, MaxSp), intent(in)				:: EEW

! BoxSize is the length of the simulation box.

real, intent(in)										:: BoxSize

! 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(Nmol(0)+1), intent(inout)			:: SPECIES
integer, dimension(Nmol(0)+1), intent(inout)			:: LENGTHlj, LENGTHion
integer, dimension(Nmol(0)+1), intent(inout)			:: STARTlj, STARTion

! LENLJ contains the number of LJ beads of each species.
! LENION contains the number of ionic beads of each species.

integer, dimension(MaxSp), intent(in)					:: LENLJ, LENION

! GROUPTYPE_LJ contains the group identity of the LJ beads for a given species.
! GROUPTYPE_ION contains the group identity of the ionic beads for a given species.

integer, dimension(MaxBeads,MaxSp), intent(in)			:: GROUPTYPE_LJ, GROUPTYPE_ION

! Nham is the number of hamiltonians being used.

integer, intent(in)										:: Nham

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

! LNW contains the weight of each hamiltonian.

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

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

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

! LnPi contains the log(pi).

real, intent(inout)										:: LnPi

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

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

! Nljgrs is the number of LJ groups in the system.
! EPS is a rank 3 array containing the eps_ij parameters for each hamiltonian.
! SIG is a rank 3 array containing the sigma_ij parameters for each hamiltonian.
! CP is a rank 3 array containing the C_ij parameters for each hamiltonian.
! ALP is a rank 3 array containing the alpha_ij parameters for each hamiltonian.
! RMAX is a rank 3 array containing the Rmax_ij parameters for each hamiltonian.
									
integer, intent(in)										:: Nljgrs
real, dimension(Nljgrs, Nljgrs, Nham), intent(in)		:: EPS, SIG, CP, ALP, RMAX

! Niongrs is the number of ionic groups in the system.
! CHARGE is a rank 2 array containing the charge of group i for each hamiltonian.
									
integer, intent(in)										:: Niongrs
real, dimension(Niongrs, Nham), intent(in)				:: CHARGE

! 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, Nion+MaxBeads), intent(inout)		:: EXPX
complex, dimension(-Kmax:Kmax, Nion+MaxBeads), intent(inout)	:: 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)				:: USELF_MOL
real, dimension(MaxMol,Nham), intent(inout)				:: ULJ_MOL
real, dimension(MaxMol,Nham), intent(inout)				:: UION_MOL
real, dimension(MaxMol), intent(inout)					:: UINTRA

! SpID indicates which type of species was used in the move.

integer													:: SpID

! Success is a logical indicating whether the move was successful or not.

logical, intent(out)									:: Success

integer													:: tag_val
integer													:: 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

integer, dimension(MaxSteps+1,MaxSp,2), intent(inout)	:: ER_HIST

! Seed is the current random number generator seed value.

integer, intent(inout)									:: Seed

real, external											:: ran2

! Local Variables

integer										:: i, j
integer										:: llj, lion
integer										:: ljbk, ionbk
integer, allocatable, dimension(:) 			:: TYPElj_new
integer, allocatable, dimension(:) 			:: TYPEion_new

logical										:: New

real										:: Largest
real										:: LnPi_tmp, LnPi_old
real										:: Uintra_new
real										:: CoulCombo
real, dimension(Nham)						:: BIGW_new
real, dimension(Nham)						:: TMP, dU
real, dimension(Nham)						:: ULJBOX_new, ULJLR_new
real, dimension(Nham)						:: UREAL_new, UFOURIER_new
real, dimension(Nham)						:: USURF_new, USELF_CH_new
real, dimension(Nham)						:: USELF_MOL_new
real, dimension(Nham)						:: ULJ_MOL_new
real, dimension(Nham)						:: UION_MOL_new
real, dimension(Nham)						:: SUMQX_new, SUMQY_new, SUMQZ_new
real, dimension(Nham)						:: ZETA_tmp

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

real, parameter								:: Pi = 3.14159265359
real, parameter								:: ec = 1.60217733e-19
real, parameter								:: eps0 = 8.854187817e-12
real, parameter								:: kB = 1.380658e-23
										
complex, dimension(Nkvec, Nham)				:: SUMQEXPV_new
complex, allocatable, dimension(:,:)		:: EXPX_new, EXPY_new, EXPZ_new


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

Success = .False.

tag_val = 0
tag_mol = 0

llj = LENLJ(SpID)

lion = LENION(SpID)

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

allocate( TYPElj_new(llj) )

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

Xlj_new = 0.0
Ylj_new = 0.0
Zlj_new = 0.0

TYPElj_new( 1:llj ) = GROUPTYPE_LJ( 1:llj, SpID)

ULJBOX_new = ULJBOX
ULJLR_new = ULJLR
ULJ_MOL_new = 0.0

Uintra_new = 0.0

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

	allocate( DAMPion_new(lion) )

	Xion_new = 0.0
	Yion_new = 0.0
	Zion_new = 0.0

	TYPEion_new( 1:lion ) = GROUPTYPE_ION( 1:lion, SpID)

end if

UREAL_new = UREAL
USURF_new = USURF
UFOURIER_new = UFOURIER
USELF_CH_new = USELF_CH
USELF_MOL_new = 0.0
UION_MOL_new = 0.0

SUMQX_new = SUMQX
SUMQY_new = SUMQY
SUMQZ_new = SUMQZ

SUMQEXPV_new = SUMQEXPV

TMP = LNPSI

LnPi_old = LnPi

ljbk = 0
ionbk = 0

do i = 1, lion + llj
	
	if(	BEADTYPE(i,SpID) == 'LJ' ) then

		ljbk = ljbk + 1

		DAMPlj2_new( ljbk ) = sqrt( BEADDAMP( i, 1, SpID ) )
		DAMPlj3_new( ljbk ) = BEADDAMP( i, 1, SpID ) ** (1.0/3.0)

	else if( BEADTYPE(i,SpID) == 'ION' ) then

		ionbk = ionbk + 1

		DAMPion_new( ionbk ) = sqrt( BEADDAMP( i, 1, SpID ) )

	end if

end do

New = .True.

ZETA_tmp(:) = ZETA(SpID,:)

do i = 1, ERSTEPS(SpID)
	
	call Grow( EREND(i,SpID), STEPSTART(:,SpID), STEPLENGTH(:,SpID), &
			   NTRIALS(:,SpID), llj, lion, MaxBeads, METHOD(:,SpID), &
			   MaxInt, MaxReal, INTPARAM(:,:,SpID), REALPARAM(:,:,SpID), &
			   BEADTYPE(:,SpID), New, ERSTART(i,SpID), &
			   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, 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_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(:,:,SpID), f14_lj, f14_ion, Seed )

	do j = ERSTART(i,SpID), EREND(i,SpID)
	
		BIGW_new = BIGW_new / real( NTRIALS(j,SpID) )

	end do

	if( minval( BIGW_new ) < 1.0e-250 ) then

		LnPi_tmp = -1.0e10

	else

		TMP = TMP + log( BIGW_new ) + &
				log( ZETA_tmp * ( BoxSize**3.0 ) / real( Nmol(SpID) + 1 ) ) / real( ERSTEPS(SpID) )
		
		Largest = maxval( LNW + TMP ) 

		LnPi_tmp = log( sum( exp( LNW + TMP - Largest ) ) ) + Largest

		LnPi_tmp = LnPi_tmp + ( EEW(1,SpID) - EEW(EESTEPS(SpID),SpID) ) / real( ERSTEPS(SpID) )

	end if

	if( ran2(Seed) < 1.0 - exp(LnPi_tmp - LnPi_old) ) exit

	LnPi_old = LnPi_tmp

	if( i == ERSTEPS(SpID) ) then

		Success = .true.

		dU = ULJBOX_new + ULJLR_new + ULJ_MOL_new + &
			 UREAL_new + USURF_new + UFOURIER_new - &
			 USELF_CH_new - USELF_MOL_new + UION_MOL_new - &
			 ( ULJBOX + ULJLR + 0.0 + &
			   UREAL + USURF + UFOURIER - &
			   USELF_CH - 0.0 + 0.0 )

		dU = -dU * BETA

		LNPSI = LNPSI + dU + &
				log( ZETA_tmp * ( BoxSize**3.0 ) / real( Nmol(SpID) + 1 ) )

		Largest = maxval( LNW + LNPSI )

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


		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

		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 )

		end if

		tag_val = tag_val + 1
		tag_mol = Nmol(0) + 1

		if( EESTEPS(SpID) == tag_val ) then
			
			tag_val = 0
			tag_mol = 0

		end if

		SPECIES( Nmol(0)+1 ) = SpID

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

		if( Nmol(0) == 0 ) then

			STARTlj(1) = 1
			STARTion(1) = 1

		else

			STARTlj( Nmol(0)+1 ) = STARTlj( Nmol(0) ) + LENGTHlj( Nmol(0) )
			STARTion( Nmol(0)+1 ) = STARTion( Nmol(0) ) + LENGTHion( Nmol(0) )

		end if

		USELF_MOL( Nmol(0)+1,: ) = USELF_MOL_new(:)

		ULJ_MOL( Nmol(0)+1,: ) = ULJ_MOL_new(:)

		UION_MOL( Nmol(0)+1,: ) = UION_MOL_new(:)

		UINTRA( Nmol(0)+1 ) = Uintra_new

		Nmol(0) = Nmol(0) + 1
		Nmol(SpID) = Nmol(SpID) + 1

		Nlj = Nlj + llj

		Nion = Nion + lion

	end if

end do

ER_HIST(i, SpID, 1) = ER_HIST(i, SpID, 1) + 1

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

return

end subroutine Add