e6molecule.f90

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

subroutine e6molecule( Nc, Xc, Yc, Zc, TYPEc, DAMP2c, DAMP3c, &
					   Np, Xp, Yp, Zp, TYPEp, DAMP2p, DAMP3p, &
					   Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
					   BoxSize, ENERGY )

implicit none

! This routine calculates the energy of Nc LJ beads interacting
! with Np LJ beads.

! Nc is the number of LJ 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)						:: DAMP2c, DAMP3c

! Np is the number of LJ 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)						:: DAMP2p, DAMP3p

! Nham is the number of hamiltonians.
! Nljgrs is the number of LJ groups in the system.
! 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.
! CP is a rank 3 array containing the C_ij parameters for each hamiltonian.
! ALP is a rank 3 array containing the alpha_ij parameters for each hamiltonian.
! RMAX is a rank 3 array containing the Rmax_ij parameters for each hamiltonian.
									
integer, intent(in)									:: Nham
integer, intent(in)									:: Nljgrs
real, dimension(Nljgrs, Nljgrs, Nham), intent(in)	:: EPS, SIG, CP, ALP, RMAX

! BoxSize is the length of the simulation box.

real, intent(in)									:: BoxSize

! 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												:: rij2, rij
real												:: sigor2, sigor6
real												:: sigij, cpij
real												:: alphaij, epsij
real												:: damp2ci, damp3ci
real												:: damp2pj, damp3pj
real												:: factor


ENERGY = 0.0

do i = 1, Nc

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

	typeci = TYPEc(i)

	damp2ci = DAMP2c(i)
	damp3ci = DAMP3c(i)

	do j = 1, Np

		typepj = TYPEp(j)

		damp2pj = DAMP2p(j)
		damp3pj = DAMP3p(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

		rij2 = xij*xij + yij*yij + zij*zij

		do h = 1, Nham

			cpij = CP( typepj, typeci, h )

			if( cpij > 0.0 ) then ! Exp-6 fluid

				rij = sqrt( rij2 )

				sigij = 0.5 * ( damp3ci + damp3pj ) * SIG( typepj, typeci, h )

				epsij = damp2ci * damp2pj * EPS( typepj, typeci, h )

				if( rij  < RMAX( typepj, typeci, h ) * sigij ) then

					ENERGY(h) = 1.0e300	* epsij

				else
				
					alphaij = ALP( typepj, typeci, h )

					sigor2 = ( sigij * sigij ) / rij2

					ENERGY(h) = ENERGY(h) + epsij / ( alphaij - 6.0 ) * &
						( 6.0 * exp( alphaij * ( 1.0 - rij / sigij ) ) - &
						  cpij * alphaij * sigor2 * sigor2 * sigor2 )

				end if

			else ! LJ fluid

!				alphaij = ALP( typepj, typeci, h )
				alphaij = 12

				epsij = damp2ci * damp2pj * EPS( typepj, typeci, h )

				sigij = 0.5 * ( damp3ci + damp3pj ) * SIG( typepj, typeci, h )

				sigor2 = ( sigij * sigij ) / rij2

				sigor6 = sigor2 * sigor2 * sigor2

!				factor = alphaij * epsij / ( alphaij - 6.0 ) * &
!						 ( alphaij / 6.0 ) ** ( 6.0 / ( alphaij -6.0 )  )
				factor = 4.0 * epsij 

				ENERGY(h) = ENERGY(h) + factor * &
!							( sigor2 ** ( alphaij / 2.0 ) - sigor6 )
							( sigor6 * sigor6 - sigor6 )
								
			end if
				  
		end do

	end do

end do

return

end	subroutine e6molecule


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

subroutine e6molecule2( Nc, Xc, Yc, Zc, TYPEc, &
						DAMP2c_1, DAMP3c_1, DAMP2c_2, DAMP3c_2, &
						Np, Xp, Yp, Zp, TYPEp, &
						DAMP2p_1, DAMP3p_1, DAMP2p_2, DAMP3p_2, &
						Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
						BoxSize, dENERGY )

implicit none

! This routine calculates the energy of Nc LJ beads interacting
! with Np LJ beads.

! Nc is the number of LJ 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)						:: DAMP2c_1, DAMP3c_1
real, dimension(Nc), intent(in)						:: DAMP2c_2, DAMP3c_2

! Np is the number of LJ 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)						:: DAMP2p_1, DAMP3p_1
real, dimension(Np), intent(in)						:: DAMP2p_2, DAMP3p_2

! Nham is the number of hamiltonians.
! Nljgrs is the number of LJ groups in the system.
! 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.
! CP is a rank 3 array containing the C_ij parameters for each hamiltonian.
! ALP is a rank 3 array containing the alpha_ij parameters for each hamiltonian.
! RMAX is a rank 3 array containing the Rmax_ij parameters for each hamiltonian.
									
