surf_move.f90

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

subroutine Surf_Move( Ns, DX, DY, DZ, TYPEs, DAMPs, &
					  Nham, Niongrs, CHARGE, BoxSize, &
				      SUMQX, SUMQY, SUMQZ, SUMQX_NEW, &
					  SUMQY_NEW, SUMQZ_NEW, DU_SURF)

implicit none

! This routine calculates the change in the surface potential term of the 
! Ewald sum for a displacement, rotation, transfer, or volume change.

! Ns is the number of ionic beads in a molecule that is displaced, rotated,
! or transferred.  For a volume change Ns is the total number of ionic beads
! in the system.

integer, intent(in)									:: Ns

! The values of DX, DY, and DZ depends on the move type.
!	Displacement				DX(i) = Xnew(i) - Xold(i)
!	Rotation					DX(i) = Xnew(i) - Xold(i)
!	Transfer, Donor phase		DX(i) = - Xold(i)
!	Transfer, Receiver phase	DX(i) = Xnew(i)
!	Volume Change				Recalculate	from Scatch
!	Initialization				DX(i) = X(i) , SUMQX = 0 , and Ns = total number of ions

real, dimension(Ns), intent(in)						:: DX, DY, DZ

! TYPEs contains the group identity of the Ns beads.

integer, dimension(Ns), intent(in)					:: TYPEs

real, dimension(Ns), intent(in)						:: DAMPs

! Nham is the number of hamiltonians.
! Niongrs is the number of ion groups.

integer, intent(in)									:: Nham
integer, intent(in)									:: Niongrs

! CHARGE is an array containing the charge of the Ns sites for each hamiltonian.

real, dimension(Niongrs, Nham), intent(in)			:: CHARGE

! BoxSize is the length of the simulation box.

real, intent(in)									:: BoxSize

! SUMQX is the summation of qi * xi for all ions in the box before a move.

real, dimension(Nham), intent(in)					:: SUMQX, SUMQY, SUMQZ

! SUMQX_NEW is the summation of qi * xi for all ions in the box after a move.

real, dimension(Nham), intent(out)					:: SUMQX_NEW, SUMQY_NEW, SUMQZ_NEW

! DU_SURF is the change in the surface potential after a move.

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

! Local variables

integer												:: i, h
real												:: ch_damp
real, parameter										:: Pi = 3.14159265359




do h = 1, Nham
	
	SUMQX_NEW(h) = SUMQX(h)
	SUMQY_NEW(h) = SUMQY(h)
	SUMQZ_NEW(h) = SUMQZ(h)
	
	do i = 1, Ns

		ch_damp = DAMPs(i) * CHARGE(TYPEs(i), h)
	
		SUMQX_NEW(h) = SUMQX_NEW(h) + ch_damp * DX(i)
		SUMQY_NEW(h) = SUMQY_NEW(h) + ch_damp * DY(i)
		SUMQZ_NEW(h) = SUMQZ_NEW(h) + ch_damp * DZ(i)
	
	end do

	DU_SURF(h) = 2.0 * Pi / ( 3.0 * BoxSize	* BoxSize * BoxSize ) * &
				 ( SUMQX_NEW(h) * SUMQX_NEW(h) - SUMQX(h) * SUMQX(h) + &
				   SUMQY_NEW(h) * SUMQY_NEW(h) - SUMQY(h) * SUMQY(h) + &
				   SUMQZ_NEW(h) * SUMQZ_NEW(h) - SUMQZ(h) * SUMQZ(h) )

end do

end subroutine Surf_Move