realinteract.f90

巨正则系综蒙特卡罗算法的源程序；可以用来进行吸附等分子模拟；最大的好处在于可以插入或删除原子

subroutine RealInteract( Nb, Xb, Yb, Zb, TYPEb, DAMPb, &
						 Nmol, LENGTH, Nham, Niongrs, &
						 CHARGE, BoxSize, Alpha, ENERGY )

implicit none

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

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

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

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

! 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

! ENERGY contains the total real coulombic 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

	Nc = LENGTH(i)

	if( Nc == 0 ) cycle

	Start = Finish + 1
	Finish = Start + LENGTH(i) - 1

	if( Finish == Nb ) cycle

	call Realmolecule( Nc, Xb(Start:Finish), Yb(Start:Finish), &
					   Zb(Start:Finish), TYPEb(Start:Finish), &
					   DAMPb(Start:Finish), &
					   Nb - Finish, Xb(Finish+1:Nb), Yb(Finish+1:Nb), &
					   Zb(Finish+1:Nb), TYPEb(Finish+1:Nb), &
					   DAMPb(Finish+1:Nb), &
					   Nham, Niongrs, CHARGE, BoxSize, &
					   Alpha, ENERGY_part)

	ENERGY = ENERGY + ENERGY_part

end do

return

end subroutine RealInteract