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 IMPLICIT NONE CHARACTER(*), PARAMETER :: NameOutput = 'InitialConditionsHOD.dat' INTEGER :: Ntraj, Nframe, NSurf INTEGER :: i, iAtom, iX, iNormalMode LOGICAL :: PrintTitleResults = .FALSE. INTEGER, DIMENSION(3) :: InitialState REAL, DIMENSION(3) :: VCOM REAL :: XCart, YCart, Phi, Theta, Psi, Etrans_traj REAL, DIMENSION(3,3) :: RotationalMatrix REAL :: Veltrans,RotationalEnergy, TranslationalEnergy, FinalMagnitudeJ, & FinalProjJ_M, FinalProjJ_Ka, FinalProjJ_Kc REAL :: AverageRotationalEnergy=0., AverageTranslationalEnergy=0., & AverageSqRotationalEnergy=0., AverageSqTranslationalEnergy=0. REAL :: AverageMagnitudeJ=0., AverageProjJ_M=0., AverageProjJ_Ka=0., AverageProjJ_Kc=0., & AverageSqMagnitudeJ=0., AverageSqProjJ_M=0., AverageSqProjJ_Ka=0., AverageSqProjJ_Kc=0. REAL :: MaxTranslationalEnergy, MaxRotationalEnergy, MinTranslationalEnergy, MinRotationalEnergy REAL :: SdTranslationalEnergy=0., SdRotationalEnergy=0., SdMagnitudeJ=0., & SdProjJ_M=0., SdProjJ_Ka=0., SdProjJ_Kc=0. REAL, DIMENSION(3) :: FinalInternalEnergyPerMode REAL, DIMENSION(3) :: AverageInternalEnergyPerMode=0. , AverageSqInternalEnergyPerMode=0., SdInternalEnergyPerMode REAL, DIMENSION(3) :: 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, 3 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 !!! OPEN(2612,FILE='water.xyz',STATUS='unknown') 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, 3 IF( FinalInternalEnergyPerMode( iNormalMode ) > MaxInternalEnergyPerMode( iNormalMode )) & MaxInternalEnergyPerMode( iNormalMode ) = FinalInternalEnergyPerMode( iNormalMode ) IF( FinalInternalEnergyPerMode( iNormalMode ) < MinInternalEnergyPerMode( iNormalMode )) & MinInternalEnergyPerMode( iNormalMode ) = FinalInternalEnergyPerMode( iNormalMode ) ENDDO !!! WRITE(2612,*) 6 WRITE(2612,*) " " ! 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,3 ! 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, Ka, Kc CALL PrepareRotationalState( J, M, Ka, Kc, 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. ) DO iAtom =1 ,3 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 ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional IF ( ShiftVelocityDistribution ) THEN ! 30-05-2016 ! Etrans = Etrans + EtransShift Etrans_traj = Etrans + EtransShift ELSE Etrans_traj = Etrans ENDIF ! Convert translational energy in COM velocity ! VelTrans = SQRT( 2. * Etrans / J2eV / ( TotMass * amu2kg) ) * 1.E-5 VelTrans = SQRT( 2. * Etrans_traj / J2eV / ( TotMass * amu2kg) ) * 1.E-5 ! Normal velocity directed towards Z negative axis DO iAtom=1,3 Velocities(iAtom,3) = Velocities(iAtom,3) - VelTrans ENDDO ! Check rotational energy and evaluate averages CALL CheckRotationalEnergy( 3, Masses, Coordinates, Velocities, RotationalEnergy, FinalMagnitudeJ, & FinalProjJ_M, FinalProjJ_Ka, FinalProjJ_Kc ) 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_Ka = AverageProjJ_Ka + FinalProjJ_Ka / REAL( Numbtrajs ) AverageSqProjJ_Ka = AverageSqProjJ_Ka + FinalProjJ_Ka**2. / REAL( Numbtrajs ) AverageProjJ_Kc = AverageProjJ_Kc + FinalProjJ_Kc / REAL( Numbtrajs ) AverageSqProjJ_Kc = AverageSqProjJ_Kc + FinalProjJ_Kc**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','etrans(eV)','phi','theta','psi','X','Y' PrintTitleResults = .TRUE. ENDIF PRINT '(1x,i6,4x,3(i1,4x),3x,f7.4,1x,3(f6.2,1x),2(f7.4,1x))', Ntraj, ( v( i), i=1,3 ), 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,3 WRITE(10,'(3(1x,f12.9))') (Coordinates(iAtom,iX),iX=1,3) ENDDO DO iAtom=1,3 WRITE(10,'(3(1x,f13.9))') (Velocities(iAtom,iX),iX=1,3) ENDDO !!! WRITE(2612,'(a,3(1x,f12.9))') 'O',(Coordinates(1,iX),iX=1,3) WRITE(2612,'(a,3(1x,f12.9))') 'H',(Coordinates(2,iX),iX=1,3) WRITE(2612,'(a,3(1x,f12.9))') 'H',(Coordinates(3,iX),iX=1,3) ! 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) !!! CLOSE(2612) 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_Ka = SQRT( AverageSqProjJ_Ka - AverageProjJ_Ka**2. ) SdProjJ_Kc = SQRT( AverageSqProjJ_Kc - AverageProjJ_Kc**2. ) DO i = 1,3 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)', & AverageRotationalEnergy, SdRotationalEnergy, MinRotationalEnergy, MaxRotationalEnergy PRINT'(" Average |J| = ",f10.6," sd = ",f10.6,", M = ",f10.6," sd = ",f10.6)', & AverageMagnitudeJ, SdMagnitudeJ, AverageProjJ_M, SdProjJ_M PRINT'(" Average Ka = ",f10.6," sd = ",f10.6,", Kc = ",f10.6," sd = ",f10.6,1/)', & AverageProjJ_Ka, SdProjJ_Ka, AverageProjJ_Kc, SdProjJ_Kc PRINT '(1/,A,100("*"),1/)',' ' CALL WriteSeedForNextSeries STOP END PROGRAM