resread.f90

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

subroutine ResRead( InFile, MaxBeads, reslen, &
					Xb, Yb, Zb, &
					nresmv, maxresmv, maxresp, &
					PROB_MV, RESMOVE, RESPARAM, &
					nvib, nbend, ntor, &
					maxvib, maxbend, maxtor, &
					IVIB, IBEND, ITOR, &
					RVIB, RBEND, RTOR, &
					BEADTYPE, GROUPTYPE, &
					BONDSAPART, Nham, Nljgrs, &
					EPS, SIG, CP, ALP, RMAX, &
					Niongrs, CHARGE, &
					f14_lj, f14_ion, &
					BoxSize, Uintra, Ulj, Uion )

implicit none

character*30, intent(in)				:: InFile

integer									:: MaxBeads
integer									:: reslen

real, dimension(MaxBeads)				:: Xb, Yb, Zb

integer									:: nresmv, maxresmv, maxresp

real, dimension(maxresmv)				:: PROB_MV

character*10, dimension(maxresmv)		:: RESMOVE

integer, dimension(maxresp, maxresmv)	:: RESPARAM

integer									:: nvib, nbend, ntor
integer, intent(in)						:: maxvib, maxbend, maxtor

integer, dimension(2,maxvib)			:: IVIB
integer, dimension(3,maxbend)			:: IBEND
integer, dimension(4,maxtor)			:: ITOR

real, dimension(2,maxvib)				:: RVIB
real, dimension(2,maxbend)				:: RBEND
real, dimension(4,maxtor)				:: RTOR

character*5, dimension(MaxBeads)		:: BEADTYPE
integer, dimension(MaxBeads)			:: GROUPTYPE 
integer, dimension(MaxBeads,MaxBeads)	:: BONDSAPART

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

real, dimension(Nljgrs,Nljgrs,Nham)		:: EPS
real, dimension(Nljgrs,Nljgrs,Nham)		:: SIG
real, dimension(Nljgrs,Nljgrs,Nham)		:: CP
real, dimension(Nljgrs,Nljgrs,Nham)		:: ALP
real, dimension(Nljgrs,Nljgrs,Nham)		:: RMAX

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

real, intent(in)						:: f14_lj
real, intent(in)						:: f14_ion

! BoxSize is the length of the simulation box.

real, intent(in)						:: BoxSize

real									:: Uintra

real, dimension(Nham)					:: Ulj, Uion

! Local stuff

integer							:: i, j, bd

real, parameter					:: Pi = 3.14159265359
real, parameter					:: ec = 1.60217733e-19
real, parameter					:: eps0 = 8.854187817e-12
real, parameter					:: kB = 1.380658e-23

real							:: r, utemp
real							:: uvib, ubend, utor
real							:: CoulCombo
real, dimension(1)				:: ZERO2 = 0.0
real, dimension(1)				:: ONE = 1.0
real, dimension(4)				:: Xc, Yc, Zc
real, dimension(Nham)			:: U_part


CoulCombo = ec * ec * 1.0e10 / ( 4.0 * Pi * eps0 * kB )

open(77, file = InFile)

read(77,*)
read(77,*) reslen
read(77,*)
read(77,*) nvib, nbend, ntor

read(77,*)

do i = 1, nvib

	read(77,*) bd, IVIB(1:2,i), RVIB(1:2,i)

end do

read(77,*)

do i = 1, nbend

	read(77,*) bd, IBEND(1:3,i), RBEND(1:2,i)

end do

RBEND(2,:) = RBEND(2,:) * Pi / 180.0

read(77,*)

do i = 1, ntor

	read(77,*) bd, ITOR(1:4,i), RTOR(1:4,i)

end do

read(77,*)
read(77,*) nresmv

read(77,*)

do i = 1, nresmv

	read(77,*) PROB_MV(i), RESMOVE(i)
	
	select case ( RESMOVE(i) )

		case( 'BeadDisp' )

		case( 'CrankSh' )

			read(77,*) RESPARAM(1:4,i)
			read(77,*) RESPARAM(5:4+RESPARAM(1,i),i)

		case( 'BendSh' )

			read(77,*) RESPARAM(1:3,i)
			read(77,*) RESPARAM(4:3+RESPARAM(1,i),i)

		case( 'ConRot4' )

			read(77,*) RESPARAM(1:7,i)

		case( 'ConRot7' )

		case( 'Pass' )

	end select
	
