create.f90

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

		call Fourier_Move( Nion, Xion, Yion, Zion, &
						   TYPEion, DAMPion, &
						   Nham, Niongrs, CHARGE, BoxSize, &
						   Kmax, Nkvec, KX, KY, KZ, CONST, &
						   EXPX, EXPY, EXPZ, &
						   EXPX_tmp, EXPY_tmp, EXPZ_tmp, &
						   SUMQEXPV, SUMQEXPV_new, UFOURIER )

		UFOURIER = UFOURIER * CoulCombo

		EXPX = EXPX_tmp
		EXPY = EXPY_tmp
		EXPZ = EXPZ_tmp
		
		SUMQEXPV = SUMQEXPV_new

		call RealInteract( Nion, Xion, Yion, Zion, &
						   TYPEion, DAMPion, &
						   Nmol(0), LENGTHion, Nham, &
						   Niongrs, CHARGE, BoxSize, &
						   Alpha, UREAL)
						
		UREAL = UREAL * CoulCombo

		do j = 1, Nmol(0)

			sp = SPECIES(j)
			stion = STARTion(j)
			lion = LENGTHion(j)
			endion = stion + lion - 1
		
			call SelfMolecule( lion, Xion(stion:endion), &
							   Yion(stion:endion), &
							   Zion(stion:endion), &
							   TYPEion(stion:endion), &
							   DAMPion(stion:endion), &
							   Nham, Niongrs, CHARGE, BoxSize, &
							   Alpha, USELF_MOL(j,:) )

			USELF_MOL(j,:) = USELF_MOL(j,:) * CoulCombo

			do k = 1, LENGTHlj(j) + LENGTHion(j) - 1
				
				do m = k + 1, LENGTHlj(j) + LENGTHion(j)

					if( BONDSAPART(k,m,sp) >= 3 .AND. BEADTYPE(k,sp) == 'ION' .AND. &
						BEADTYPE(m,sp) == 'ION') then

						ionbk = stion + IONBEAD(k,sp) - 1
						ionbm = stion + IONBEAD(m,sp) - 1

						call IonMolecule( 1, Xion(ionbm), Yion(ionbm), &
									   	  Zion(ionbm), TYPEion(ionbm), &
										  DAMPion(ionbm), &
									 	  1, Xion(ionbk), Yion(ionbk), &
										  Zion(ionbk), TYPEion(ionbk), &
										  DAMPion(ionbk), &
										  Nham, Niongrs, CHARGE, &
										  BoxSize, UION_MOL_part )

						if( BONDSAPART(k,m,sp) == 3 ) UION_MOL_part = f14_ion * UION_MOL_part 

						UION_MOL(j,:) = UION_MOL(j,:) + UION_MOL_part(:) * CoulCombo

					end if

				end do

			end do

		end do

		do h = 1, Nham
				
			do j = 1, Nion
		
				USELF_CH(h) = USELF_CH(h) + DAMPion(j) * DAMPion(j) * &
											CHARGE( TYPEion(j), h ) * &
											CHARGE( TYPEion(j), h )
				
			end do

			USELF_CH(h) = USELF_CH(h) * &
							Alpha / sqrt(Pi) / BoxSize * CoulCombo

		end do

	end if

	do j = 1, Nmol(0)

		stlj = STARTlj(j)
		stion = STARTion(j)
		sp = SPECIES(j)

		do n = 1, NSTEPS(sp)

			do k = STEPSTART(n,sp), STEPSTART(n,sp) + STEPLENGTH(n,sp) - 1

				ljb = LJBEAD(k,sp)
				ionb = IONBEAD(k,sp)

				select case ( METHOD(k,sp) )

						case( 'ConeTor' )
					
							do m = 1, 3

								if( BEADTYPE( INTPARAM(m,k,sp), sp ) == 'LJ' )	then
						
									Xc(m) = Xlj( stlj + LJBEAD( INTPARAM(m,k,sp), sp ) - 1 )
									Yc(m) = Ylj( stlj + LJBEAD( INTPARAM(m,k,sp), sp ) - 1 )
									Zc(m) = Zlj( stlj + LJBEAD( INTPARAM(m,k,sp), sp ) - 1 )

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

									Xc(m) = Xion( stion + IONBEAD( INTPARAM(m,k,sp), sp ) - 1 )
									Yc(m) = Yion( stion + IONBEAD( INTPARAM(m,k,sp), sp ) - 1 )
									Zc(m) = Zion( stion + IONBEAD( INTPARAM(m,k,sp), sp ) - 1 )
						
								end if

							end do

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(4) = Xlj( stlj + ljb - 1 )
								Yc(4) = Ylj( stlj + ljb - 1 )
								Zc(4) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(4) = Xion( stion + ionb - 1 )
								Yc(4) = Yion( stion + ionb - 1 )
								Zc(4) = Zion( stion + ionb - 1 )

							end if

							call Outfold( 4, 0, BoxSize, Xc, Yc, Zc, &
										  ZERO2, ZERO2, ZERO2 )

							call IntraTorsion( Xc, Yc, Zc, INTPARAM(4,k,sp), &
											   REALPARAM(3:8,k,sp), Utor )

						case( 'Stretch' )
					
							if( BEADTYPE( INTPARAM(1,k,sp), sp ) == 'LJ' )	then
						
								Xc(1) = Xlj( stlj + LJBEAD( INTPARAM(1,k,sp), sp ) - 1 )
								Yc(1) = Ylj( stlj + LJBEAD( INTPARAM(1,k,sp), sp ) - 1 )
								Zc(1) = Zlj( stlj + LJBEAD( INTPARAM(1,k,sp), sp ) - 1 )

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

								Xc(1) = Xion( stion + IONBEAD( INTPARAM(1,k,sp), sp ) - 1 )
								Yc(1) = Yion( stion + IONBEAD( INTPARAM(1,k,sp), sp ) - 1 )
								Zc(1) = Zion( stion + IONBEAD( INTPARAM(1,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(2) = Xlj( stlj + ljb - 1 )
								Yc(2) = Ylj( stlj + ljb - 1 )
								Zc(2) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(2) = Xion( stion + ionb - 1 )
								Yc(2) = Yion( stion + ionb - 1 )
								Zc(2) = Zion( stion + ionb - 1 )

							end if

							call Outfold( 2, 0, BoxSize, Xc, Yc, Zc, &
										  ZERO2, ZERO2, ZERO2 )

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

							Uvib = 0.5 * REALPARAM(1,k,sp) * ( rtrial - REALPARAM(2,k,sp) ) ** 2.0


						case( 'FxBend' )

							! INTPARAM(1,k) = # of bending potentials evaluated
							! INTPARAM(2,k) = bead # from which the current bead is to be grown
						
							if( BEADTYPE( INTPARAM(2,k,sp), sp ) == 'LJ' )	then

								Xc(2) = Xlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Yc(2) = Ylj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Zc(2) = Zlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )

							else if( BEADTYPE( INTPARAM(2,k,sp), sp ) == 'ION' ) then
		
								Xc(2) = Xion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Yc(2) = Yion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Zc(2) = Zion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(3) = Xlj( stlj + ljb - 1 )
								Yc(3) = Ylj( stlj + ljb - 1 )
								Zc(3) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(3) = Xion( stion + ionb - 1 )
								Yc(3) = Yion( stion + ionb - 1 )
								Zc(3) = Zion( stion + ionb - 1 )

							end if

							do m = 1, INTPARAM(1,k,sp)
	
								if( BEADTYPE( INTPARAM(m+2,k,sp), sp ) == 'LJ' )	then
						
									Xc(m+3) = Xlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Yc(m+3) = Ylj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Zc(m+3) = Zlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m+2,k,sp), sp ) == 'ION' ) then

									Xc(m+3) = Xion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Yc(m+3) = Yion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Zc(m+3) = Zion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
						
								end if

							end do

							ntemp1 = 2+INTPARAM(1,k,sp)

							call Outfold( ntemp1, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), &
										  ZERO2, ZERO2, ZERO2 )

							do m = 1, INTPARAM(1,k,sp)

								Xc(1) = Xc(m+3)
								Yc(1) = Yc(m+3)
								Zc(1) = Zc(m+3)

								ntemp1 = 2*m
								ntemp2 = 1+2*m
		
								call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1,k,sp), &
												REALPARAM(ntemp2,k,sp), Utemp )

								Ubend = Ubend + Utemp

							end do

						case( 'FxBendTor' )

							if( BEADTYPE( INTPARAM(4,k,sp), sp ) == 'LJ' )	then
						
								Xc(3) = Xlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Yc(3) = Ylj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Zc(3) = Zlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )

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

								Xc(3) = Xion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Yc(3) = Yion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Zc(3) = Zion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(4) = Xlj( stlj + ljb - 1 )
								Yc(4) = Ylj( stlj + ljb - 1 )
								Zc(4) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(4) = Xion( stion + ionb - 1 )
								Yc(4) = Yion( stion + ionb - 1 )
								Zc(4) = Zion( stion + ionb - 1 )

							end if
					
							do m = 1, INTPARAM(1,k,sp)	+ 2 * INTPARAM(2,k,sp)

								if( BEADTYPE( INTPARAM(m+4,k,sp), sp ) == 'LJ' )	then
						
									Xc(m+4) = Xlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Yc(m+4) = Ylj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Zc(m+4) = Zlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m+4,k,sp), sp ) == 'ION' ) then

									Xc(m+4) = Xion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Yc(m+4) = Yion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Zc(m+4) = Zion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
						
								end if

							end do

							ntemp1 = 2 + INTPARAM(1,k,sp)	+ 2 * INTPARAM(2,k,sp)

							call Outfold( ntemp1, 0, BoxSize, &
										  Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 )

							do m = 1, INTPARAM(1,k,sp)

								Xc(2) = Xc(m+4)
								Yc(2) = Yc(m+4)
								Zc(2) = Zc(m+4)

								ntemp1 = 2*m
								ntemp2 = 1+2*m
		
								call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1,k,sp), &
												REALPARAM(ntemp2,k,sp), Utemp )

								Ubend = Ubend + Utemp

							end do

							nbendp = INTPARAM(1,k,sp)

							do m = 1, INTPARAM(2,k,sp)

								Xc(1) = Xc(2*m+nbendp+3)
								Yc(1) = Yc(2*m+nbendp+3)
								Zc(1) = Zc(2*m+nbendp+3)

								Xc(2) = Xc(2*m+nbendp+4)
								Yc(2) = Yc(2*m+nbendp+4)
								Zc(2) = Zc(2*m+nbendp+4)

								if( INTPARAM(3,k,sp) == 1 ) then

									ntemp1 = 2*nbendp-4+6*m
									ntemp2 = 2*nbendp+1+6*m
		
									call IntraTorsion( Xc, Yc, Zc, 1, &
													   REALPARAM(ntemp1:ntemp2,k,sp), &
													   Utemp )

								else if( INTPARAM(3,k,sp) == 2 ) then

									ntemp1 = 2*nbendp-2+4*m
									ntemp2 = 2*nbendp+1+4*m

									call IntraTorsion( Xc, Yc, Zc, 2, &
													   REALPARAM(ntemp1:ntemp2,k,sp), &
													   Utemp )

								end if

								Utor = Utor + Utemp

							end do


						case( 'StBend' )

							! INTPARAM(1,k) = # of bending potentials evaluated
							! INTPARAM(2,k) = bead # from which the current bead is to be grown
						
							if( BEADTYPE( INTPARAM(2,k,sp), sp ) == 'LJ' )	then

								Xc(2) = Xlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Yc(2) = Ylj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Zc(2) = Zlj( stlj + LJBEAD( INTPARAM(2,k,sp), sp ) - 1 )

							else if( BEADTYPE( INTPARAM(2,k,sp), sp ) == 'ION' ) then
		
								Xc(2) = Xion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Yc(2) = Yion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
								Zc(2) = Zion( stion + IONBEAD( INTPARAM(2,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(3) = Xlj( stlj + ljb - 1 )
								Yc(3) = Ylj( stlj + ljb - 1 )
								Zc(3) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(3) = Xion( stion + ionb - 1 )
								Yc(3) = Yion( stion + ionb - 1 )
								Zc(3) = Zion( stion + ionb - 1 )

							end if

							do m = 1, INTPARAM(1,k,sp)
	
								if( BEADTYPE( INTPARAM(m+2,k,sp), sp ) == 'LJ' )	then
						
									Xc(m+3) = Xlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Yc(m+3) = Ylj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Zc(m+3) = Zlj( stlj + LJBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m+2,k,sp), sp ) == 'ION' ) then

									Xc(m+3) = Xion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Yc(m+3) = Yion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
									Zc(m+3) = Zion( stion + IONBEAD( INTPARAM(m+2,k,sp), sp ) - 1 )
						
								end if

							end do

							ntemp1 = 2+INTPARAM(1,k,sp)

							call Outfold( ntemp1, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), &
										  ZERO2, ZERO2, ZERO2 )

							rtrial = sqrt( ( Xc(3) - Xc(2) ) ** 2 + ( Yc(3) - Yc(2) ) ** 2 + &
										   ( Zc(3) - Zc(2) ) ** 2 )

							Uvib = 0.5 * REALPARAM(1,k,sp) * ( rtrial - REALPARAM(2,k,sp) ) ** 2.0

							do m = 1, INTPARAM(1,k,sp)

								Xc(1) = Xc(m+3)
								Yc(1) = Yc(m+3)
								Zc(1) = Zc(m+3)

								ntemp1 = 1+2*m
								ntemp2 = 2+2*m
		
								call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1,k,sp), &
												REALPARAM(ntemp2,k,sp), Utemp )

								Ubend = Ubend + Utemp

							end do

						case( 'StBendTor' )

							if( BEADTYPE( INTPARAM(4,k,sp), sp ) == 'LJ' )	then
						
								Xc(3) = Xlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Yc(3) = Ylj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Zc(3) = Zlj( stlj + LJBEAD( INTPARAM(4,k,sp), sp ) - 1 )

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

								Xc(3) = Xion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Yc(3) = Yion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
								Zc(3) = Zion( stion + IONBEAD( INTPARAM(4,k,sp), sp ) - 1 )
						
							end if

							if( BEADTYPE(k,sp) == 'LJ' ) then
			
								Xc(4) = Xlj( stlj + ljb - 1 )
								Yc(4) = Ylj( stlj + ljb - 1 )
								Zc(4) = Zlj( stlj + ljb - 1 )
						
							else if( BEADTYPE(k,sp) == 'ION' ) then
					
								Xc(4) = Xion( stion + ionb - 1 )
								Yc(4) = Yion( stion + ionb - 1 )
								Zc(4) = Zion( stion + ionb - 1 )

							end if
					
							do m = 1, INTPARAM(1,k,sp)	+ 2 * INTPARAM(2,k,sp)

								if( BEADTYPE( INTPARAM(m+4,k,sp), sp ) == 'LJ' )	then
						
									Xc(m+4) = Xlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Yc(m+4) = Ylj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Zc(m+4) = Zlj( stlj + LJBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )

								else if( BEADTYPE( INTPARAM(m+4,k,sp), sp ) == 'ION' ) then

									Xc(m+4) = Xion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Yc(m+4) = Yion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
									Zc(m+4) = Zion( stion + IONBEAD( INTPARAM(m+4,k,sp), sp ) - 1 )
						
								end if

							end do

							ntemp1 = 2 + INTPARAM(1,k,sp)	+ 2 * INTPARAM(2,k,sp)

							call Outfold( ntemp1, 0, BoxSize, &
										  Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 )

							rtrial = sqrt( ( Xc(4) - Xc(3) ) ** 2 + ( Yc(4) - Yc(3) ) ** 2 + &
										   ( Zc(4) - Zc(3) ) ** 2 )

							Uvib = 0.5 * REALPARAM(1,k,sp) * ( rtrial - REALPARAM(2,k,sp) ) ** 2.0

							do m = 1, INTPARAM(1,k,sp)

								Xc(2) = Xc(m+4)
								Yc(2) = Yc(m+4)
								Zc(2) = Zc(m+4)

								ntemp1 = 1+2*m
								ntemp2 = 2+2*m
		
								call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1,k,sp), &
												REALPARAM(ntemp2,k,sp), Utemp )

								Ubend = Ubend + Utemp

							end do

							nbendp = INTPARAM(1,k,sp)

							do m = 1, INTPARAM(2,k,sp)

								Xc(1) = Xc(2*m+nbendp+3)
								Yc(1) = Yc(2*m+nbendp+3)
								Zc(1) = Zc(2*m+nbendp+3)

								Xc(2) = Xc(2*m+nbendp+4)
								Yc(2) = Yc(2*m+nbendp+4)
								Zc(2) = Zc(2*m+nbendp+4)

								if( INTPARAM(3,k,sp) == 1 ) then

									ntemp1 = 2*nbendp-3+6*m
									ntemp2 = 2*nbendp+2+6*m
		
									call IntraTorsion( Xc, Yc, Zc, 1, &
													   REALPARAM(ntemp1:ntemp2,k,sp), &
													   Utemp )

								else if( INTPARAM(3,k,sp) == 2 ) then

									ntemp1 = 2*nbendp-1+4*m
									ntemp2 = 2*nbendp+2+4*m

									call IntraTorsion( Xc, Yc, Zc, 2, &
													   REALPARAM(ntemp1:ntemp2,k,sp), &
													   Utemp )

								end if

								Utor = Utor + Utemp

							end do

						case default

							Uvib = 0.0
							Ubend = 0.0
							Utor = 0.0

					end select