subroutine IonMolecule( Nc, Xc, Yc, Zc, TYPEc, DAMPc, &
						Np, Xp, Yp, Zp, TYPEp, DAMPp, &
						Nham, Niongrs, CHARGE, BoxSize, ENERGY )

implicit none

! This routine calculates the coulombic energy between two groups of beads. 

! Nc is the number of ionic 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
real, dimension(Nc), intent(in)				:: DAMPc

! Np is the number of ionic 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
integer, dimension(Np), intent(in)			:: TYPEp
real, dimension(Np), intent(in)				:: Xp, Yp, Zp
real, dimension(Np), intent(in)				:: DAMPp

! Nham is the number of hamiltonians.
! 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)							:: Nham
integer, intent(in)							:: Niongrs
real, dimension(Niongrs, Nham), intent(in)	:: CHARGE

! 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										:: rij
real										:: dampci, damppj


ENERGY = 0.0

do i = 1, Nc

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

	typeci = TYPEc(i)

	dampci = DAMPc(i)

	do j = 1, Np

		typepj = TYPEp(j) 
		
		damppj = DAMPp(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

		rij = sqrt(xij*xij + yij*yij + zij*zij)

		do h = 1, Nham

			ENERGY(h) = ENERGY(h) + CHARGE( typeci, h	) * CHARGE( typepj, h ) &
						* dampci * damppj / rij

		end do

	end do

end do

return

end	subroutine IonMolecule