lrcorr.f90

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

	  
subroutine lrcorr( Nljgrs, Nham, BoxSize, NGROUPS, EPS, SIG, XLRCORR, &
				   ELRCORR, ENERGY )

implicit none

!  This routine calculates the Theodorou and Suter corrections to
!  the LJ energy.  XLRCORR and ELRCORR are read from file lrcorr

! Nljgrs is the number of LJ groups in the simulation.
! Nham is the number of hamiltonians being used.
! NGROUPS contains the number of groups for each Ngljgrs in the simulation box.

integer, intent(in)										:: Nljgrs  
integer, intent(in)										:: Nham
integer, dimension(Nljgrs)								:: NGROUPS
									
! BoxSize is the length of the simulation box.

real, intent(in)										:: BoxSize

! EPS is a rank 3 array containing the eps_ij parameters for each hamiltonian.
! SIG is a rank 3 array containing the sigma_ij parameters for each hamiltonian.

real, dimension(Nljgrs, Nljgrs, Nham), Intent(in)		:: EPS
real, dimension(Nljgrs, Nljgrs, Nham), Intent(in)		:: SIG

! XLRCORR and ELRCORR contain correction factors for the long range LJ energy.

real, dimension(605), intent(in)						:: XLRCORR  
real, dimension(605), intent(in)						:: ELRCORR

! ENERGY contains the long range LJ energy for each hamiltonian.

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

! Local Variables.

integer													:: h, i, j, ref
integer													:: Ncount
real													:: R
real													:: RhoStar
real													:: LStar
real													:: DeltaLStar
real													:: DeltaU
real, dimension(Nljgrs, Nljgrs, Nham)					:: EPSStar
real, dimension(Nljgrs, Nljgrs, Nham)					:: SIGStar
real													:: Pi = 3.14159265359

ref = 0
i = 0

do while ( ref == 0 )
	
	i = i + 1
	if( EPS(i,i,1) /= 0.0 ) ref = i
	
end do
	
EPSStar = EPS / EPS(ref,ref,1)
SIGStar = SIG / SIG(ref,ref,1)

LStar = BoxSize / SIG(ref,ref,1)

ENERGY = 0.0

do i = 1, Nljgrs
	do j = 1, Nljgrs
	
		RhoStar = NGROUPS(i) / ( LStar * LStar * LStar )
	
		do h = 1, Nham

			R = LStar / SIGStar(i, j, h)

			Ncount = int( min( dim(R, 2.0), 3.0 ) * 100 + min( dim(R, 5.0), 2.0 ) * 50 +  &
						  min( dim(R, 7.0), 3.0 ) * 20 + min( dim(R, 10.0), 5.0 ) * 10 +	   &
						  min( dim(R, 15.0), 15.0 ) * 5 + min( dim(R, 30.0), 20.0 ) ) + 1

			DeltaU = 0.0  !  default value

			if (Ncount < 604) then
			
				DeltaLStar = ( R - XLRCORR(Ncount) ) / ( XLRCORR(Ncount+1) - XLRCORR(Ncount) )
			
			    DeltaU = 8.0 * Pi * RhoStar * EPSStar(i,j,h) * &
						 SigStar(i,j,h) * SigStar(i,j,h) * SigStar(i,j,h) *  &
						 ( ELRCORR(Ncount) + ( ELRCORR(Ncount+1) - &
						   ELRCORR(Ncount) ) * DeltaLStar )

				ENERGY(h) = ENERGY(h) + DeltaU * NGROUPS(j) * EPS(ref,ref,1)
				
			end if 

		end do
	end do	
end do

return

end	subroutine lrcorr