subroutine e6interact( Nb, Xb, Yb, Zb, TYPEb, Nmol, LENGTH, &
					   Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					   BoxSize, ENERGY )

implicit none

! This routine calculates the total energy of a group of molecules.

! Nb is the total number of LJ beads in a system.
! TYPEb contains the group identity of each bead.
! Xb, Yb, Zb are the coordinates of each bead.

integer, intent(in)									:: Nb
integer, dimension(Nb), intent(in)					:: TYPEb
real, dimension(Nb), intent(in)						:: Xb, Yb, Zb

! Nmol is the number of molecules in the system.
! LENGTH contains the number of LJ beads in each molecule.

integer, intent(in)									:: Nmol
integer, dimension(Nmol), intent(in)				:: LENGTH

! 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 total energy of the system for each hamiltonian.

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

! Local variables

integer												:: i, Nc
integer												:: Start, Finish
real, dimension(Nham)								:: ENERGY_part



ENERGY = 0.0

Finish = 0

do i = 1, Nmol - 1
	Start = Finish + 1
	Finish = Start + LENGTH(i) - 1

	Nc = LENGTH(i)

	call e6molecule( Nc, Xb(Start:Finish), Yb(Start:Finish), &
					 Zb(Start:Finish), TYPEb(Start:Finish), &
					 Nb - Finish, Xb(Finish+1:Nb), Yb(Finish+1:Nb), &
					 Zb(Finish+1:Nb), TYPEb(Finish+1:Nb), Nham, &
					 Nljgrs, EPS, SIG, CP, ALP, RMAX, BoxSize, ENERGY_part)

	ENERGY = ENERGY + ENERGY_part

end do

return

end subroutine e6interact