grow.f90

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

							xtr = Xion_new(ionb)
							ytr = Yion_new(ionb)
							ztr = Zion_new(ionb)

						end if

						call Bend( lenlj, lenion, &
									Xlj_new, Ylj_new, Zlj_new, &
									Xion_new, Yion_new, Zion_new, &
									LJBEAD, IONBEAD, &
									BEADTYPE, MaxInt, MaxReal, &
									INTPARAM(:,k), REALPARAM(:,k), &
									BoxSize, Nham, BETA, LNW, &
									FxSt, newold, &
									ncount, nDoOver, &
									xtr, ytr, ztr, &
									UVIB(j,k), UBEND(j,k), Seed )  

					case( 'FxBendTor' )

						FxSt = 1
						newold = 2

						if( BEADTYPE(k) == 'LJ' ) then
			
							xtr = Xlj_new(ljb)
							ytr = Ylj_new(ljb)
							ztr = Zlj_new(ljb)
						
						else if( BEADTYPE(k) == 'ION' ) then
					
							xtr = Xion_new(ionb)
							ytr = Yion_new(ionb)
							ztr = Zion_new(ionb)

						end if

						call BendTor( lenlj, lenion, &
									  Xlj_new, Ylj_new, Zlj_new, &
									  Xion_new, Yion_new, Zion_new, &
									  LJBEAD, IONBEAD, &
									  BEADTYPE, MaxInt, MaxReal, &
									  INTPARAM(:,k), REALPARAM(:,k), &
									  BoxSize, Nham, BETA, LNW, &
									  FxSt, newold, &
									  ncount, nDoOver, &
									  xtr, ytr, ztr, &
									  UVIB(j,k), UBEND(j,k), UTOR(j,k), &
									  Seed )
									  
					case( 'StBend' )

						FxSt = 2
						newold = 2

						if( BEADTYPE(k) == 'LJ' ) then
			
							xtr = Xlj_new(ljb)
							ytr = Ylj_new(ljb)
							ztr = Zlj_new(ljb)
						
						else if( BEADTYPE(k) == 'ION' ) then
					
							xtr = Xion_new(ionb)
							ytr = Yion_new(ionb)
							ztr = Zion_new(ionb)

						end if

						call Bend( lenlj, lenion, &
									Xlj_new, Ylj_new, Zlj_new, &
									Xion_new, Yion_new, Zion_new, &
									LJBEAD, IONBEAD, &
									BEADTYPE, MaxInt, MaxReal, &
									INTPARAM(:,k), REALPARAM(:,k), &
									BoxSize, Nham, BETA, LNW, &
									FxSt, newold, &
									ncount, nDoOver, &
									xtr, ytr, ztr, &
									UVIB(j,k), UBEND(j,k), Seed )  

					case( 'StBendTor' )

						FxSt = 2
						newold = 2

						if( BEADTYPE(k) == 'LJ' ) then
			
							xtr = Xlj_new(ljb)
							ytr = Ylj_new(ljb)
							ztr = Zlj_new(ljb)
						
						else if( BEADTYPE(k) == 'ION' ) then
					
							xtr = Xion_new(ionb)
							ytr = Yion_new(ionb)
							ztr = Zion_new(ionb)

						end if

						call BendTor( lenlj, lenion, &
									  Xlj_new, Ylj_new, Zlj_new, &
									  Xion_new, Yion_new, Zion_new, &
									  LJBEAD, IONBEAD, &
									  BEADTYPE, MaxInt, MaxReal, &
									  INTPARAM(:,k), REALPARAM(:,k), &
									  BoxSize, Nham, BETA, LNW, &
									  FxSt, newold, &
									  ncount, nDoOver, &
									  xtr, ytr, ztr, &
									  UVIB(j,k), UBEND(j,k), UTOR(j,k), &
									  Seed )
									  
					case( 'ResRand' )

						n = INTPARAM(1,k)

						if( resstart == 0 ) resstart = k

						resl = reslen(n)

						do m = 1, resl

							if( BEADTYPE(m+resstart-1) == 'LJ' ) then
		
								Xres(m) = Xlj_new( LJBEAD(m+resstart-1) )
								Yres(m) = Ylj_new( LJBEAD(m+resstart-1) )
								Zres(m) = Zlj_new( LJBEAD(m+resstart-1) )

							else if( BEADTYPE(m+resstart-1) == 'ION' ) then
							
								Xres(m) = Xion_new( IONBEAD(m+resstart-1) )
								Yres(m) = Yion_new( IONBEAD(m+resstart-1) )
								Zres(m) = Zion_new( IONBEAD(m+resstart-1) )

							end if

						end do

						do m = 1, resl

							Xres_tr(j,m) = Xres(m)
							Yres_tr(j,m) = Yres(m)
							Zres_tr(j,m) = Zres(m)

						end do												

					case default

						UVIB(j,k) = 0.0
						UBEND(j,k) = 0.0
						UTOR(j,k) = 0.0

				end select

			end do

		else
						
			if( j == 1 ) then
				
				do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1

					if( BEADTYPE(k) == 'LJ' ) then
				
						NLJGROUPS( TYPElj_new(LJBEAD(k)) ) = NLJGROUPS( TYPElj_new(LJBEAD(k)) ) + 1

					end if

				end do

			end if

			k = STEPSTART(i) - 1

			cb_step_loop: do while ( k < STEPSTART(i) + STEPLENGTH(i) - 1 )

				k = k + 1
			
				ljb = LJBEAD(k)
				ionb = IONBEAD(k)

				select case ( METHOD(k) )

					case( 'Random' )

						if( BEADTYPE(k) == 'LJ' )	then
												
							Xlj_tr(j,ljb) = ran2(Seed) * BoxSize
							Ylj_tr(j,ljb) = ran2(Seed) * BoxSize
							Zlj_tr(j,ljb) = ran2(Seed) * BoxSize

						else if( BEADTYPE(k) == 'ION' ) then

							Xion_tr(j,ionb) = ran2(Seed) * BoxSize
							Yion_tr(j,ionb) = ran2(Seed) * BoxSize
							Zion_tr(j,ionb) = ran2(Seed) * BoxSize

						end if

					case( 'Sphere' )
					
						FxSt = 1
						newold = 1

						call Stretch( lenlj, lenion, &
										Xlj_new, Ylj_new, Zlj_new, &
										Xion_new, Yion_new, Zion_new, &
										LJBEAD, IONBEAD, &
										BEADTYPE, MaxInt, MaxReal, &
										INTPARAM(:,k), REALPARAM(:,k), &
										BoxSize, Nham, BETA, &
										FxSt, newold, &
										xtr, ytr, ztr, &
										UVIB(j,k), Seed )  

						if( BEADTYPE(k) == 'LJ' ) then

							Xlj_tr(j,ljb) = xtr
							Ylj_tr(j,ljb) = ytr
							Zlj_tr(j,ljb) = ztr
						
						else if( BEADTYPE(k) == 'ION' ) then
					
							Xion_tr(j,ionb)	= xtr
							Yion_tr(j,ionb)	= ytr
							Zion_tr(j,ionb)	= ztr

						end if

					case( 'Cone' )

						Nb = INTPARAM(3,k)

						allocate( TempX( Nb ) )
						allocate( TempY( Nb ) )
						allocate( TempZ( Nb ) )

						call cone2( lenlj, lenion, Nb, &
									Xlj_new, Ylj_new, Zlj_new, &
									Xion_new, Yion_new, Zion_new, &
									LJBEAD, IONBEAD, &
									BEADTYPE, MaxInt, MaxReal, &
									INTPARAM(:,k), REALPARAM(:,k), &
									BoxSize, &
									TempX, TempY, TempZ, &
									Seed )  

						if( BEADTYPE(k) == 'LJ' ) then
									 
							ljb_temp = ljb
							ionb_temp = ionb + 1

						else if( BEADTYPE(k) == 'ION' ) then

							ljb_temp = ljb + 1
							ionb_temp = ionb 

						end if
					
						do m = 1, Nb

							if( BEADTYPE( k + m - 1 ) == 'LJ' )	then

								Xlj_tr(j,ljb_temp) = TempX(m)
								Ylj_tr(j,ljb_temp) = TempY(m)
								Zlj_tr(j,ljb_temp) = TempZ(m)

								ljb_temp = ljb_temp + 1

							else if( BEADTYPE( k + m - 1 ) == 'ION'	) then

								Xion_tr(j,ionb_temp) = TempX(m)
								Yion_tr(j,ionb_temp) = TempY(m)
								Zion_tr(j,ionb_temp) = TempZ(m)

								ionb_temp = ionb_temp + 1

							end if
					
						end do

						deallocate( TempX )
						deallocate( TempY )
						deallocate( TempZ )

					case( 'FxBend' )

						FxSt = 1
						newold = 1

						call Bend( lenlj, lenion, &
									Xlj_new, Ylj_new, Zlj_new, &
									Xion_new, Yion_new, Zion_new, &
									LJBEAD, IONBEAD, &
									BEADTYPE, MaxInt, MaxReal, &
									INTPARAM(:,k), REALPARAM(:,k), &
									BoxSize, Nham, BETA, LNW, &
									FxSt, newold, &
									ncount, nDoOver, &
									xtr, ytr, ztr, &
									UVIB(j,k), UBEND(j,k), Seed )  

						if( ncount /= 0 ) then
						  
							k = STEPSTART(i) - 1

							cycle cb_step_loop
						
						else

							if( BEADTYPE(k) == 'LJ' ) then

								Xlj_tr(j,ljb) = xtr
								Ylj_tr(j,ljb) = ytr
								Zlj_tr(j,ljb) = ztr

							else if( BEADTYPE(k) == 'ION' ) then

								Xion_tr(j,ionb) = xtr
								Yion_tr(j,ionb) = ytr
								Zion_tr(j,ionb) = ztr

							end if

						end if

					case( 'FxBendTor' )

						FxSt = 1
						newold = 1

						call BendTor( lenlj, lenion, &
									  Xlj_new, Ylj_new, Zlj_new, &
									  Xion_new, Yion_new, Zion_new, &
									  LJBEAD, IONBEAD, &
									  BEADTYPE, MaxInt, MaxReal, &
									  INTPARAM(:,k), REALPARAM(:,k), &
									  BoxSize, Nham, BETA, LNW, &
									  FxSt, newold, &
									  ncount, nDoOver, &
									  xtr, ytr, ztr, &
									  UVIB(j,k), UBEND(j,k), UTOR(j,k), &
									  Seed )
									  
						if( ncount /= 0 ) then
						  
							k = STEPSTART(i) - 1

							cycle cb_step_loop
						
						else

							if( BEADTYPE(k) == 'LJ' ) then

								Xlj_tr(j,ljb) = xtr
								Ylj_tr(j,ljb) = ytr
								Zlj_tr(j,ljb) = ztr

							else if( BEADTYPE(k) == 'ION' ) then

								Xion_tr(j,ionb) = xtr
								Yion_tr(j,ionb) = ytr
								Zion_tr(j,ionb) = ztr

							end if

						end if

					case( 'Stretch' )
					
						FxSt = 2
						newold = 1

						call Stretch( lenlj, lenion, &
										Xlj_new, Ylj_new, Zlj_new, &
										Xion_new, Yion_new, Zion_new, &
										LJBEAD, IONBEAD, &
										BEADTYPE, MaxInt, MaxReal, &
										INTPARAM(:,k), REALPARAM(:,k), &
										BoxSize, Nham, BETA, &
										FxSt, newold, &
										xtr, ytr, ztr, &
										UVIB(j,k), Seed )  

						if( BEADTYPE(k) == 'LJ' ) then

							Xlj_tr(j,ljb) = xtr
							Ylj_tr(j,ljb) = ytr
							Zlj_tr(j,ljb) = ztr
						
						else if( BEADTYPE(k) == 'ION' ) then
					
							Xion_tr(j,ionb)	= xtr
							Yion_tr(j,ionb)	= ytr
							Zion_tr(j,ionb)	= ztr

						end if

					case( 'StBend' )

						FxSt = 2
						newold = 1

						call Bend( lenlj, lenion, &
									Xlj_new, Ylj_new, Zlj_new, &
									Xion_new, Yion_new, Zion_new, &
									LJBEAD, IONBEAD, &
									BEADTYPE, MaxInt, MaxReal, &
									INTPARAM(:,k), REALPARAM(:,k), &
									BoxSize, Nham, BETA, LNW, &
									FxSt, newold, &
									ncount, nDoOver, &
									xtr, ytr, ztr, &
									UVIB(j,k), UBEND(j,k), Seed )  

						if( ncount /= 0 ) then
						  
							k = STEPSTART(i) - 1

							cycle cb_step_loop
						
						else

							if( BEADTYPE(k) == 'LJ' ) then

								Xlj_tr(j,ljb) = xtr
								Ylj_tr(j,ljb) = ytr
								Zlj_tr(j,ljb) = ztr

							else if( BEADTYPE(k) == 'ION' ) then

								Xion_tr(j,ionb) = xtr
								Yion_tr(j,ionb) = ytr
								Zion_tr(j,ionb) = ztr

							end if

						end if

					case( 'StBendTor' )

						FxSt = 2
						newold = 1

						call BendTor( lenlj, lenion, &
									  Xlj_new, Ylj_new, Zlj_new, &
									  Xion_new, Yion_new, Zion_new, &
									  LJBEAD, IONBEAD, &
									  BEADTYPE, MaxInt, MaxReal, &
									  INTPARAM(:,k), REALPARAM(:,k), &
									  BoxSize, Nham, BETA, LNW, &
									  FxSt, newold, &
									  ncount, nDoOver, &
									  xtr, ytr, ztr, &
									  UVIB(j,k), UBEND(j,k), UTOR(j,k), &
									  Seed )
									  
						if( ncount /= 0 ) then
						  
							k = STEPSTART(i) - 1

							cycle cb_step_loop
						
						else

							if( BEADTYPE(k) == 'LJ' ) then

								Xlj_tr(j,ljb) = xtr
								Ylj_tr(j,ljb) = ytr
								Zlj_tr(j,ljb) = ztr

							else if( BEADTYPE(k) == 'ION' ) then

								Xion_tr(j,ionb) = xtr
								Yion_tr(j,ionb) = ytr
								Zion_tr(j,ionb) = ztr

							end if

						end if
		
					case( 'Match' )

						if( BEADTYPE( INTPARAM(1,k) ) == 'LJ' )	then
						
							xtr = Xlj_new( LJBEAD( INTPARAM(1,k) ) )
							ytr = Ylj_new( LJBEAD( INTPARAM(1,k) ) )
							ztr = Zlj_new( LJBEAD( INTPARAM(1,k) ) )

						else if( BEADTYPE( INTPARAM(1,k) ) == 'ION' ) then

							xtr = Xion_new( IONBEAD( INTPARAM(1,k) ) )
							ytr = Yion_new( IONBEAD( INTPARAM(1,k) ) )
							ztr = Zion_new( IONBEAD( INTPARAM(1,k) ) )
						
						end if

						if( BEADTYPE(k) == 'LJ' ) then

							Xlj_tr(j,ljb) = xtr
							Ylj_tr(j,ljb) = ytr
							Zlj_tr(j,ljb) = ztr
						
						else if( BEADTYPE(k) == 'ION' ) then
					
							Xion_tr(j,ionb)	= xtr
							Yion_tr(j,ionb)	= ytr
							Zion_tr(j,ionb)	= ztr

						end if
							
					case( 'Complete' )

					case( 'ResRand' )

						n = INTPARAM(1,k)

						if( j == 1 ) then
						
							if( resstart == 0 ) resstart = k
	
							resl = reslen(n)

							resmol = int( ran2(Seed) * Nresmol ) + 1

							do m = 1, resl

								Xres(m) = Xresmols(resmol,m,n)
								Yres(m) = Yresmols(resmol,m,n)
								Zres(m) = Zresmols(resmol,m,n)

							end do												
							
						end if

						if(new) UTOR(j,k) = Uint_resm(resmol,n)
						if(new) dULJ_MOL_tr(j,:) = Ulj_resm(resmol,:,n)
						if(new) dUION_MOL_tr(j,:) = Uion_resm(resmol,:,n)

						call ResRand( resl, Xres, Yres, Zres, &
										BoxSize, Seed )