subroutine alpha_ch( MaxSp, MaxMol, MaxBeads, MaxEE, Nmol, &
					 Nlj, Nion, Xlj, Ylj, Zlj, &
					 TYPElj, DAMPlj2, DAMPlj3, &
					 Xion, Yion, Zion, &
					 TYPEion, DAMPion, BoxSize, &
					 LENGTHlj, LENGTHion, SPECIES, &
					 STARTlj, STARTion, &
					 BEADTYPE, BEADDAMP, BONDSAPART, &
					 EESTEPS, EEW, &
					 Nham, BETA, ZETA, LNW, LNPSI, LnPi, &
					 Nljgrs, EPS, SIG, CP, ALP, RMAX, Niongrs, &
					 CHARGE, Alpha, Kmax, Nkvec, KX, KY, KZ, &
					 CONST, EXPX, EXPY, EXPZ, SUMQEXPV, &
					 SUMQX, SUMQY, SUMQZ, ULJBOX, ULJ_MOL, & 
					 UFOURIER, UREAL, USURF, USELF_CH, &
					 USELF_MOL, UION_MOL, &
					 tag_val, tag_mol, updown, &
					 f14_lj, f14_ion, &
					 SpID, Success, Seed )

implicit none

! Nmol is the number of molecules in the simulation box.
! MaxSp is the maximum number of species in the system.
! Nlj is the number of LJ beads in the simulation box.
! Nion is the number of ionic beads in the simulation box.

integer, intent(in)									:: MaxSp
integer, intent(in)									:: MaxMol
integer, intent(in)									:: MaxBeads
integer, intent(in)									:: MaxEE
integer, dimension(0:MaxSp), intent(in)				:: Nmol
integer, intent(in)									:: Nlj
integer, intent(in)									:: Nion

! 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), intent(inout)					:: Xlj, Ylj, Zlj
real, dimension(Nion), 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), intent(in)					:: TYPElj
integer, dimension(Nion), intent(in)				:: TYPEion

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

! BoxSize is the length of the simulation box.

real, intent(in)									:: BoxSize

! LENGTHlj contains the number of LJ beads in each molecule.
! LENGTHion	contains the number of ionic beads in each molecule.
! SPECIES contains the species identity of 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)), intent(in)				:: LENGTHlj, LENGTHion
integer, dimension(Nmol(0)), intent(in)				:: SPECIES
integer, dimension(Nmol(0)), intent(in)				:: STARTlj, STARTion

character*5, dimension(MaxBeads, MaxSp)				:: BEADTYPE
real, dimension(MaxBeads, MaxEE, MaxSp)				:: BEADDAMP

! BONDSAPART gives the bondlengths any two LJ beads are away from each other.

integer, dimension(MaxBeads,MaxBeads,MaxSp)			:: BONDSAPART

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

! Nham is the number of hamiltonians being used.

integer, intent(in)									:: Nham

! beta contains the reciprical temperature.

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

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

! 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(inout)				:: LNPSI

! LnPi contains the log(pi).

real, intent(inout)									:: LnPi

! 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), intent(inout)		:: EXPX
complex, dimension(-Kmax:Kmax, Nion), 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

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

real, dimension(Nham), intent(inout)				:: ULJBOX
real, dimension(MaxMol, Nham), intent(inout)		:: ULJ_MOL

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

integer												:: tag_val, tag_mol

integer												:: updown

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

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

! SpID gives the identity of the species that was displaced or rotated.

integer												:: SpID

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

logical, intent(out)								:: Success

! Seed is the current random number generator seed value.

integer, intent(inout)								:: Seed 

real, external										:: ran2

! Local Variables.

integer												:: h, i, j, k, m, Count
integer												:: lenlj, stlj, endlj
integer												:: lenion, stion, endion
integer												:: ljbk, ljbm, ionbk, ionbm
integer, dimension( Nlj )							:: temp4
integer, dimension( Nion )							:: temp8
integer, dimension( MaxBeads )						:: LJBEAD, IONBEAD

real												:: CoulCombo
real												:: Largest, LnPi_new
real												:: LnPi_tmp
real												:: max_damp

