subroutine e6molecule( Nc, Xc, Yc, Zc, TYPEc, Np, Xp, Yp, Zp, &
					   TYPEp, Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					   BoxSize, ENERGY )

implicit none

! This routine calculates the energy of Nc LJ beads interacting
! with Np LJ beads.

! Nc is the number of LJ beads in one group.
! TYPEc contains the group identity of the Nc beads.
! Xc, Yc, Zc are the coordinates of the Nc beads.

integer, intent(in)									:: Nc
integer, dimension(Nc), intent(in)					:: TYPEc
real, dimension(Nc), intent(in)						:: Xc, Yc, Zc

! Np is the number of LJ beads in another group.
! TYPEp contains the group identity of the Np beads.
! Xp, Yp, Zp are the coordinates of the Np beads.
									
integer, intent(in)									:: Np
real, dimension(Np), intent(in)						:: Xp, Yp, Zp
integer, dimension(Np), intent(in)					:: TYPEp

! Nham is the number of hamiltonians.
! 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)									:: Nham
integer, intent(in)									:: Nljgrs
real, dimension(Nljgrs, Nljgrs, Nham), intent(in)	:: EPS, SIG, CP, ALP, RMAX

! BoxSize is the length of the simulation box.

real, intent(in)									:: BoxSize

! ENERGY contains the energy of interaction between group c
! and group p for each hamiltonian.

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

! Local variables

integer												:: i, j, h
integer												:: typeci, typepj
real												:: xij, yij, zij
real												:: xi, yi, zi
real												:: rij2, rij
real												:: sigor2, sigor6
real												:: sigij, cpij
real												:: alphaij, epsij
real												:: factor


ENERGY = 0.0

do i = 1, Nc

	xi = Xc(i)
	yi = Yc(i)
	zi = Zc(i)

	typeci = TYPEc(i)

	do j = 1, Np

		typepj = TYPEp(j)

		xij = abs( Xp(j) - xi )
		yij = abs( Yp(j) - yi )
		zij = abs( Zp(j) - zi )

		if( xij > BoxSize - xij ) xij = xij - BoxSize
		if( yij > BoxSize - yij ) yij = yij - BoxSize
		if( zij > BoxSize - zij ) zij = zij - BoxSize

		rij2 = xij*xij + yij*yij + zij*zij

		do h = 1, Nham

			cpij = CP( typepj, typeci, h )

			if( cpij > 0.0 ) then ! Exp-6 fluid

				rij = sqrt( rij2 )

				if( rij  < RMAX( typepj, typeci, h ) ) then

					ENERGY(h) = 1.0e300
					return

				else
				
					alphaij = ALP( typepj, typeci, h )

					epsij = EPS( typepj, typeci, h )

					sigij = SIG( typepj, typeci, h )

					sigor2 = ( sigij * sigij ) / rij2

					ENERGY(h) = ENERGY(h) + epsij / ( alphaij - 6.0 ) * &
						( 6.0 * exp( alphaij * ( 1.0 - rij / sigij ) ) - &
						  cpij * alphaij * sigor2 * sigor2 * sigor2 )

				end if

			else ! LJ fluid

!					alphaij = ALP( typepj, typeci, h )
					alphaij = 12

					epsij = EPS( typepj, typeci, h )

					sigij = SIG( typepj, typeci, h )

					sigor2 = ( sigij * sigij ) / rij2

					sigor6 = sigor2 * sigor2 * sigor2

!					factor = alphaij * epsij / ( alphaij - 6.0 ) * &
!							 ( alphaij / 6.0 ) ** ( 6.0 / ( alphaij -6.0 )  )
					factor = 4.0 * epsij 

					ENERGY(h) = ENERGY(h) + factor * &
!								( sigor2 ** ( alphaij / 2.0 ) - sigor6 )
								( sigor6 * sigor6 - sigor6 )
								
			end if
				  
		end do

	end do

end do

return

end	subroutine e6molecule