end do

do i = 2, nresmv

	PROB_MV(i) = PROB_MV(i) + PROB_MV(i-1)

end do

PROB_MV = PROB_MV / PROB_MV(nresmv)

read(77,*)	
read(77,*)	
read(77,*)	
read(77,*)	
read(77,*)	

do i = 1, reslen

	read(77,*) bd, Xb(bd), Yb(bd), Zb(bd)

end do

close(77)

call Outfold( reslen, 0, BoxSize, &
			  Xb, Yb, Zb, ZERO2, ZERO2, ZERO2 )

uvib = 0.0
ubend = 0.0
utor = 0.0

do i = 1, nvib

	do j = 1, 2

		Xc(j) = Xb( IVIB(j,i) )
		Yc(j) = Yb( IVIB(j,i) )
		Zc(j) = Zb( IVIB(j,i) )

	end do

	r = sqrt( ( Xc(2) - Xc(1) ) ** 2 + ( Yc(2) - Yc(1) ) ** 2 + &
			  ( Zc(2) - Zc(1) ) ** 2 )

	uvib = uvib + 0.5 * RVIB(1,i) * ( r - RVIB(2,i) ) ** 2.0

end do

do i = 1, nbend

	do j = 1, 3

		Xc(j) = Xb( IBEND(j,i) )
		Yc(j) = Yb( IBEND(j,i) )
		Zc(j) = Zb( IBEND(j,i) )

	end do

	call IntraBend( Xc, Yc, Zc, RBEND(1,i), &
					RBEND(2,i), utemp )

	ubend = ubend + utemp

end do

do i = 1, ntor

	do j = 1, 4

		Xc(j) = Xb( ITOR(j,i) )
		Yc(j) = Yb( ITOR(j,i) )
		Zc(j) = Zb( ITOR(j,i) )

	end do

	call IntraTorsion( Xc, Yc, Zc, 2, &
					   RTOR(1:4,i), utemp )

	utor = utor + utemp

end do

uintra = uvib + ubend + utor

Ulj = 0.0

! If the intramolecular nonbonded interactions are included
! as part of the reservoir energy, then uncomment the 
! following lines.  Also change the appropriate lines in 
! cranksh.f90, alpha_ch.f90 and grow.f90. 

!do i = 1, reslen - 1
!
!	do j = i, reslen
!
!		if( BEADTYPE(i) == 'LJ' .AND. BEADTYPE(j) == 'LJ' .AND. &
!			BONDSAPART(i,j) >= 3 ) then
!
!			call e6molecule( 1, Xb(i), Yb(i), &
!							 Zb(i), GROUPTYPE(i), &
!							 ONE, ONE, &
!							 1, Xb(j), Yb(j), &
!							 Zb(j), GROUPTYPE(j), &
!							 ONE, ONE, &
!							 Nham, Nljgrs, EPS, SIG, CP, &
!							 ALP, RMAX, BoxSize, U_part )
!
!			if( BONDSAPART(i,j) == 3 ) U_part = f14_lj * U_part 
!
!			Ulj = Ulj + U_part
!
!		end if
!
!	end do
!
!end do

Uion = 0.0

!do i = 1, reslen - 1
!
!	do j = i, reslen
!
!		if( BEADTYPE(i) == 'ION' .AND. BEADTYPE(j) == 'ION' .AND. &
!			BONDSAPART(i,j) >= 3 ) then
!
!			call IonMolecule( 1, Xb(i), Yb(i), &
!							  Zb(i), GROUPTYPE(i), ONE, &
!							  1, Xb(j), Yb(j), &
!							  Zb(j), GROUPTYPE(j), ONE, &
!							  Nham, Niongrs, CHARGE, &
!							  BoxSize, U_part )
!
!			if( BONDSAPART(i,j) == 3 ) U_part = f14_ion * U_part 
!
!			Uion = Uion + U_part * CoulCombo
!
!		end if
!
!	end do
!
!end do

return

end subroutine ResRead