real, allocatable, dimension( : )					:: Xlj_old, Ylj_old, Zlj_old
real, allocatable, dimension( : )					:: DAMPlj2_old, DAMPlj3_old
real, allocatable, dimension( : )					:: DAMPlj2_new, DAMPlj3_new
real, dimension( Nlj )								:: temp1, temp2, temp3
real, dimension( Nion )								:: temp5, temp6, temp7
real, dimension( Nlj )								:: temp9, temp10
real, dimension( Nion )								:: temp11
real, allocatable, dimension( : )					:: Xion_old, Yion_old, Zion_old
real, allocatable, dimension( : )					:: DAMPion_old, DAMPion_new
real, allocatable, dimension( : )					:: dDAMPion
real, dimension(Nham)								:: LNPSI_new
real, dimension(Nham)								:: dULJBOX
real, dimension(Nham)								:: dULJ_MOL
real, dimension(Nham)								:: dUREAL
real, dimension(Nham)								:: dUSURF
real, dimension(Nham)								:: dUFOURIER
real, dimension(Nham)								:: dUSELF_CH
real, dimension(Nham)								:: dUSELF_MOL
real, dimension(Nham)								:: dUION_MOL
real, dimension(Nham)								:: dU
real, dimension(Nham)								:: U_part
real, dimension(Nham)								:: ZETA_tmp
real, dimension(Nham)								:: SUMQX_NEW, SUMQY_NEW, SUMQZ_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


Success = .false.

if( tag_mol == 0 ) then

	if( Nmol(SpID) == 0 ) then
	
		tag_val = 0
		return

	end if

	tag_mol = int( Nmol(SpID) * ran2(Seed) ) + 1

	i = 0
	Count = 0

	do while ( Count < tag_mol )
		
		i = i + 1
		
		if( SPECIES(i) == SpID ) Count = Count + 1

	end do

	tag_mol = i

end if

lenlj = LENGTHlj( tag_mol )
stlj = STARTlj( tag_mol )
endlj = stlj + lenlj - 1

lenion = LENGTHion( tag_mol )
stion = STARTion( tag_mol )
endion = stion + lenion - 1

allocate( Xlj_old(lenlj) )
allocate( Ylj_old(lenlj) )
allocate( Zlj_old(lenlj) )

allocate( DAMPlj2_old(lenlj) )
allocate( DAMPlj3_old(lenlj) )

allocate( DAMPlj2_new(lenlj) )
allocate( DAMPlj3_new(lenlj) )

Xlj_old( 1:lenlj ) = Xlj( stlj:endlj )
Ylj_old( 1:lenlj ) = Ylj( stlj:endlj )
Zlj_old( 1:lenlj ) = Zlj( stlj:endlj )

DAMPlj2_old( 1:lenlj ) = DAMPlj2( stlj:endlj )
DAMPlj3_old( 1:lenlj ) = DAMPlj3( stlj:endlj )

! Set DAMPlj2_new values later

if( lenion > 0 ) then

	allocate( Xion_old(lenion) )
	allocate( Yion_old(lenion) )
	allocate( Zion_old(lenion) )

	allocate( DAMPion_old(lenion) )
	allocate( DAMPion_new(lenion) )
	allocate( dDAMPion(lenion) )

	Xion_old( 1:lenion ) = Xion( stion:endion )
	Yion_old( 1:lenion ) = Yion( stion:endion )
	Zion_old( 1:lenion ) = Zion( stion:endion )

	DAMPion_old( 1:lenion ) = DAMPion( stion:endion )

	! Set DAMPion_new values later

end if

LJBEAD = 0
IONBEAD = 0

ljbk = 0
ionbk = 0

do i = 1, lenion + lenlj
	
	if(	BEADTYPE(i,SpID) == 'LJ' ) then

		ljbk = ljbk + 1
		
		LJBEAD(i) = ljbk

		IONBEAD(i) = ionbk
		
		DAMPlj2_new( ljbk ) = sqrt( BEADDAMP( i, tag_val + updown, SpID ) )
		DAMPlj3_new( ljbk ) = BEADDAMP( i, tag_val + updown, SpID )	** (1.0/3.0)

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

		ionbk = ionbk + 1

		IONBEAD(i) = ionbk

		LJBEAD(i) = ljbk

		DAMPion_new( ionbk ) = sqrt( BEADDAMP( i, tag_val + updown, SpID ) )

	end if

end do

ljbk = 0
ionbk = 0

! Calculation of ULJBOX.

