MODULE RotationalState USE CoordinatesAndVelocities ! Use Coordinates, Velocities, Masses, TotMass USE Constants USE Random IMPLICIT NONE REAL, DIMENSION(3), SAVE :: aAxis, cAxis REAL, DIMENSION(3,3), SAVE :: RotationalAxes CONTAINS ! ----------------------------------------- SUBROUTINE PrepareRotationalState( J, M, Ka, Kc, RotationalMatrix) IMPLICIT NONE INTEGER :: i,i1,i2,i3,n1,n2 INTEGER , INTENT(IN) :: J, M, Ka, Kc REAL, DIMENSION(3) :: AngularMom, z, AngularVelocity, InertiaMoments REAL, DIMENSION(3,3) :: Inertia, InvInertia REAL, DIMENSION(3,3) :: RotationalMatrix1, RotationalMatrix2 REAL, DIMENSION(3,3), INTENT(OUT) :: RotationalMatrix REAL, DIMENSION(3,3) :: RotationalVelocities REAL, DIMENSION(3,3) :: Rotated REAL :: phi, theta, psi REAL :: alph, bet, gamm REAL :: RotationalEnergyRef, RotationalEnergy LOGICAL , SAVE :: FirstLine = .FALSE. IF ( .NOT. FirstLine) THEN ! Compute rotational energy for oblate symmetric top rotor in state (J, M, K) (in atomic units) CALL ComputeInertiaTensor( EquilibriumCOMCoordinates, Inertia) PRINT '(1/,32("*"),A,32("*")1/)',' R O T A T I O N A L S T A T E ' PRINT '(A,4(i2,A))', 'The rotational state of the molecule is: ( J =',J,', M =',M,', Ka =',Ka,', Kc =',Kc,' ).' IF (Ka*Ka + Kc*Kc > J*(J+1)) THEN PRINT '(A)', "J projection in ac-plane is longer than J. Rotational state does not exist." STOP ENDIF CALL CalculatePrincipalRotationalAxes( Inertia, RotationalAxes, InertiaMoments) PRINT '(1/,A)', 'Rotational constants for HOD in its equilibrium geometry (cm-1):' PRINT '(3(1x,f8.5))', 1./( 2. * InertiaMoments) *au2eV *eV2wn ! Experimental values available at DOI: 10.1006/jmsp.1993.1266 (see also CCCBDB) DO i = 1,3 Rotated( i, :) = MATMUL( TRANSPOSE(RotationalAxes), EquilibriumCOMCoordinates( i, :)) ENDDO PRINT '(1/,A)', 'Equilibrium coordinates of molecule with principal axes (a,b,c) along (x,y,z) :' PRINT '(a,3(1x,f12.9))', 'O', (Rotated(1,i),i=1,3) PRINT '(a,3(1x,f12.9))', 'H', (Rotated(2,i),i=1,3) PRINT '(a,3(1x,f12.9))', 'D', (Rotated(3,i),i=1,3) FirstLine = .TRUE. ENDIF ! The transpose of RotationalAxes (RA^T) represents the FIRST rotation, bringing HOD ! to a reference frame with (a,b,c) along (x,y,z), respectively. ! Setup SECOND rotation: orient the molecule with J along z-axis, and projection ! of J on a-axis equal to Ka and projection of J on c-axis equal to Kc alph = ACOS( REAL( Ka ) / SQRT( REAL( J *( J + 1 ) - Kc*Kc) ) ) ! Projection of J on a-axis equal to Ka bet = ACOS( REAL( Kc ) / SQRT( REAL( J *( J + 1 )) ) ) ! Projection of J on c-axis equal to Kc gamm = 0. ! The third projection is automatically determined! RotationalMatrix1 = RotationalMatrix_ZYZ( alph, bet, gamm ) ! Setup THIRD rotation: orient the molecule in such a way that the projection ! of J on z-axis is equal to M phi = 2.*PI*RAN0() theta = ACOS( REAL( M) / SQRT( REAL( J *( J + 1 )) ) ) psi = 2.*PI*RAN0() ! Random orientation of figure axis around J RotationalMatrix2 = RotationalMatrix_ZYZ( phi, theta, psi ) ! Compose the three rotations in one Rotational Matrix: R = R2 * R1 * RA^T RotationalMatrix = MATMUL( RotationalMatrix2, MATMUL( RotationalMatrix1, TRANSPOSE(RotationalAxes) )) DO i = 1,3 Coordinates( i, :) = MATMUL( RotationalMatrix, Coordinates( i, :)) Velocities( i, :) = MATMUL( RotationalMatrix, Velocities( i, :)) ENDDO ! OPEN(1, FILE='ROTATIONAL_MATRIX.dat') ! DO i=1,3 ! WRITE (1,*) RotationalMatrix(i,1), RotationalMatrix(i,2), RotationalMatrix(i,3) ! ENDDO ! CLOSE(1) ! Keep track of the a- and c-axis to chek for Ka and Kc (projection of L on the axes ) ! After applying RA^T the a- and c-axis are oriented along x and z, respectively aAxis = (/ 1., 0., 0. /) ! a is along x after RA^T cAxis = (/ 0., 0., 1. /) ! c is along z after RA^T aAxis = MATMUL( RotationalMatrix2, MATMUL( RotationalMatrix1, aAxis )) ! only apply second and third rotations cAxis = MATMUL( RotationalMatrix2, MATMUL( RotationalMatrix1, cAxis )) ! Now the Angular Momentum (in atomic units) AngularMom = (/ 0., 0., SQRT( REAL( J*( J + 1 )) ) /) AngularMom = MATMUL( RotationalMatrix2, AngularMom ) ! only apply third rotation WRITE(2612,'(a,3(1x,f12.6))') 'F',(AngularMom(i),i=1,3) WRITE(2612,'(a,3(1x,f12.6))') 'Cl',(aAxis(i),i=1,3) WRITE(2612,'(a,3(1x,f12.6))') 'Cl',(cAxis(i),i=1,3) ! Molecule is oriented and angular momentum is set: now determine rotational velcocities. ! Compute Inertia tensor (in atomic units) CALL ComputeInertiaTensor( Coordinates, Inertia) ! Invert the Inertia tensor CALL Invert( Inertia, InvInertia) ! Calculate angular velocity: omega = Inertia-1 * AngularMom (in atomic units) AngularVelocity = 0. DO i1 =1, 3 DO i2 = 1, 3 AngularVelocity( i1) = AngularVelocity( i1) + InvInertia( i1, i2)*AngularMom( i2) ENDDO ENDDO ! Calculate Rotational velocity per each atom and add it to the molecular velocity (in Ang/fs) DO i = 1, 3 RotationalVelocities(i,1) = ( AngularVelocity( 2)*Coordinates( i, 3) - AngularVelocity( 3)*Coordinates( i, 2) ) * fs2au RotationalVelocities(i,2) = ( AngularVelocity( 3)*Coordinates( i, 1) - AngularVelocity( 1)*Coordinates( i, 3) ) * fs2au RotationalVelocities(i,3) = ( AngularVelocity( 1)*Coordinates( i, 2) - AngularVelocity( 2)*Coordinates( i, 1) ) * fs2au ENDDO DO i1 = 1, 3 DO i2 = 1, 3 Velocities(i1,i2) = Velocities(i1,i2) + RotationalVelocities(i1,i2) ENDDO ENDDO RETURN END SUBROUTINE ! ----------------------------------------- SUBROUTINE CheckRotationalEnergy( N, M, X, V, RotationalEnergy, MagnitudeJ, ProjJ_M, ProjJ_Ka, ProjJ_Kc ) IMPLICIT NONE INTEGER , INTENT(IN) :: N REAL, INTENT(OUT) :: MagnitudeJ, ProjJ_M, ProjJ_Ka, ProjJ_Kc REAL, DIMENSION(N) , INTENT(IN) :: M REAL, DIMENSION(N,3), INTENT(IN) :: X,V REAL, INTENT(OUT) :: RotationalEnergy REAL, DIMENSION(N,3) :: XCOM REAL, DIMENSION(3) :: AngularVelocity, AngularMom, COM REAL, DIMENSION(3,3) :: Inertia, InvInertia INTEGER :: iAtom, iX ! Move to COM frame CALL FindCOM( X, COM) XCOM = 0. DO iAtom = 1, N DO iX = 1, 3 XCOM(iAtom, iX) = X(iAtom, iX) - COM( iX ) ENDDO ENDDO ! Compute Angular Momentum (in atomic units) CALL ComputeAngularMomentum( 3, M, XCOM, V, AngularMom) ! Compute the inverse of Inertia tensor (in atomic units) CALL ComputeInertiaTensor( XCOM, Inertia) CALL Invert( Inertia, InvInertia) ! Now the angular velocity (in atomic units) AngularVelocity = MATMUL( InvInertia, AngularMom) ! Compute magnitude of angular momentum,.. MagnitudeJ = SQRT( SUM ( AngularMom**2.) ) ! .. the projection on the a-axis,.. ProjJ_Ka = DOT_PRODUCT( aAxis , AngularMom) ! / SQRT( SUM ( aAxis**2.)) ! Not needed, since magnitude 1. ! .. the projection on the c-axis,.. ProjJ_Kc = DOT_PRODUCT( cAxis , AngularMom) ! .. and the projection on the z-axis ProjJ_M = AngularMom(3) ! Now the rotational energy of the generated initial conditions: ! Erot = 1/2 (omega)T * Inertia * omega = 1/2 (omega)T * AngularMom (in atomic units) RotationalEnergy = 0.5 * DOT_PRODUCT( AngularVelocity, AngularMom) RotationalEnergy = RotationalEnergy * au2eV RETURN END SUBROUTINE ! ----------------------------------------- SUBROUTINE ComputeAngularMomentum( N, M, X, V, AngularMom) ! Compute angular momentum (au) for a system of N atoms ! Positions (X) in Ang, Velocities (V) in Ang/fs and Masses(M) in amu IMPLICIT NONE INTEGER , INTENT(IN) :: N REAL, DIMENSION(N) , INTENT(IN) :: M REAL, DIMENSION(N,3), INTENT(IN) :: X,V REAL, DIMENSION(3), INTENT(OUT) :: AngularMom INTEGER :: i AngularMom = 0. DO i=1,N AngularMom( 1) = AngularMom( 1) + M( i)*( X( i, 2)*V( i, 3) - X( i, 3)*V( i, 2) ) AngularMom( 2) = AngularMom( 2) + M( i)*( X( i, 3)*V( i, 1) - X( i, 1)*V( i, 3) ) AngularMom( 3) = AngularMom( 3) + M( i)*( X( i, 1)*V( i, 2) - X( i, 2)*V( i, 1) ) ENDDO AngularMom = AngularMom * amu2au / au2Ang**2. / fs2au RETURN END SUBROUTINE ! ----------------------------------------- SUBROUTINE ComputeInertiaTensor( Positions, Inertia) ! Compute Inertia tensor (in atomic units) from Coordinates (in Ang) and Masses (in amu) IMPLICIT NONE INTEGER :: i1, i2, i3, n1, n2 REAL, DIMENSION(3,3), INTENT(IN) :: Positions REAL, DIMENSION(3,3), INTENT(OUT) :: Inertia Inertia= 0. DO i1 = 1,3 DO i2 = 1,3 DO i3 = 1,3 IF (i1==i2) THEN ! Diagonal elements n1 = MOD( i1 , 3) + 1 ! n1 and n2 are 2,3 for element (1,1) n2 = MOD( i1 + 1, 3) + 1 ! 1,3 for element (2,2) ! 1,2 for element (3,3) Inertia(i1,i2) = Inertia( i1, i2) + Masses( i3)*( Positions( i3, n1)**2.0 + Positions( i3, n2)**2.0 ) & *amu2au / au2Ang**2.0 ELSE ! Off-diagonal elements Inertia( i1, i2) = Inertia( i1, i2) - Masses( i3)*Positions( i3, i1)*Positions( i3, i2) *amu2au / au2Ang**2.0 ENDIF ENDDO ENDDO ENDDO RETURN END SUBROUTINE ! -------------------------------------------- SUBROUTINE Invert( Matrix, InvMatrix) ! Invert 3x3 Matrix IMPLICIT NONE INTEGER :: i1, i2 REAL :: Det REAL, DIMENSION(3,3), INTENT(IN) :: Matrix REAL, DIMENSION(3,3), INTENT(OUT) :: InvMatrix Det = Matrix( 1, 1)*( Matrix( 2, 2)*Matrix( 3, 3) - Matrix( 2, 3)*Matrix( 3, 2) ) & -Matrix( 1, 2)*( Matrix( 2, 1)*Matrix( 3, 3) - Matrix( 2, 3)*Matrix( 3, 1) ) & +Matrix( 1, 3)*( Matrix( 2, 1)*Matrix( 3, 2) - Matrix( 2, 2)*Matrix( 3, 1) ) DO i1 = 1, 3 DO i2 = 1, 3 InvMatrix( i2, i1 ) = 1.0 / Det * (Matrix( MOD( i1 , 3) + 1, MOD( i2 , 3) + 1)* & Matrix( MOD( i1 + 1, 3) + 1, MOD( i2 + 1, 3) + 1)- & Matrix( MOD( i1 , 3) + 1, MOD( i2 + 1, 3) + 1)* & Matrix( MOD( i1 + 1, 3) + 1, MOD( i2 , 3) + 1) ) ENDDO ENDDO RETURN END SUBROUTINE ! -------------------------------------------- FUNCTION Rotation_ZXZ(x, phi, theta, psi) RESULT(y) ! Rotate a vector using three euler angles. ! according to definitions in Goldstein, classical mechanics (ZXZ convention): ! conterclockwise rotation about z axis by phi (xyz -> x'y'z'); ! conterclockwise rotation about x' axis by theta (x'y'z' -> x''y''z''); ! conterclockwise rotation about z'' axis by psi (x''y''z'' -> XYZ ). IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: x REAL, DIMENSION(3) :: y REAL, DIMENSION(3,3) :: rot REAL, INTENT(IN) :: phi, theta, psi y=0.0 rot(1,1) = cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi) rot(1,2) = cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi) rot(1,3) = sin(psi)*sin(theta) rot(2,1) = -1.*sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi) rot(2,2) = -1.*sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi) rot(2,3) = cos(psi)*sin(theta) rot(3,1) = sin(theta)*sin(phi) rot(3,2) = -1.*sin(theta)*cos(phi) rot(3,3) = cos(theta) y = matmul(rot,x) RETURN END FUNCTION ! -------------------------------------------- FUNCTION RotationalMatrix_ZYZ( phi, theta, psi) RESULT(rot) ! Rotate a vector using three euler angles ! according to ZYZ convention: ! conterclockwise rotation about z axis by phi (xyz -> x'y'z'); ! conterclockwise rotation about y' axis by theta (x'y'z' -> x''y''z''); ! conterclockwise rotation about z'' axis by psi (x''y''z'' -> XYZ ). IMPLICIT NONE REAL, DIMENSION(3,3) :: rot REAL, INTENT(IN) :: phi, theta, psi rot(1,1) = cos(theta)*cos(phi)*cos(psi) - sin(phi)*sin(psi) rot(1,2) = cos(phi)*sin(psi) + cos(theta)*cos(psi)*sin(phi) rot(1,3) = -1.*cos(psi)*sin(theta) rot(2,1) = -1.*sin(phi)*cos(psi) - cos(theta)*sin(psi)*cos(phi) rot(2,2) = cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi) rot(2,3) = sin(psi)*sin(theta) rot(3,1) = sin(theta)*cos(phi) rot(3,2) = sin(theta)*sin(phi) rot(3,3) = cos(theta) RETURN END FUNCTION ! ----------------------------------------------- SUBROUTINE CalculatePrincipalRotationalAxes( Inertia, RotationalAxes, InertiaMoments) REAL, DIMENSION(3,3), INTENT(IN) :: Inertia REAL, DIMENSION(3,3), INTENT(OUT) :: RotationalAxes REAL, DIMENSION(3), INTENT(OUT) :: InertiaMoments INTEGER, PARAMETER :: LWMAX = 1000 REAL, DIMENSION(LWMAX) :: Work INTEGER :: LWork=102, INFO ! Optimal LWork determided from querying DSYEV RotationalAxes = Inertia CALL DSYEV( 'V', 'U', 3, RotationalAxes, 3, InertiaMoments, Work, LWork, INFO ) IF ( INFO /= 0 ) THEN PRINT*, 'ERROR in DSYEV. Stop' STOP ENDIF END SUBROUTINE ! END MODULE