MODULE VibrationalState USE Constants USE CoordinatesAndVelocities ! Use Coordinates, Velocities, Masses, TotMass USE RotationalState ! Use !!!!!!! USE Random IMPLICIT NONE TYPE NormalMode ! Type containing all the informations about one Normal coordinate LOGICAL :: Active ! TRUE = the mode is 'Active' LOGICAL :: Interpolated ! TRUE = the potential has been interpolated ! Data concerning the Potential INTEGER :: NPoints ! Number of Points REAL, DIMENSION(:), POINTER :: Q ! Mass scaled coordinates [mass weighted au] REAL, DIMENSION(:), POINTER :: Pot ! Value of Potential in Q [au] REAL, DIMENSION(:), POINTER :: Coeff ! Spline coefficients REAL :: Qmin, PotQmin ! Q corresponding to minimum Pot and Pot(Qmin) [same units as above] LOGICAL :: Minimum ! TRUE = minimum computed ! Data concerning the eigenvalues of the Normal Mode Potential INTEGER :: NEigenvalues ! Number of eigenvalues REAL, DIMENSION(:), ALLOCATABLE :: Eigenvalues ! Energy of the i-th eigenvalue [meV] REAL, DIMENSION(:), ALLOCATABLE :: Qinner ! Inner turning points of the i-th eigenvalue REAL, DIMENSION(:), ALLOCATABLE :: Qouter ! Outer turning points of the i-th eigenvalue LOGICAL, DIMENSION(:), ALLOCATABLE :: Qturning ! TRUE = turning points computed REAL :: HarmonicFreq ! Harmonic Frequency of the Normal Mode, in eV REAL, DIMENSION(5,3) :: Eigenstate ! Cartesian representation of Normal coordinate END TYPE NormalMode TYPE( NormalMode), DIMENSION(9), SAVE :: NormalModes ! Array containing informations of all Normal Modes REAL, DIMENSION(:), POINTER :: pointQ, pointPot, pointCoeff ! Pointers used by FUNCTION: Potential, DPotential REAL :: InternalEnergy ! Vibrational energy corrected by PotQmin; ! Arrays for SelectInitialConditions INTEGER, PARAMETER :: idt = 999 REAL, DIMENSION(idt) :: Time, Q, CoeffQ, VQ, CoeffVQ INTEGER :: Nt LOGICAL :: DiscardTrajectory ! Flag that determine wheter initial conditions do not conserve normal mode energy conservation INTERFACE FUNCTION DPL ( N, Xi, Yi, YCoeff, X ) ! Interface to F77 function DPL REAL :: DPL INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN) :: Xi, Yi, YCoeff REAL :: X END FUNCTION DPL FUNCTION DPLP ( N, Xi, Yi, YCoeff, X ) ! Interface to F77 function DPLP REAL :: DPLP INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN) :: Xi, Yi, YCoeff REAL :: X END FUNCTION DPLP SUBROUTINE DSPLIN( N, Xi, Yi, YCoeff, Xa, Na, Xb, Nb, alpha, beta, b ) INTEGER, INTENT(IN) :: N, Na, Nb REAL, DIMENSION(N), INTENT(IN) :: Xi, Yi, alpha, beta, b REAL, INTENT(IN) :: Xa, Xb REAL, DIMENSION(N), INTENT(OUT) :: YCoeff END SUBROUTINE DSPLIN END INTERFACE CONTAINS ! --------------------------------------------------------------------------------------------------------------- SUBROUTINE PrepareVibrationalState( v , FinalInternalEnergyPerMode ) IMPLICIT NONE INTEGER :: i1, i2, i3, n1, n2 INTEGER, DIMENSION(9), INTENT(IN) :: v REAL, PARAMETER :: Precis = 1.E-10, dQ = 1. ! Parameters for searching for minima REAL :: Q, VQ, QCart, VQCart REAL, DIMENSION(9) :: PotentialEnergyIn, KineticEnergyIn ! Arrays with values of Energy before composition REAL, DIMENSION(9) :: PotentialEnergyFin, KineticEnergyFin ! Arrays with values of Energy after composition REAL, DIMENSION(9), INTENT(OUT) :: FinalInternalEnergyPerMode LOGICAL, SAVE :: FirstLine = .FALSE. IF ( .NOT. FirstLine) THEN PRINT '(1/,32("*"),A,32("*")1/)',' V I B R A T I O N A L S T A T E ' FirstLine = .TRUE. ENDIF 100 CONTINUE ! If trajectory is discarded, repeat the procedure until accepted CALL InitializeMoleculeInCOMCoordinates DO i1 = 1, 9 IF ( .NOT. NormalModes(i1)%Interpolated ) THEN ! Interpolate Potential i1 ( if not done before ) CALL InterpolatePotential( i1 ) NormalModes(i1)%Interpolated = .TRUE. PRINT*, 'Potential for Normal Coordinate ',i1,' has been interpolated. ' ENDIF CALL SelectMode( i1 ) ! Select i1 Normal Coordinate - needed by functions Potential and DPotential IF ( .NOT. NormalModes(i1)%Minimum ) THEN ! Find minimum in Potential i1 ( if not done before ) CALL SAKEN( NormalModes(i1)%Q(1), NormalModes(i1)%Q( NormalModes(i1)%NPoints ), & ! Inf. and Sup. limits dQ, Precis, DPotential, NormalModes(i1)%Qmin, i2, 6) IF( i2 /= 0 ) THEN PRINT*, 'ERROR in finding minimum in ',i1,' Potential. Stop!' STOP ENDIF NormalModes(i1)%PotQmin = Potential( NormalModes(i1)%Qmin ) NormalModes(i1)%Minimum = .TRUE. PRINT*, 'Minimum found in Normal Coordinate Potential ',i1,' at: ',NormalModes(i1)%Qmin,& ' with value (au): ',NormalModes(i1)%PotQmin ENDIF IF ( NormalModes(i1)%Active ) THEN !--------- Generate initial conditions for Active Normal Coordinates IF ( .NOT. NormalModes(i1)%Qturning( v( i1 ) ) ) THEN ! Find turning points of the corresponding Vibrational Energy ( if not done before ) ! PRINT*, 'Quanta of Vibration in Normal Coordinate ',i1,' is : ', v( i1 ) ! PRINT*, 'Vibrational Energy in Normal Coordinate ',i1,' is (eV): ',NormalModes(i1)%Eigenvalues( v( i1 ) ) / 1000. InternalEnergy = NormalModes(i1)%Eigenvalues( v( i1 ) ) / au2eV / 1000. + NormalModes(i1)%PotQmin ! Looking for outer turning point CALL SAKEN( NormalModes(i1)%Qmin, NormalModes(i1)%Q( NormalModes(i1)%NPoints ), & ! Inf. and Sup. limits dQ, Precis, PotentialMinusInternalEnergy, NormalModes(i1)%Qouter(v(i1)), i2, 6) IF( i2 /= 0 ) THEN PRINT*, 'ERROR in finding OUTER t.p. in ',i1,' Potential. Stop!' STOP ENDIF ! Looking for inner turning point CALL SAKEN( NormalModes(i1)%Q(1), NormalModes(i1)%Qmin, & ! Inf. and Sup. limits dQ, Precis, PotentialMinusInternalEnergy, NormalModes(i1)%Qinner(v(i1)), i2, 6) IF( i2 /= 0 ) THEN PRINT*, 'ERROR in finding INNER t.p. in ',i1,' Potential. Stop!' STOP ENDIF ! PRINT*, 'Outer turning point in Normal Coordinate ',i1,' is: ',NormalModes(i1)%Qouter(v(i1)) ! PRINT*, 'Inner turning point in Normal Coordinate ',i1,' is: ',NormalModes(i1)%Qinner(v(i1)) NormalModes(i1)%Qturning( v( i1 ) ) = .TRUE. ENDIF CALL SelectInitialConditions( NormalModes(i1)%Qouter(v(i1)), NormalModes(i1)%Qmin, Precis, & Q, VQ ) ! Given starting point for integrating equations of motion, give back random Q and VQ ELSE !----------- Non-active coordinates initialized in equilibrium position Q = NormalModes(i1)%Qmin VQ = 0. ENDIF !---------------------------------------------------------------------- ! PRINT*, 'Q selected in Normal Coordinate ',i1,' is: ', Q ! PRINT*, 'VQ selected in Normal Coordinate ',i1,' is: ', VQ ! Evaluate Potential and Kinetic Energy in the i1 Normal Coordinate ( in eV ) PotentialEnergyIn( i1) = ( Potential( Q ) - NormalModes(i1)%PotQmin ) * au2eV KineticEnergyIn( i1) = 0.5 * VQ**2. * au2eV ! PRINT*, 'Initial Potential Energy in Normal Coordinate ',i1,' is: (eV) ', PotentialEnergyIn( i1) ! PRINT*, 'Initial Kinetic Energy in Normal Coordinate ',i1,' is: (eV) ', KineticEnergyIn( i1) ! PRINT*, 'Initial Total Energy in Normal Coordinate ',i1,' is: (eV) ', KineticEnergyIn( i1)+PotentialEnergyIn( i1) DO i2 = 1, 5 DO i3 = 1, 3 ! Evaluate Cartesian displacement and velocities corresponding to Q and VQ QCart = Q * NormalModes(i1)%Eigenstate( i2, i3) * ( au2Ang / SQRT( Masses( i2 ) * amu2au )) VQCart = VQ * NormalModes(i1)%Eigenstate( i2, i3) * ( au2Ang * fs2au / SQRT( Masses( i2 ) * amu2au )) ! Sum displacements and velocities to the Coordinates and Velocities Coordinates( i2, i3) = Coordinates( i2, i3) + QCart Velocities( i2, i3) = Velocities( i2, i3) + VQCart ENDDO ENDDO CALL DeselectMode ! Deselect Normal Coordinate ENDDO ! Eliminate center of mass motion or spurius angular momentum if present CALL Corrections ! Evaluate the energy drift in each normal coordinate due to corrections CALL CheckVibrationalState( PotentialEnergyIn, KineticEnergyIn, PotentialEnergyFin, KineticEnergyFin ) ! If total energy drift in at least one normal coordinate > 1 meV -----> discard initial condition IF (DiscardTrajectory) GOTO 100 FinalInternalEnergyPerMode = PotentialEnergyFin + KineticEnergyFin PRINT*, '' RETURN END SUBROUTINE PrepareVibrationalState ! -------------------------------------------------------------------------------------------------------------- SUBROUTINE InterpolatePotential( i ) ! Interpolate the potential points for the i-th Normal Coordinates IMPLICIT NONE INTEGER, INTENT(IN) :: i REAL, DIMENSION(:), ALLOCATABLE :: alpha, beta, b ! working arrays for DSPLIN of length equal to NPoints ! Allocate working arrays ALLOCATE( alpha( NormalModes(i)%NPoints ), beta( NormalModes(i)%NPoints ), & b( NormalModes(i)%NPoints ) ) ! Interpolate i-th Normal Mode Points CALL DSPLIN( NormalModes(i)%NPoints, NormalModes(i)%Q, NormalModes(i)%Pot, NormalModes(i)%Coeff, & 0.0, 0, 0.0, 0, alpha, beta, b ) ! Deallocate working arrays DEALLOCATE( alpha, beta, b ) RETURN END SUBROUTINE InterpolatePotential ! -------------------------------------------------------------------------------------------------------------- SUBROUTINE SelectInitialConditions( Qstart, Qeq, Precis, Qsel, VQsel ) IMPLICIT NONE INTEGER :: i REAL, INTENT(IN) :: Qstart, Qeq, Precis REAL, INTENT(OUT) :: Qsel, VQsel REAL :: dVQin, dVQfin ! Initial and final derivatives of VQ REAL :: Period, Delta REAL :: InitialEnergy, FinalEnergy, Drift REAL, DIMENSION(2) :: z, s ! REAL, DIMENSION(:), ALLOCATABLE :: a1, a2, a3, a4, a5 ! Working variables LOGICAL :: fin ! INTEGER :: kod ! REAL :: t ! INTEGER :: nc ! Number of time in which Pot. min. is passed INTEGER :: Ntr ! Number of timestep for completing a Period REAL :: ht, hin = 5. ! Time-step for integration i = SIZE( Time ) ! Ok, not very smart solution for arrays Time, Q, CoeffQ, VQ, CoeffVQ ALLOCATE( a1(i), a2(i), a3(i), a4(i), a5(i)) ! they are large enough for all the timesteps nc = 0 Nt = 1 ! First step Time( Nt ) = 0. ! set initial time Q( Nt ) = Qstart ! set initial position VQ( Nt ) = 0. ! set initial velocity t = 0. ! z( 1 ) = Q( Nt ) ! set working variables z( 2 ) = VQ( Nt ) ! ! Compute initial Energy ( for later Check ). NB: this energy is not corrected by VQmin (no need)! InitialEnergy = 0.5 * z( 2 )**2 + Potential( z( 1 ) ) ht = hin 11 CONTINUE s = 0. IF( ht .GT. hin) ht = hin ! Always start with the same guess for integration step Nt = Nt + 1 IF(nt.GT.idt) THEN write(6,*) 'IDT too small in integration' STOP END IF CALL DIFSY2(2, t, z, ht, Precis, s, IntegrationStep, fin) ! PRINT *, t, z(1), z(2) Time( Nt ) = t Q( Nt ) = z( 1 ) VQ( Nt ) = z( 2 ) IF(Nt.EQ.2) GO TO 11 ! IF( (Q(Nt) - Qeq) * (Q(Nt-1) - Qeq) < 0. ) nc = nc + 1 ! Potential minimum passed IF( nc == 3 ) GO TO 12 ! Stop after passing three times potential minimum IF( (Q(Nt) - Qeq) > 0 .AND. VQ(Nt)*VQ(Nt-1) < 0. ) THEN ! Time interval in which the first period ends Ntr = Nt ! PRINT*, 'Number of step for completing one period: ', Ntr END IF GO TO 11 12 CONTINUE ! Check conservation of Energy FinalEnergy = 0.5 * z( 2 )**2 + Potential( z( 1 ) ) Drift = ABS( FinalEnergy - InitialEnergy ) ! PRINT*, 'Energy drift in Integration (eV): ', Drift * au2eV ! Calculate Period dVQin = -1.*DPotential( Qstart ) ! Set initial and final derivatives of VQ dVQfin = -1.*DPotential( Q( Nt ) ) ! for interpolate CALL DSPLIN( Nt, Time, VQ, CoeffVQ, dVQin, 2, dVQfin, 2, a1, a2, a3 ) ! Interpolate VQ( t ) CALL SAKEN( Time( Ntr-1 ), Time( Ntr ), 1., Precis, Velocity, Period, kod, 6 ) ! The Period is the zero of the function VQ ( t ) ! Interpolate Q(t) and VQ(t) over one Period using TRISPL (spline for periodic functions) Nt = Ntr Time( Nt ) = Period Q( Nt ) = Qstart VQ( Nt ) = 0. CALL TRISPL( Nt, Time, Q, CoeffQ, a1, a2, a3, a4, a5 ) CALL TRISPL( Nt, Time, VQ, CoeffVQ, a1, a2, a3, a4, a5 ) ! Sample uniformly period for choosing initial conditions ! RAN0 gives random number in range [0;1] Delta = Period * RAN0() ! initial conditions are Q( Delta ) and VQ( Delta ) Qsel = DPL( Nt, Time, Q, CoeffQ, Delta ) VQsel = DPL( Nt, Time, VQ, CoeffVQ, Delta ) RETURN END SUBROUTINE SelectInitialConditions ! ------------------------------------------------------------------------------------------------------------------------ SUBROUTINE SelectMode( i ) ! Direct pointers to i-th Potential data IMPLICIT NONE INTEGER, INTENT(IN) :: i pointQ => NormalModes(i)%Q pointPot => NormalModes(i)%Pot pointCoeff => NormalModes(i)%Coeff RETURN END SUBROUTINE SelectMode ! ------------------------------------------------------------------------------------------------------------------------ SUBROUTINE DeselectMode ! Nullify pointers IMPLICIT NONE NULLIFY( pointQ, pointPot, pointCoeff ) RETURN END SUBROUTINE DeselectMode ! ------------------------------------------------------------------------------------------------------------------------ FUNCTION Potential( Q ) RESULT( Pot ) ! Evaluate interpolated Potential in Q IMPLICIT NONE REAL, INTENT(IN) :: Q REAL :: Pot Pot = DPL( SIZE(pointQ), pointQ, pointPot, pointCoeff, Q ) RETURN END FUNCTION Potential ! ------------------------------------------------------------------------------------------------------------------------ FUNCTION DPotential( Q ) RESULT( dPot ) ! Evaluate derivative of interpolated Potential in Q IMPLICIT NONE REAL, INTENT(IN) :: Q REAL :: dPot dPot = DPLP( SIZE(pointQ), pointQ, pointPot, pointCoeff, Q ) RETURN END FUNCTION DPotential ! ------------------------------------------------------------------------------------------------------------------------ FUNCTION PotentialMinusInternalEnergy( Q ) RESULT( Pot ) ! This function is equal to Potential - InternalEnergy ; the zeros of this function are the turning points IMPLICIT NONE REAL, INTENT(IN) :: Q REAL :: Pot Pot = DPL( SIZE(pointQ), pointQ, pointPot, pointCoeff, Q ) - InternalEnergy RETURN END FUNCTION PotentialMinusInternalEnergy ! ------------------------------------------------------------------------------------------------------------------------- SUBROUTINE IntegrationStep( t, z, dz ) IMPLICIT NONE REAL, DIMENSION(2), INTENT(IN) :: z REAL, DIMENSION(2), INTENT(OUT) :: dz REAL, INTENT(IN) :: t dz( 1 ) = z( 2 ) ! dQ = VQ( t ) dz( 2 ) = -1. * DPotential( z(1) ) ! dVQ = -dV/dQ( Q(t) ) RETURN END SUBROUTINE IntegrationStep ! ------------------------------------------------------------------------------------------------------------------------- FUNCTION Velocity( t ) RESULT( VelQ ) ! Evaluate interpolated velocity at time t IMPLICIT NONE REAL, INTENT(IN) :: t REAL :: VelQ VelQ = DPL( NT, Time, VQ, CoeffVQ, t ) RETURN END FUNCTION Velocity ! ------------------------------------------------------------------------------------------------------------------------- SUBROUTINE Corrections IMPLICIT NONE INTEGER :: i, i1, i2,i3 REAL, DIMENSION(3) :: COM, VCOM, AngularMom, AngularVelocity REAL, DIMENSION(3,3) :: Inertia, InvInertia REAL :: TranslationalEnergy, RotationalEnergy ! 1- Find out whether the COM has been moved and if there is COM velocity. In case, correct them. CALL FindCOM( Coordinates, COM ) CALL FindVCOM( Velocities, VCOM ) DO i1 = 1,5 DO i2 = 1,3 Coordinates( i1, i2) = Coordinates( i1, i2) - COM( i2 ) Velocities( i1, i2) = Velocities( i1, i2) - VCOM( i2 ) ENDDO ENDDO ! PRINT *, 'Center of Mass coordinates: ',COM(1), COM(2), COM(3) ! PRINT *, 'Velocity of Center of Mass : ',VCOM(1), VCOM(2), VCOM(3) TranslationalEnergy = 0.5 * TotMass * ( VCOM(1)**2. + VCOM(2)**2. + VCOM(3)**2. ) & * amu2au / (fs2au*au2Ang)**2. * au2eV ! PRINT*, 'Translational Energy corrected is (eV) : ', TranslationalEnergy ! 2- Find out whether rotational velocity has been introduced. In case, correct it. ! 2.1 - Compute spurius angular momentum CALL ComputeAngularMomentum( 5, Masses, Coordinates, Velocities, AngularMom) ! PRINT*, 'Spurius AngularMom: ', AngularMom ! 2.2 - Compute Inertia tensor CALL ComputeInertiaTensor( Coordinates, Inertia) CALL Invert( Inertia, InvInertia) ! 2.3 - Now the angular velocity AngularVelocity = MATMUL( InvInertia , AngularMom ) ! 2.4 - Correct cartesian velocities with the rotational velocities DO i = 1, 5 Velocities(i,1) = Velocities(i,1) -(AngularVelocity( 2)*Coordinates( i, 3) -AngularVelocity( 3)*Coordinates( i, 2)) * fs2au Velocities(i,2) = Velocities(i,2) -(AngularVelocity( 3)*Coordinates( i, 1) -AngularVelocity( 1)*Coordinates( i, 3)) * fs2au Velocities(i,3) = Velocities(i,3) -(AngularVelocity( 1)*Coordinates( i, 2) -AngularVelocity( 2)*Coordinates( i, 1)) * fs2au ENDDO ! RotationalEnergy = 0.5 * DOT_PRODUCT( AngularVelocity, AngularMom) * au2eV ! PRINT*, 'Rotational Energy corrected is (eV) : ', RotationalEnergy ! ! 2.5 - Recompute Angular momentum for check ! AngularMom = 0. ! DO i=1,5 ! AngularMom( 1) = AngularMom( 1) + Masses( i)*( Coordinates( i, 2)*Velocities( i, 3) - Coordinates( i, 3)*Velocities( i, 2) ) ! AngularMom( 2) = AngularMom( 2) + Masses( i)*( Coordinates( i, 3)*Velocities( i, 1) - Coordinates( i, 1)*Velocities( i, 3) ) ! AngularMom( 3) = AngularMom( 3) + Masses( i)*( Coordinates( i, 1)*Velocities( i, 2) - Coordinates( i, 2)*Velocities( i, 1) ) ! ENDDO ! PRINT*, AngularMom* amu2au / au2Ang**2. / fs2au, 'AngularMom after corrections' RETURN END SUBROUTINE Corrections ! ----------------------------------------------------------------------------------------------------------------------------- SUBROUTINE CheckVibrationalState( PotentialEnergyIn, KineticEnergyIn, PotentialEnergyFin, KineticEnergyFin ) IMPLICIT NONE INTEGER :: i1, i2, i3 REAL, DIMENSION(9), INTENT(IN) :: PotentialEnergyIn, KineticEnergyIn REAL, DIMENSION(9), INTENT(OUT):: PotentialEnergyFin, KineticEnergyFin REAL :: Q,VQ, TotalDrift DiscardTrajectory = .FALSE. DO i1 = 1, 9 CALL SelectMode( i1 ) ! Select i1 Normal Coordinate to be used in functions Potential and DPotential ! Project Coordinates and Velocities on the Normal Coordinates Q's and VQ's Q = 0. VQ = 0. DO i2 = 1, 5 DO i3 = 1, 3 Q = Q + NormalModes(i1)%Eigenstate( i2, i3) * ( Coordinates( i2, i3) - EquilibriumCOMCoordinates( i2, i3)) & / au2Ang * SQRT( Masses( i2 ) * amu2au ) VQ = VQ + NormalModes(i1)%Eigenstate( i2, i3) * Velocities( i2, i3) & / ( au2Ang * fs2au ) * SQRT( Masses( i2 ) * amu2au ) ENDDO ENDDO ! Evaluate Potential and Kinetic energy after corrections PotentialEnergyFin( i1 ) = (Potential( Q ) - NormalModes(i1)%PotQmin ) * au2eV KineticEnergyFin( i1 ) = 0.5 * VQ**2. * au2eV ! Compute the total energy drift after correction TotalDrift = PotentialEnergyIn( i1 ) + KineticEnergyIn( i1 ) - (PotentialEnergyFin( i1 ) + KineticEnergyFin( i1 )) ! PRINT*, 'Drift in Potential Energy in Normal Mode ',i1,' eqaul to (eV): ', & ! PotentialEnergyIn( i1 ) - PotentialEnergyFin( i1 ) ! PRINT*, 'Drift in Kinetic Energy in Normal Mode ',i1,' eqaul to (eV): ', & ! KineticEnergyIn( i1 ) - KineticEnergyFin( i1 ) ! PRINT*, 'Drift in Total Energy in Normal Mode ',i1,' eqaul to (eV): ',TotalDrift ! If one of the Normal coordinate energy has drifted by more than 1 meV ----> discard trajectory IF (abs(TotalDrift) > 1.E-3 ) DiscardTrajectory = .TRUE. CALL DeselectMode ! Deselect coordinate ENDDO ! IF (DiscardTrajectory) PRINT*,' TRAJECTORY DISCARDED!' RETURN END SUBROUTINE CheckVibrationalState END MODULE VibrationalState