if( stlj == 1 )	then

	if( Nmol(0) == 1 ) then
	
		dULJBOX	= 0.0

	else

		call e6molecule2( lenlj, Xlj_old, Ylj_old, Zlj_old, TYPElj(1:lenlj), &
						  DAMPlj2_old, DAMPlj3_old, DAMPlj2_new, DAMPlj3_new, &
						  Nlj - lenlj, Xlj(endlj+1:Nlj), Ylj(endlj+1:Nlj), &						  
						  Zlj(endlj+1:Nlj), TYPElj(endlj+1:Nlj), &
						  DAMPlj2(endlj+1:Nlj),	DAMPlj3(endlj+1:Nlj), &
						  DAMPlj2(endlj+1:Nlj),	DAMPlj3(endlj+1:Nlj), &
						  Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
						  BoxSize, dULJBOX )

	end if

else if( stlj + lenlj - 1 == Nlj ) then

	call e6molecule2( lenlj, Xlj_old, Ylj_old, Zlj_old, TYPElj(stlj:endlj), &
					  DAMPlj2_old, DAMPlj3_old, DAMPlj2_new, DAMPlj3_new, &
					  Nlj - lenlj, Xlj(1:stlj-1), Ylj(1:stlj-1), &						  
					  Zlj(1:stlj-1), TYPElj(1:stlj-1), &
					  DAMPlj2(1:stlj-1), DAMPlj3(1:stlj-1), &
					  DAMPlj2(1:stlj-1), DAMPlj3(1:stlj-1), &
					  Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					  BoxSize, dULJBOX )
	
else

	temp1( 1:stlj-1 ) = Xlj( 1:stlj-1 )
	temp1( stlj:Nlj-lenlj ) = Xlj( endlj+1:Nlj )
	
	temp2( 1:stlj-1 ) = Ylj( 1:stlj-1)
	temp2( stlj:Nlj-lenlj ) = Ylj( endlj+1:Nlj)
	
	temp3( 1:stlj-1 ) = Zlj( 1:stlj-1)
	temp3( stlj:Nlj-lenlj ) = Zlj( endlj+1:Nlj)
	
	temp4( 1:stlj-1 ) = TYPElj( 1:stlj-1)
	temp4( stlj:Nlj-lenlj ) = TYPElj( endlj+1:Nlj)

	temp9( 1:stlj-1 ) = DAMPlj2( 1:stlj-1)
	temp9( stlj:Nlj-lenlj ) = DAMPlj2( endlj+1:Nlj)

	temp10( 1:stlj-1 ) = DAMPlj3( 1:stlj-1)
	temp10( stlj:Nlj-lenlj ) = DAMPlj3( endlj+1:Nlj)

	call e6molecule2( lenlj, Xlj_old, Ylj_old, Zlj_old, TYPElj(stlj:endlj), &
					  DAMPlj2_old, DAMPlj3_old, DAMPlj2_new, DAMPlj3_new, &
					  Nlj - lenlj, temp1, temp2, temp3, temp4, &
					  temp9, temp10, temp9, temp10, &
					  Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					  BoxSize, dULJBOX )

end if

! Calculation of ULJ_MOL.

dULJ_MOL(:) = 0.0

! If the intramolecular nonbonded interactions are included
! as part of the reservoir energy, then comment the 
! following lines.  Also change the appropriate lines in 
! resread.f90, cranksh.f90 and grow.f90. 

do k = 1, lenlj + lenion - 1
		
	do m = k + 1, lenlj + lenion

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

			ljbk = LJBEAD(k)
			ljbm = LJBEAD(m)

			call e6molecule( 1, Xlj_old(ljbm), Ylj_old(ljbm), &
							 Zlj_old(ljbm), TYPElj(stlj+ljbm-1), &
							 DAMPlj2_new(ljbm), DAMPlj3_new(ljbm), &
							 1, Xlj_old(ljbk), Ylj_old(ljbk), &
							 Zlj_old(ljbk), TYPElj(stlj+ljbk-1), &
							 DAMPlj2_new(ljbk), DAMPlj3_new(ljbk), &
							 Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
							 BoxSize, U_part )

			if( BONDSAPART(k,m,SpID) == 3 ) U_part = f14_lj * U_part 

			dULJ_MOL = dULJ_MOL + U_part

		end if

	end do

end do

dULJ_MOL(:) = dULJ_MOL(:) - ULJ_MOL(tag_mol,:)

! Stop commenting lines here.

! Coulombic Energies

