subroutine SelfMolecule( Nc, Xc, Yc, Zc, TYPEc, DAMPc, &
						 Nham, Niongrs, CHARGE, &
						 BoxSize, Alpha, ENERGY)

implicit none

! This routine calculates the self Ewald sum energy of a molecule 
! with Nc ionic beads.

! The second term in A + T, equation 5.24 

! Nc is the number of ionic beads in a molecule.
! 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

! 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

! Alpha is an Ewald sum parameter, Alpha = kappa * L, for kappa in A + T.

real, intent(in)							:: Alpha

! erf is the complementary error function.

real, external								:: erf

! ENERGY contains the self energy of the molecule.

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

! Local variables

integer										:: i, j, h
integer										:: typeci, typecj

real										:: xij, yij, zij
real										:: rij
real										:: tmp
real										:: dampci, dampcj

ENERGY = 0.0

do i = 1, Nc - 1

	typeci = TYPEc(i)
	dampci = DAMPc(i)

	do j = i + 1, Nc
		
		typecj = TYPEc(j)
		dampcj = DAMPc(j)

		xij = abs( Xc(j) - Xc(i) )
		yij = abs( Yc(j) - Yc(i) )
		zij = abs( Zc(j) - Zc(i) )

		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)

		tmp = erf( Alpha * rij / BoxSize ) / rij

		do h = 1, Nham

			ENERGY(h) = ENERGY(h) + dampci * dampcj * tmp * &
						CHARGE( typeci, h	) * CHARGE( typecj, h ) 

		end do

	end do

end do

return

end	subroutine SelfMolecule