PROGRAM Main USE Input USE Random USE CoordinatesAndVelocities USE VibrationalState USE RotationalState USE VelocityDistribution USE BoltzmannDistribution ! ------------------------------------------------------------------------------------------------- IMPLICIT NONE CHARACTER(*), PARAMETER :: NameOutput = 'InitialConditionsCH4.dat' INTEGER :: Ntraj, Nframe, NSurf INTEGER :: i, iAtom, iX, iNormalMode LOGICAL :: PrintTitleResults = .FALSE. INTEGER, DIMENSION(9) :: InitialState REAL, DIMENSION(3) :: VCOM REAL :: XCart, YCart, Phi, Theta, Psi REAL :: Veltrans,RotationalEnergy, TranslationalEnergy REAL :: AverageRotationalEnergy=0., AverageTranslationalEnergy=0., AverageInternalEnergy=0., & AverageSqRotationalEnergy=0., AverageSqTranslationalEnergy=0., AverageSqInternalEnergy=0. REAL :: MaxTranslationalEnergy, MaxRotationalEnergy, MinTranslationalEnergy, MinRotationalEnergy, & MaxInternalEnergy, MinInternalEnergy REAL :: SdTranslationalEnergy, SdRotationalEnergy, SdInternalEnergy REAL, DIMENSION(9) :: FinalInternalEnergyPerMode REAL, DIMENSION(9) :: AverageInternalEnergyPerMode=0. , AverageSqInternalEnergyPerMode=0., SdInternalEnergyPerMode REAL, DIMENSION(9) :: MinInternalEnergyPerMode, MaxInternalEnergyPerMode ! ------------------------------------------------------------------------------------------------- ! Read input files CALL ReadInput ! Read seed from Input CALL ReadSeedForRandomNumberGeneration ! If velocity distribution setup the parameters for that IF ( UseVelocityDistribution ) THEN CALL SetupVelocityDistribution( TotMass*amu2kg , Vs, DelVs, Etransmax ) ENDIF ! If Boltzmann distribution, compute fractional populations of energy levels IF ( UseBoltzmannDistribution ) THEN CALL ComputeBoltzmannFactors( Tn, Evibmax ) ENDIF ! Initialize some Max e Min values for final checks DO iNormalMode =1, 9 MinInternalEnergyPerMode(iNormalMode) = NormalModes(iNormalMode)%Eigenvalues( NormalModes(iNormalMode)%NEigenvalues ) * 1000. MaxInternalEnergyPerMode(iNormalMode) = 0. ENDDO MaxTranslationalEnergy = 0. MinTranslationalEnergy = Etransmax MaxRotationalEnergy = 0. MinRotationalEnergy = 1. MaxInternalEnergy = 0. MinInternalEnergy = SUM( MinInternalEnergyPerMode ) ! Open file where initial conditions will be written OPEN(10,FILE=NameOutput,STATUS='new') WRITE(10,*) TRIM(ADJUSTL(Title)) WRITE(10,*) ' Seed = ',Seed WRITE(10,*) Numbtrajs, NtrajMin, NtrajMax DO Ntraj = NtrajMin, NtrajMax ! --------------------------------Loop around initial conditions--- ! Assign vibrational state according to distribution ( if not according to Boltz. Distrib, as from input ) IF ( UseBoltzmannDistribution ) THEN CALL SelectVibrationalState( InitialState ) v = InitialState ENDIF CALL PrepareVibrationalState( v, FinalInternalEnergyPerMode ) AverageInternalEnergy = AverageInternalEnergy + SUM( FinalInternalEnergyPerMode ) / REAL( Numbtrajs ) AverageSqInternalEnergy = AverageSqInternalEnergy + SUM( FinalInternalEnergyPerMode )**2 / REAL( Numbtrajs ) IF (SUM( FinalInternalEnergyPerMode ) > MaxInternalEnergy ) MaxInternalEnergy = SUM( FinalInternalEnergyPerMode ) IF (SUM( FinalInternalEnergyPerMode ) < MinInternalEnergy ) MinInternalEnergy = SUM( FinalInternalEnergyPerMode ) AverageInternalEnergyPerMode = AverageInternalEnergyPerMode + FinalInternalEnergyPerMode / REAL( Numbtrajs ) AverageSqInternalEnergyPerMode = AverageSqInternalEnergyPerMode + FinalInternalEnergyPerMode **2 / REAL( Numbtrajs ) DO iNormalMode =1, 9 IF( FinalInternalEnergyPerMode( iNormalMode ) > MaxInternalEnergyPerMode( iNormalMode )) & MaxInternalEnergyPerMode( iNormalMode ) = FinalInternalEnergyPerMode( iNormalMode ) IF( FinalInternalEnergyPerMode( iNormalMode ) < MinInternalEnergyPerMode( iNormalMode )) & MinInternalEnergyPerMode( iNormalMode ) = FinalInternalEnergyPerMode( iNormalMode ) ENDDO ! CALL PrepareRotationalState( J, M, K) ! Assign orientation to the molecule (if not random, as from input) IF ( RandomMolecularOrientation ) THEN Phi = 2. * PI * RAN0() ! Psi = 2. * PI * RAN0() ! Euler angles Theta = ACOS( 2. * RAN0() - 1. ) ! ELSE Phi = Phiin * PI / 180. Psi = Psiin * PI / 180. Theta = Thetain * PI / 180. ENDIF DO iAtom=1,5 ! Rotate the molecule Coordinates(iAtom, : ) = Rotation_ZXZ( Coordinates(iAtom, : ), Phi, Theta, Psi ) Velocities(iAtom, : ) = Rotation_ZXZ( Velocities(iAtom, : ), Phi, Theta, Psi ) ENDDO ! Assign impact site to the molecule (if not random, as from input) IF ( RandomImpactSite ) THEN Xin = RAN0() Yin = RAN0() ENDIF XCart = SurfaceLatticeConstant * ( Xin + Yin * COS( PI * 2./3. )) YCart = SurfaceLatticeConstant * Yin * SIN( PI * 2./3. ) DO iAtom =1 ,5 Coordinates(iAtom, 1 ) = Coordinates(iAtom, 1 ) + XCart Coordinates(iAtom, 2 ) = Coordinates(iAtom, 2 ) + YCart ! Put the COM of mass of the molecule at a distance InitialZ from the surface Coordinates(iAtom, 3 ) = Coordinates(iAtom, 3 ) + InitialZ ENDDO ! Assign translational energy to the molecule (if not according to Vel Distrib, as from input) IF ( UseVelocityDistribution ) THEN CALL SelectEnergyFromDistribution( Etransmax, Etrans) ENDIF ! Convert translational energy in COM velocity VelTrans = SQRT( 2. * Etrans / J2eV / ( TotMass * amu2kg) ) * 1.E-5 ! Normal velocity directed towards Z negative axis DO iAtom=1,5 Velocities(iAtom,3) = Velocities(iAtom,3) - VelTrans ENDDO ! Check rotational energy and evaluate averages CALL CheckRotationalEnergy( 5, Masses, Coordinates, Velocities, RotationalEnergy ) AverageRotationalEnergy = AverageRotationalEnergy + RotationalEnergy / REAL( Numbtrajs ) AverageSqRotationalEnergy = AverageSqRotationalEnergy + RotationalEnergy**2. / REAL( Numbtrajs ) IF (RotationalEnergy > MaxRotationalEnergy ) MaxRotationalEnergy = RotationalEnergy IF (RotationalEnergy < MinRotationalEnergy ) MinRotationalEnergy = RotationalEnergy ! Check translational energy and evaluate averages CALL FindVCOM( Velocities, VCOM) TranslationalEnergy = 0.5 * TotMass * amu2kg * SUM(( VCOM * 1.E5 )**2.) * J2eV AverageTranslationalEnergy = AverageTranslationalEnergy + TranslationalEnergy / REAL( Numbtrajs ) AverageSqTranslationalEnergy = AverageSqTranslationalEnergy + TranslationalEnergy**2. / REAL( Numbtrajs ) IF (TranslationalEnergy > MaxTranslationalEnergy ) MaxTranslationalEnergy = TranslationalEnergy IF (TranslationalEnergy < MinTranslationalEnergy ) MinTranslationalEnergy = TranslationalEnergy ! Print information about the Ntraj-th molecule in stout IF (.NOT. PrintTitleResults) THEN PRINT '(1/,41("*"),A,42("*")1/)',' R E S U L T S ' PRINT '(1/,4x,a,4x,a,1x,a,3x,a,4x,a,2x,a,4x,2(3x,a,3x),1/)', 'ntr', & 'v1 v2 v3 v4_1 v4_2 v5_1 v5_2 v6_1 v6_2','etrans(eV)','phi','theta','psi','X','Y' PrintTitleResults = .TRUE. ENDIF PRINT '(1x,i6,4x,9(i1,4x),3x,f7.4,1x,3(f6.2,1x),2(f7.4,1x))', Ntraj, ( v( i), i=1,9 ), TranslationalEnergy, & Phi*180./PI, Theta*180./PI, Psi*180./PI, XCart, YCart ! Write Initial conditon for the Ntraj-th molecule in Output file WRITE(10,*) Ntraj DO iAtom=1,5 WRITE(10,'(3(1x,f12.9))') (Coordinates(iAtom,iX),iX=1,3) ENDDO DO iAtom=1,5 WRITE(10,'(3(1x,f13.9))') (Velocities(iAtom,iX),iX=1,3) ENDDO ! Generating two random number in order to sample surface initial conditions Nsurf = IRAN ( Nsurfa, Nsurfb ) Nframe = IRAN ( Nframea, Nframeb ) WRITE(10,*) NSurf, Nframe ENDDO ! ----------------------------------------------------------------------------------------- CLOSE(10) PRINT '(1/,40("*"),A,41("*")1/)',' A V E R A G E S ' SdInternalEnergy = SQRT( AverageSqInternalEnergy - AverageInternalEnergy**2. ) SdInternalEnergyPerMode = SQRT( AverageSqInternalEnergyPerMode - AverageInternalEnergyPerMode**2. ) SdRotationalEnergy = SQRT( AverageSqRotationalEnergy - AverageRotationalEnergy**2. ) SdTranslationalEnergy = SQRT( AverageSqTranslationalEnergy - AverageTranslationalEnergy**2. ) DO i = 1,9 PRINT '(" For Normal Coordinate : ",i2,1x," (eV) = ",f13.11," sd = ",f13.11," E min = ",f13.11," E max = ",f13.11,1/)',& i,AverageInternalEnergyPerMode ( i ), SdInternalEnergyPerMode(i), MinInternalEnergyPerMode( i ), MaxInternalEnergyPerMode( i ) ENDDO PRINT'(1/," Average Vibrational Energy (eV) = ",f13.11," sd = ",f13.11," E min = ",f13.11," E max = ",f13.11)', & AverageInternalEnergy, SdInternalEnergy, MinInternalEnergy, MaxInternalEnergy PRINT'(1/," Average Translational Energy (eV) = ",f13.11," sd = ",f13.11," E min = ",f13.11," E max = ",f13.11)', & AverageTranslationalEnergy, SdTranslationalEnergy, MinTranslationalEnergy, MaxTranslationalEnergy PRINT'(1/," Average Rotational Energy (eV) = ",f13.11," sd = ",f13.11," E min = ",f13.11," E max = ",f13.11,1/)', & AverageRotationalEnergy, SdRotationalEnergy, MinRotationalEnergy, MaxRotationalEnergy PRINT '(1/,A,100("*"),1/)',' ' CALL WriteSeedForNextSeries STOP END PROGRAM