integer, intent(in)									:: Nham
integer, intent(in)									:: Nljgrs
real, dimension(Nljgrs, Nljgrs, Nham), intent(in)	:: EPS, SIG, CP, ALP, RMAX

! BoxSize is the length of the simulation box.

real, intent(in)									:: BoxSize

! 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												:: rij2, rij
real												:: sigor2, sigor6
real												:: sigij, cpij
real												:: alphaij, epsij
real												:: damp2ci_1, damp3ci_1
real												:: damp2ci_2, damp3ci_2
real												:: damp2pj_1, damp3pj_1
real												:: damp2pj_2, damp3pj_2
real												:: factor

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)

	damp2ci_1 = DAMP2c_1(i)
	damp3ci_1 = DAMP3c_1(i)

	damp2ci_2 = DAMP2c_2(i)
	damp3ci_2 = DAMP3c_2(i)

	do j = 1, Np

		typepj = TYPEp(j)

		damp2pj_1 = DAMP2p_1(j)
		damp3pj_1 = DAMP3p_1(j)

		damp2pj_2 = DAMP2p_2(j)
		damp3pj_2 = DAMP3p_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

		rij2 = xij*xij + yij*yij + zij*zij

		do h = 1, Nham

			cpij = CP( typepj, typeci, h )

			if( cpij > 0.0 ) then ! Exp-6 fluid

				rij = sqrt( rij2 )

				! calculate both energies

				! Energy 1

				sigij = 0.5 * ( damp3ci_1 + damp3pj_1 ) * SIG( typepj, typeci, h )

				epsij = damp2ci_1 * damp2pj_1 * EPS( typepj, typeci, h )

				if( rij  < RMAX( typepj, typeci, h ) * sigij ) then

					ENERGY_1(h) = 1.0e300 * epsij

				else
				
					alphaij = ALP( typepj, typeci, h )

					sigor2 = ( sigij * sigij ) / rij2

					ENERGY_1(h) = ENERGY_1(h) + epsij / ( alphaij - 6.0 ) * &
						( 6.0 * exp( alphaij * ( 1.0 - rij / sigij ) ) - &
						  cpij * alphaij * sigor2 * sigor2 * sigor2 )

				end if

				! Energy 2

				sigij = 0.5 * ( damp3ci_2 + damp3pj_2 ) * SIG( typepj, typeci, h )

				epsij = damp2ci_2 * damp2pj_2 * EPS( typepj, typeci, h )

				if( rij  < RMAX( typepj, typeci, h ) * sigij ) then

					ENERGY_2(h) = 1.0e300 * epsij

				else
				
					alphaij = ALP( typepj, typeci, h )

					sigor2 = ( sigij * sigij ) / rij2

					ENERGY_2(h) = ENERGY_2(h) + epsij / ( alphaij - 6.0 ) * &
						( 6.0 * exp( alphaij * ( 1.0 - rij / sigij ) ) - &
						  cpij * alphaij * sigor2 * sigor2 * sigor2 )

				end if

			else ! LJ fluid

				! Energy 1

!				alphaij = ALP( typepj, typeci, h )
				alphaij = 12.0

				epsij = damp2ci_1 * damp2pj_1 * EPS( typepj, typeci, h )

				sigij = 0.5 * ( damp3ci_1 + damp3pj_1 ) * SIG( typepj, typeci, h )

				sigor2 = ( sigij * sigij ) / rij2

				sigor6 = sigor2 * sigor2 * sigor2

!				factor = alphaij * epsij / ( alphaij - 6.0 ) * &
!						 ( alphaij / 6.0 ) ** ( 6.0 / ( alphaij -6.0 )  )
				factor = 4.0 * epsij 

				ENERGY_1(h) = ENERGY_1(h) + factor * &
!							( sigor2 ** ( alphaij / 2.0 ) - sigor6 )
							( sigor6 * sigor6 - sigor6 )
								
				! Energy 2

!				alphaij = ALP( typepj, typeci, h )
				alphaij = 12.0

				epsij = damp2ci_2 * damp2pj_2 * EPS( typepj, typeci, h )

				sigij = 0.5 * ( damp3ci_2 + damp3pj_2 ) * SIG( typepj, typeci, h )

				sigor2 = ( sigij * sigij ) / rij2

				sigor6 = sigor2 * sigor2 * sigor2

!				factor = alphaij * epsij / ( alphaij - 6.0 ) * &
!						 ( alphaij / 6.0 ) ** ( 6.0 / ( alphaij -6.0 )  )
				factor = 4.0 * epsij 

				ENERGY_2(h) = ENERGY_2(h) + factor * &
!							( sigor2 ** ( alphaij / 2.0 ) - sigor6 )
							( sigor6 * sigor6 - sigor6 )
								
			end if
				  
		end do

	end do

end do

dENERGY = ENERGY_2 - ENERGY_1

return

end	subroutine e6molecule2