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

implicit none

! This routine calculates the real Ewald sum energy of Nc ionic beads 
! interacting with Np ionic 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

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

real, intent(in)							:: Alpha

! erfc is the complementary error function.

real, external								:: erfc

! 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
real										:: tmp


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)

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

		do h = 1, Nham

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

		end do

	end do

end do

return

end	subroutine RealMolecule

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine RealMolecule2( Nc, Xc, Yc, Zc, TYPEc, DAMPc_1, DAMPc_2, &
						  Np, Xp, Yp, Zp, TYPEp, DAMPp_1, DAMPp_2, &
						  Nham, Niongrs, CHARGE, BoxSize, Alpha, dENERGY )

implicit none

! This routine calculates the real Ewald sum energy of Nc ionic beads 
! interacting with Np ionic 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_1,	DAMPc_2

! 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_1,	DAMPp_2

! 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

! erfc is the complementary error function.

real, external								:: erfc

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

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

! Local variables

integer										:: i, j, h
integer										:: typeci, typepj
real										:: xij, yij, zij
real										:: xi, yi, zi
real										:: rij
real										:: dampci_1, damppj_1
real										:: dampci_2, damppj_2
real										:: tmp1, tmp2

real, dimension(Nham)						:: ENERGY_1, ENERGY_2


ENERGY_1 = 0.0
ENERGY_2 = 0.0

do i = 1, Nc

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

	typeci = TYPEc(i)

	dampci_1 = DAMPc_1(i)
	dampci_2 = DAMPc_2(i)

	do j = 1, Np

		typepj = TYPEp(j) 

		damppj_1 = DAMPp_1(j) 
		damppj_2 = DAMPp_2(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)

		tmp1 = erfc( Alpha * rij / BoxSize ) / rij
			  
		do h = 1, Nham

			tmp2 = CHARGE( typeci, h ) * CHARGE( typepj, h ) * tmp1

			ENERGY_1(h) = ENERGY_1(h) + tmp2 * dampci_1 * damppj_1

			ENERGY_2(h) = ENERGY_2(h) + tmp2 * dampci_2 * damppj_2

		end do

	end do

end do

dENERGY = ENERGY_2 - ENERGY_1

return

end	subroutine RealMolecule2