if( lenion > 0 ) then

	dDAMPion = DAMPion_new - DAMPion_old

	max_damp = maxval( abs(dDAMPion) )

end if

if( lenion > 0 .AND. max_damp > 1.0e-7 ) then

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

	if( stion == 1 ) then

		if( Nmol(0) == 1 ) then

			dUREAL = 0.0

		else

			call RealMolecule2( lenion, Xion_old, Yion_old, Zion_old, &
								TYPEion(1:lenion), &
								DAMPion_old, DAMPion_new, &
								Nion - lenion, Xion(endion+1:Nion), &
								Yion(endion+1:Nion), Zion(endion+1:Nion), &
								TYPEion(endion+1:Nion), &
								DAMPion(endion+1:Nion), DAMPion(endion+1:Nion), &
								Nham, Niongrs, CHARGE, BoxSize, Alpha, dUREAL )
	
		end if

	else if( stion + lenion - 1 == Nion ) then

		call RealMolecule2( lenion, Xion_old, Yion_old, Zion_old, &
							TYPEion(stion:endion), &
							DAMPion_old, DAMPion_new, &
							Nion - lenion, Xion(1:stion-1), &
							Yion(1:stion-1), Zion(1:stion-1), &
							TYPEion(1:stion-1), &
							DAMPion(1:stion-1), DAMPion(1:stion-1), &
							Nham, Niongrs, CHARGE, BoxSize, Alpha, dUREAL )
	
	else

		temp5( 1:stion-1 ) = Xion( 1:stion-1 )
		temp5( stion:Nion-lenion ) = Xion( endion+1:Nion )
		
		temp6( 1:stion-1 ) = Yion( 1:stion-1 )
		temp6( stion:Nion-lenion ) = Yion( endion+1:Nion )
		
		temp7( 1:stion-1 ) = Zion( 1:stion-1 )
		temp7( stion:Nion-lenion ) = Zion( endion+1:Nion )
		
		temp8( 1:stion-1 ) = TYPEion( 1:stion-1 )
		temp8( stion:Nion-lenion ) = TYPEion( endion+1:Nion )

		temp11( 1:stion-1 ) = DAMPion( 1:stion-1 )
		temp11( stion:Nion-lenion ) = DAMPion( endion+1:Nion )

		call RealMolecule2( lenion, Xion_old, Yion_old, Zion_old, &
							TYPEion(stion:endion), &
							DAMPion_old, DAMPion_new, &
							Nion - lenion, temp5, temp6, temp7, &
							temp8, temp11, temp11, &
							Nham, Niongrs, CHARGE, BoxSize, Alpha, dUREAL )

	end if

	dUREAL = dUREAL	* CoulCombo

	call Fourier_Move2( lenion, TYPEion( stion:endion ), &
						dDAMPion, &
						Nham, Niongrs, CHARGE, BoxSize, &
						Kmax, Nkvec, KX, KY, KZ, CONST, &
						EXPX(:,stion:endion), &
						EXPY(:,stion:endion), &
						EXPZ(:,stion:endion), &
						SUMQEXPV, SUMQEXPV_NEW, &
						dUFOURIER )

	dUFOURIER = dUFOURIER * CoulCombo

	call Surf_Move( lenion, Xion_old, Yion_old, Zion_old, &
					TYPEion( stion:endion ), dDAMPion, &
					Nham, Niongrs, CHARGE, BoxSize, &
					SUMQX, SUMQY, SUMQZ, SUMQX_NEW, &
					SUMQY_NEW, SUMQZ_NEW, dUSURF )

	dUSURF = dUSURF	* CoulCombo

	call SelfMolecule( lenion, Xion_old, Yion_old, Zion_old, &
					   TYPEion(stion:endion), DAMPion_new, &
					   Nham, Niongrs, CHARGE, BoxSize, Alpha, &
					   dUSELF_MOL )

	dUSELF_MOL(:) = dUSELF_MOL(:) * CoulCombo 

	dUSELF_MOL = dUSELF_MOL - USELF_MOL(tag_mol,:)

	dUION_MOL = 0.0

	! If the intramolecular nonbonded interactions are included
	! as part of the reservoir energy, then comment the 
	! following lines.  Also change the appropriate lines in 
	! resread.f90, cranksh.f90 and grow.f90. 

	do k = 1, lenlj + lenion - 1
			
		do m = k + 1, lenlj + lenion

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

				ionbk = IONBEAD(k)
				ionbm = IONBEAD(m)

				call IonMolecule( 1, Xion_old(ionbm), Yion_old(ionbm), &
								  Zion_old(ionbm), TYPEion(stion+ionbm-1), &
								  DAMPion_new(ionbm), &
								  1, Xion_old(ionbk), Yion_old(ionbk), &
								  Zion_old(ionbk), TYPEion(stion+ionbk-1), &
								  DAMPion_new(ionbk), &
								  Nham, Niongrs, CHARGE, BoxSize, U_part )

				if( BONDSAPART(k,m,SpID) == 3 ) U_part = f14_ion * U_part 

				dUION_MOL = dUION_MOL + U_part

			end if

		end do

	end do

	dUION_MOL = dUION_MOL * CoulCombo

	dUION_MOL(:) = dUION_MOL(:) - UION_MOL(tag_mol,:)

	! Stop commenting lines here.

	dUSELF_CH = 0.0
	
	do h = 1, Nham
			
		do j = stion, endion

			dUSELF_CH(h) = dUSELF_CH(h) + &
							CHARGE( TYPEion(j), h ) * CHARGE( TYPEion(j), h ) * &
							( DAMPion_new(j-stion+1) * DAMPion_new(j-stion+1) - &
							  DAMPion_old(j-stion+1) * DAMPion_old(j-stion+1) )
			
		end do

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

	end do

