MODULE Input USE VibrationalState USE CoordinatesAndVelocities USE Constants ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional IMPLICIT NONE CHARACTER(*), PARAMETER :: PotentialsDirectory = 'Potential_v/' ! Name directory with data Potentials CHARACTER(*), DIMENSION(9), PARAMETER :: PotentialFileName = & ! Array with name files with Potentials (/ PotentialsDirectory//'Potential_v1.dat ', & PotentialsDirectory//'Potential_v2.dat ', & PotentialsDirectory//'Potential_v3.dat ', & PotentialsDirectory//'Potential_v4_1.dat', & PotentialsDirectory//'Potential_v4_2.dat', & PotentialsDirectory//'Potential_v5_1.dat', & PotentialsDirectory//'Potential_v5_2.dat', & PotentialsDirectory//'Potential_v6_1.dat', & PotentialsDirectory//'Potential_v6_2.dat' /) CHARACTER(*), PARAMETER :: EigenstatesDirectory = 'CartesianRepresentations/' ! Name directory with Eigenstates of NModes CHARACTER(*), DIMENSION(9), PARAMETER :: EigenstateFileName = & ! Array with name files with Eigenstates of NModes (/ EigenstatesDirectory//'CartesianRepresentation_v1.dat ', & EigenstatesDirectory//'CartesianRepresentation_v2.dat ', & EigenstatesDirectory//'CartesianRepresentation_v3.dat ', & EigenstatesDirectory//'CartesianRepresentation_v4_1.dat', & EigenstatesDirectory//'CartesianRepresentation_v4_2.dat', & EigenstatesDirectory//'CartesianRepresentation_v5_1.dat', & EigenstatesDirectory//'CartesianRepresentation_v5_2.dat', & EigenstatesDirectory//'CartesianRepresentation_v6_1.dat', & EigenstatesDirectory//'CartesianRepresentation_v6_2.dat' /) CHARACTER(*), PARAMETER :: NameMainInput = 'MainInput.inp' ! Main input CHARACTER(*), PARAMETER :: NameEquilibriumCoordinates = 'EquilibriumGeometryCHD3.inp' ! Equilibrium geometry CHARACTER(*), PARAMETER :: NameFileEnergyLevelsPotentials = 'EnergyLevelsPotentials.inp' ! File with eigenvalues of every Potential INTEGER, DIMENSION(9) :: v ! Quantum numbers for vibrational state INTEGER :: J, M, K ! Quantum numbers for rotational state INTEGER :: Nsurfa, Nsurfb, Nframea, Nframeb ! First and last #surf and #frame CHARACTER(50) :: Title ! Title of output file INTEGER :: Ntrajmin, Ntrajmax, NumbTrajs ! Numbers of first and last trajectory, total number of trajs LOGICAL :: UseBoltzmannDistribution, & ! Compute initial vibrational state according to Boltzmann Distribution UseVelocityDistribution, & ! Compute initial transational energy according to a specified Velocity distribution RandomImpactSite, & ! Select Impact site randomly RandomMolecularOrientation, & ! Randomly sample molecular orientation ShiftVelocityDistribution ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional REAL :: InitialZ, & ! Initial Z position of COM Tn, Evibmax, & ! Nozzle temperature and maximum vibrational energy for computing Boltzmann Distribution Vs, Delvs, Etransmax, & ! Data for velocity distribution Etrans, & ! Normal incidence energy EtransShift ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional REAL :: SurfaceLatticeConstant, Xin, Yin ! Surface lattice constant and impact site (if not random) REAL :: Phiin,Thetain,Psiin ! Euler angle for initial molecular orientation (if not random) CONTAINS ! ------------------------------------------------------------------------------------------------- SUBROUTINE ReadInput IMPLICIT NONE PRINT'(1/,43("*"),A,44("*")1/)', ' I N P U T ' CALL ReadMainInput CALL ReadEquilibriumCoordinates CALL ReadPotentialData CALL ReadEnergyLevelsPotentials CALL ReadEigenstatesPotentials END SUBROUTINE ! ------------------------------------------------------------------------------------------------- SUBROUTINE ReadMainInput ! Read main Input file IMPLICIT NONE INTEGER :: i CHARACTER(20) :: Introduction OPEN(1,file=NameMainInput,status='old') READ(1,'(a20,a50)') Introduction, Title PRINT '(A,1/,A,1/)', ' Initial conditions for CHD3:',Title READ(1,*) READ(1,*) Ntrajmin,Ntrajmax IF ( Ntrajmax < Ntrajmin ) THEN PRINT*, ' Error in ntrajmax/ntrajmin!' STOP ENDIF Numbtrajs = Ntrajmax - Ntrajmin + 1 PRINT '(A,I5,A,I5,A,I5,1/)', ' Preparing ',Numbtrajs,' initial conditions, from: ',Ntrajmin,' to ',Ntrajmax READ(1,*) InitialZ PRINT '(A,f4.2,A,1/)', ' Center of mass is ',InitialZ,' Ang from the surface.' READ(1,*) READ(1,*) UseBoltzmannDistribution, Tn, Evibmax READ(1,*) DO i = 1, 9 READ(1,*) NormalModes( i )%Active , v( i ) ENDDO IF ( UseBoltzmannDistribution ) THEN IF ( Tn == 0. ) THEN PRINT*, ' Nozzle temperature set to 0K !' STOP ENDIF PRINT *, ' Sampling initial vibrational state from Boltzmann distribution' PRINT '(a,f6.2)', ' Nozzle temperature (K): ', Tn PRINT '(a,f4.2,1/)', ' Max evib considered (eV): ',evibmax ELSE DO i= 1, 9 IF ( .NOT. NormalModes( i )%Active ) THEN PRINT '(A,I2)', ' Normal coordinate: ',i PRINT *, ' Classical Coordinate' ELSE PRINT '(A,I2)', ' Normal coordinate: ',i PRINT '(A,I2)', ' Vibrational state = ', v( i ) ENDIF ENDDO ENDIF READ(1,*) READ(1,*) UseVelocityDistribution READ(1,*) Vs, Delvs, Etransmax ! 22-06-2015 Add shift of Vel Distrib for modeling vdW corrected functional READ(1,*) ShiftVelocityDistribution, EtransShift READ(1,*) Etrans IF ( UseVelocityDistribution ) THEN PRINT'(1/,A)', ' Sampling initial Translational energy from Velocity distribution' PRINT'(2(a,f7.2))', ' Av. vel. (m/s) = ',Vs,' width of distrib (m/s)=', DelVs PRINT'(a,f4.2,1/)' , ' Max translational energy (eV): ',Etransmax ELSE PRINT '(1/,a,f9.7,1/)', ' Perp. translational energy (eV): ',Etrans ENDIF READ(1,*) READ(1,*) J, M, K READ(1,*) RandomMolecularOrientation READ(1,*) Phiin,Thetain,Psiin IF ( J == 0 ) THEN IF ( RandomMolecularOrientation ) THEN PRINT '(A,1/)', ' Random orientation of the molecules' ELSE PRINT '(a,3(f5.1,1x),a,1/)', ' Fixed euler angles: ', Phiin, Thetain, Psiin, & ' with respect to equilibrium geometry' ENDIF ELSE PRINT '(A,1/)', ' Specific rovibrational state selected - Molecules Rotationally aligned ' PRINT '(A,i3,A,i3,A,i3,1/)', ' J = ', J, ' M = ', M, ' K = ', K ENDIF READ(1,*) READ(1,*) RandomImpactSite READ(1,*) Xin READ(1,*) Yin READ(1,*) SurfaceLatticeConstant PRINT '(a,f7.4,1/)', ' Surface Lattice Constant (Ang) :', SurfaceLatticeConstant IF ( RandomImpactSite ) THEN PRINT '(a,1/)', ' Random Impact parameter' ELSE PRINT '(a,1x,2(a,f8.5),1/)', ' Fixed impact parameter; cartesian components: ', & ' x= ', SurfaceLatticeConstant* ( Xin + Yin * COS( Pi * 4./3 )), & ' y= ', Yin *SurfaceLatticeConstant * SIN( PI * 4./3. ) ENDIF READ(1,*) Nsurfa, Nsurfb, Nframea, Nframeb PRINT '(A,i5,A,i5)', ' Surface sampled between #Surf= ',Nsurfa,' and #Surf= ',Nsurfb PRINT '(A,i5,A,i5,1/)', ' Frame sampled between #Frame= ',Nframea,' and #Frame= ',NFrameb CLOSE(1) END SUBROUTINE ! ------------------------------------------------------------------------------------------------- SUBROUTINE ReadEquilibriumCoordinates ! Read equilibrium configuration for CHD3 IMPLICIT NONE INTEGER :: i,j OPEN(1,file=NameEquilibriumCoordinates,status='old') READ(1,*) (( EquilibriumCoordinates( i, j), j= 1, 3), i= 1, 5) PRINT '(1/,A)', ' Equilibrium Coordinates Read from Input: ' DO i=1,5 PRINT '(3(1x,f9.7))', EquilibriumCoordinates( i, 1),EquilibriumCoordinates( i, 2), EquilibriumCoordinates( i, 3) ENDDO CLOSE(1) RETURN END SUBROUTINE ! ------------------------------------------------------------------------------------------------- SUBROUTINE ReadPotentialData ! Read the Potential points for all the Normal Modes IMPLICIT NONE INTEGER :: i1, i2 CHARACTER(30) :: NameFile CHARACTER(1) :: a DO i1 = 1, 9 NameFile = PotentialFileName( i1 ) OPEN(1,file=NameFile,status='old') READ(1,*) a,NormalModes(i1)%NPoints ALLOCATE( NormalModes(i1)%Q( NormalModes(i1)%NPoints ), NormalModes(i1)%Pot( NormalModes(i1)%NPoints ), & NormalModes(i1)%Coeff( NormalModes(i1)%NPoints )) NormalModes(i1)%Interpolated = .FALSE. NormalModes(i1)%Minimum = .FALSE. DO i2 = 1, NormalModes(i1)%NPoints READ(1,*) NormalModes(i1)%Q( i2 ), NormalModes(i1)%Pot( i2 ) ENDDO CLOSE(1) ENDDO END SUBROUTINE ! ------------------------------------------------------------------------------------------------- SUBROUTINE ReadEnergyLevelsPotentials ! Read eigenvalues computed for each Normal mode Potential IMPLICIT NONE INTEGER :: i1, i2, i3 CHARACTER(4) :: bla OPEN(1,file=NameFileEnergyLevelsPotentials,status='old') READ(1,*) DO i1 = 1, 9 READ(1,'(a4,i1)') bla, NormalModes(i1)%NEigenvalues ALLOCATE( NormalModes(i1)%Eigenvalues( 0:NormalModes(i1)%NEigenvalues ), & NormalModes(i1)%Qinner( 0:NormalModes(i1)%NEigenvalues ), & NormalModes(i1)%Qouter( 0:NormalModes(i1)%NEigenvalues ), & NormalModes(i1)%Qturning( 0:NormalModes(i1)%NEigenvalues ) ) DO i2 = 0, NormalModes(i1)%NEigenvalues READ(1,*) i3, NormalModes(i1)%Eigenvalues( i2 ) NormalModes(i1)%Qturning( i2 ) = .FALSE. ENDDO ENDDO CLOSE(1) RETURN END SUBROUTINE ! ------------------------------------------------------------------------------------------------- SUBROUTINE ReadEigenstatesPotentials ! Read cartesian representation of each Normal Coordinates IMPLICIT NONE CHARACTER(57) :: NameFile INTEGER :: i1, i2, i3 CHARACTER(1) :: bla DO i1 = 1, 9 NameFile = EigenstateFileName( i1 ) OPEN(1,file=NameFile,status='old') READ(1,*) READ(1,*) READ(1,*) bla, NormalModes(i1)%HarmonicFreq READ(1,*) ((NormalModes(i1)%Eigenstate( i2, i3), i3 = 1, 3), i2 = 1, 5) CLOSE(1) ENDDO RETURN END SUBROUTINE END MODULE