PROGRAM Main USE Input USE Random USE CoordinatesAndVelocities USE VibrationalState USE RotationalState USE VelocityDistribution USE BoltzmannDistribution ! ------------------------------------------------------------------------------------------------- ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional ! 23-12-2016 Bug fix: allow for energy shift with monocromatic beams ! 23-12-2016 Implement general surface geometries IMPLICIT NONE CHARACTER(*), PARAMETER :: NameOutput = 'InitialConditionsCHD3.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, DIMENSION(3,3) :: RotationalMatrix REAL :: Veltrans,RotationalEnergy, TranslationalEnergy, FinalMagnitudeJ, FinalProjJ_M, FinalProjJ_K REAL :: AverageRotationalEnergy=0., AverageTranslationalEnergy=0., & AverageSqRotationalEnergy=0., AverageSqTranslationalEnergy=0. REAL :: AverageMagnitudeJ=0., AverageProjJ_M=0., AverageProjJ_K=0., & AverageSqMagnitudeJ=0., AverageSqProjJ_M=0., AverageSqProjJ_K=0. REAL :: MaxTranslationalEnergy, MaxRotationalEnergy, MinTranslationalEnergy, MinRotationalEnergy REAL :: SdTranslationalEnergy=0., SdRotationalEnergy=0., SdMagnitudeJ=0., SdProjJ_M=0., SdProjJ_K=0. 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 ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional IF (ShiftVelocityDistribution) THEN PRINT'(1/,A)'," WARNING: we are applying a Translational Energy shift to all the molecules!" PRINT'(" The Translational Energy is selected from the Velocity Distribution; then E = ",f13.11," is added. ")', EtransShift 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 = NormalModes(iNormalMode)%Eigenvalues( NormalModes(iNormalMode)%NEigenvalues ) * 1000. MaxInternalEnergyPerMode = 0. ENDDO MaxTranslationalEnergy = 0. MinTranslationalEnergy = Etransmax MaxRotationalEnergy = 0. MinRotationalEnergy = 1. ! 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 ) 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 ! Assign orientation to the molecule IF ( J == 0 ) THEN ! If J = 0, choose between isotropic distribution and specified one 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 RotationalMatrix = RotationalMatrix_ZYZ( Phi, Theta, Psi) 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 ) Coordinates( iAtom, :) = MATMUL( RotationalMatrix, Coordinates( iAtom, :)) Velocities( iAtom, :) = MATMUL( RotationalMatrix, Velocities( iAtom, :)) ENDDO ELSE ! If J /= 0, align the molecule according to M and K CALL PrepareRotationalState( J, M, K, RotationalMatrix) ! from the rotational matrix we can compute the Phi,Theta,Psi values corresponding to a ZYZ rotation Theta = ACOS( RotationalMatrix(3,3) ) ! Phi = ASIN( RotationalMatrix(3,1) / SIN(Theta) ) Phi = ACOS( RotationalMatrix(3,1) / SIN(Theta) ) ! Psi = ASIN( RotationalMatrix(1,3) / SIN(Theta) ) Psi = ASIN( RotationalMatrix(2,3) / SIN(Theta) ) ENDIF ! PRINT*, " " ! PRINT*, " Rotational matrix:" ! PRINT*, RotationalMatrix(1,1), RotationalMatrix(1,2), RotationalMatrix(1,3) ! PRINT*, RotationalMatrix(2,1), RotationalMatrix(2,2), RotationalMatrix(2,3) ! PRINT*, RotationalMatrix(3,1), RotationalMatrix(3,2), RotationalMatrix(3,3) ! PRINT*, " " ! 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. ) ! 23-12-2016 Implement general surface geometries XCart = A1 * Xin + A2 * Yin * COS( AlphaLattice ) YCart = A2 * Yin * SIN( AlphaLattice ) 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, Etrans0) ENDIF ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional IF ( ShiftVelocityDistribution ) THEN Etrans = Etrans0 + EtransShift ! 23-12-2016 Bug fix: allow for energy shift with monocromatic beams ELSE Etrans = Etrans0 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, FinalMagnitudeJ, FinalProjJ_M, FinalProjJ_K ) AverageRotationalEnergy = AverageRotationalEnergy + RotationalEnergy / REAL( Numbtrajs ) AverageSqRotationalEnergy = AverageSqRotationalEnergy + RotationalEnergy**2. / REAL( Numbtrajs ) IF (RotationalEnergy > MaxRotationalEnergy ) MaxRotationalEnergy = RotationalEnergy IF (RotationalEnergy < MinRotationalEnergy ) MinRotationalEnergy = RotationalEnergy AverageMagnitudeJ = AverageMagnitudeJ + FinalMagnitudeJ / REAL( Numbtrajs ) AverageSqMagnitudeJ = AverageSqMagnitudeJ + FinalMagnitudeJ**2. / REAL( Numbtrajs ) AverageProjJ_M = AverageProjJ_M + FinalProjJ_M / REAL( Numbtrajs ) AverageSqProjJ_M = AverageSqProjJ_M + FinalProjJ_M**2. / REAL( Numbtrajs ) AverageProjJ_K = AverageProjJ_K + FinalProjJ_K / REAL( Numbtrajs ) AverageSqProjJ_K = AverageSqProjJ_K + FinalProjJ_K**2. / REAL( Numbtrajs ) ! 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' 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 ' SdInternalEnergyPerMode = SQRT( AverageSqInternalEnergyPerMode - AverageInternalEnergyPerMode**2. ) SdRotationalEnergy = SQRT( AverageSqRotationalEnergy - AverageRotationalEnergy**2. ) SdTranslationalEnergy = SQRT( AverageSqTranslationalEnergy - AverageTranslationalEnergy**2. ) SdMagnitudeJ = SQRT( AverageSqMagnitudeJ - AverageMagnitudeJ**2. ) SdProjJ_M = SQRT( AverageSqProjJ_M - AverageProjJ_M**2. ) SdProjJ_K = SQRT( AverageSqProjJ_K - AverageProjJ_K**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 Translational Energy (eV) = ",f13.11," sd = ",f13.11," E min = ",f13.11," E max = ",f13.11)', & AverageTranslationalEnergy, SdTranslationalEnergy, MinTranslationalEnergy, MaxTranslationalEnergy ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional IF (ShiftVelocityDistribution) THEN PRINT'(1/,A)'," WARNING: we are applying a Translational Energy shift to all the molecules!" PRINT'(" The Translational Energy is selected from the Velocity Distribution; then E = ",f13.11," is added. ")', EtransShift ENDIF 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/," Average |J| = ",f8.6," sd = ",f9.6,", M = ",f9.6," sd = ",f9.6," K = ",f9.6," sd = ",f9.6,1/)', & AverageMagnitudeJ, SdMagnitudeJ, AverageProjJ_M, SdProjJ_M, AverageProjJ_K, SdProjJ_K PRINT '(1/,A,100("*"),1/)',' ' CALL WriteSeedForNextSeries STOP END PROGRAM