else

	dUREAL = 0.0
	dUSURF = 0.0
	dUFOURIER = 0.0
	dUSELF_MOL = 0.0
	dUION_MOL = 0.0
	dUSELF_CH = 0.0

end if

dU = dULJBOX + dULJ_MOL + &
	 dUREAL + dUSURF + dUFOURIER - dUSELF_CH - dUSELF_MOL + dUION_MOL

ZETA_tmp(:) = ZETA(SpID,:)

LNPSI_new = LNPSI - BETA * dU
				
Largest = maxval( LNW + LNPSI_new )

LnPi_new = log( sum( exp( LNW + LNPSI_new - Largest ) ) ) + Largest

if( updown == 1 ) then

	LnPi_tmp = LnPi_new	+ ( EEW(tag_val + 1, SpID) - EEW(tag_val, SpID)	)

else

	LnPi_tmp = LnPi_new	- ( EEW(tag_val, SpID) - EEW(tag_val - 1, SpID)	)

end if

if( log( ran2(Seed) ) < LnPi_tmp - LnPi ) then		

	Success = .true.

	LnPi = LnPi_new 

	LNPSI = LNPSI_new

	ULJBOX = ULJBOX + dULJBOX

	ULJ_MOL(tag_mol,:) = ULJ_MOL(tag_mol,:) + dULJ_MOL(:)

	DAMPlj2( stlj:endlj ) = DAMPlj2_new( 1:lenlj ) 
	DAMPlj3( stlj:endlj ) = DAMPlj3_new( 1:lenlj )

	if( lenion > 0 ) then

		SUMQEXPV = SUMQEXPV_NEW

		SUMQX = SUMQX_NEW
		SUMQY = SUMQY_NEW
		SUMQZ = SUMQZ_NEW

		UFOURIER = UFOURIER + dUFOURIER
		UREAL = UREAL + dUREAL
		USURF = USURF + dUSURF
		USELF_CH = USELF_CH + dUSELF_CH

		USELF_MOL(tag_mol,:) = USELF_MOL(tag_mol,:) + dUSELF_MOL(:)
		UION_MOL(tag_mol,:)	= UION_MOL(tag_mol,:) + dUION_MOL(:)

		DAMPion( stion:endion ) = DAMPion_new( 1:lenion )

	end if

	tag_val = tag_val + updown

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

	end if

else

	if( tag_val == EESTEPS(SpID) ) then

		tag_val = 0
		tag_mol = 0

	end if

end if

deallocate( Xlj_old )
deallocate( Ylj_old )
deallocate( Zlj_old )

deallocate( DAMPlj2_old )
deallocate( DAMPlj3_old )
deallocate( DAMPlj2_new )
deallocate( DAMPlj3_new )

if( lenion > 0 ) then

	deallocate( Xion_old )
	deallocate( Yion_old )
	deallocate( Zion_old )

	deallocate( DAMPion_old )
	deallocate( DAMPion_new )
	deallocate( dDAMPion )

end if

return

end subroutine alpha_ch