MODULE RotationalState USE CoordinatesAndVelocities ! Use Coordinates, Velocities, Masses, TotMass USE Constants USE Random IMPLICIT NONE ! REAL, DIMENSION(5,3) :: Coordinates=0.0, Velocities=0.0 ! REAL, DIMENSION(5,3) :: EquilibriumCOMCoordinates ! REAL, PARAMETER :: TotMass = MC + MH + 3.*MD ! REAL, DIMENSION(5) :: Masses = (/ MC, MH, MD, MD, MD/) ! INTEGER :: Seed = 1234 ! LOGICAL :: EquilibriumCOMComputed = .FALSE. REAL, DIMENSION(3), SAVE :: FigureAxis ! INTERFACE ! FUNCTION RAN0 ( N ) ! Interface to F77 function RAN0 ! REAL :: RAN0 ! INTEGER, INTENT(IN) :: N ! END FUNCTION RAN0 ! END INTERFACE CONTAINS ! ----------------------------------------- SUBROUTINE PrepareRotationalState( J, M, K, RotationalMatrix) IMPLICIT NONE INTEGER :: i,i1,i2,i3,n1,n2 INTEGER , INTENT(IN) :: J, M, K REAL, DIMENSION(3) :: AngularMom, z, AngularVelocity REAL, DIMENSION(3,3) :: Inertia, InvInertia REAL, DIMENSION(3,3) :: RotationalMatrix1, RotationalMatrix2 REAL, DIMENSION(3,3), INTENT(OUT) :: RotationalMatrix REAL, DIMENSION(5,3) :: RotationalVelocities 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) RotationalEnergyRef = REAL( J*( J + 1)) / (2.* Inertia( 1, 1)) + & REAL( K )**2. * ( 1./ (2.* Inertia( 3, 3)) - 1./ (2.*Inertia( 1, 1))) PRINT '(1/,32("*"),A,32("*")1/)',' R O T A T I O N A L S T A T E ' PRINT '(A,3(i2,A))', 'The rotational state of the molecule is: ( J =',J,', M =',M,', K =',K,' ).' PRINT '(A,f7.4,/)' , 'The Rotational energy of the symmetric top rotor in this state is (meV) : ', & RotationalEnergyRef * au2eV * 1000. FirstLine = .TRUE. ENDIF ! Determine inertia tensor for distorted geometry CALL ComputeInertiaTensor( Coordinates, Inertia) ! Setup first rotation: orient the molecule with J along z-axis, and projection ! of J on molecular symmetry axis equal to K alph = 2.*PI*RAN0() ! Random rotation of molecule around figure axis bet = ACOS( REAL( K) / SQRT( REAL( J *( J + 1 )) ) ) gamm = 2.*PI*RAN0() ! Random projection of J on (X,Y)-plane RotationalMatrix1 = RotationalMatrix_ZYZ( alph, bet, gamm ) ! Setup second rotation: orient the molecule in such a way that the projection ! of J on z-axis is equal to M phi = 0. !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 two rotations in one Rotational Matrix and print in in file RotationalMatrix = MATMUL( RotationalMatrix2, RotationalMatrix1 ) DO i = 1,4 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 figure axis for chek K (projection of L on this axis ) FigureAxis = (/ 0., 0., 1. /) FigureAxis = MATMUL( RotationalMatrix, FigureAxis ) ! Now the Angular Momentum (in atomic units) AngularMom = (/ 0., 0., SQRT( REAL( J*( J + 1 )) ) /) AngularMom = MATMUL( RotationalMatrix2, AngularMom ) ! 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, 4 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, 4 DO i2 = 1, 3 Velocities(i1,i2) = Velocities(i1,i2) + RotationalVelocities(i1,i2) ENDDO ENDDO ! Do some checks.. ! CALL CheckRotationalState ! CALL CheckRotationalEnergy( 5, Masses, Coordinates, Velocities, RotationalEnergy ) RETURN END SUBROUTINE ! ----------------------------------------- ! SUBROUTINE CheckRotationalState ! Check angular momentum ( magnitude and projections) and rotational energy from Coordinates and Velocities ! INTEGER :: i, i1, i2 ! REAL :: MagnitudeJ, ProjJ_M, ProjJ_K ! REAL :: RotationalEnergy ! REAL, DIMENSION(3) :: COM ! REAL, DIMENSION(3) :: AngularMom, AngularVelocity ! REAL, DIMENSION(3,3) :: Inertia, InvInertia ! Compute Angular Momentum (in atomic units) ! CALL ComputeAngularMomentum( 5, Masses, Coordinates, Velocities, AngularMom) ! Compute the inverse of Inertia tensor (in atomic units) ! CALL ComputeInertiaTensor( Coordinates, Inertia) ! CALL Invert( Inertia, InvInertia) ! Now the angular velocity (in atomic units) ! AngularVelocity = MATMUL( InvInertia, AngularMom) !!!!!!!!!!!!!!!! ! do i1=1,3 ! print*, InvInertia( i1, 1),InvInertia( i1, 2),InvInertia( i1, 3) ! enddo !!!!!!!!!!!!!!!! ! Compute magnitude of angular momentum,.. ! MagnitudeJ = SQRT( SUM ( AngularMom**2.) ) ! .. the projection on the molecular symmetry axis,.. ! ProjJ_K = DOT_PRODUCT( FigureAxis , AngularMom) ! / SQRT( SUM ( FigureAxis **2.)) ! Not needed, since magnitude 1. ! .. 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) ! PRINT '(A,f7.4)', 'The Rotational energy of the generated molecule is (meV) : ', RotationalEnergy * au2eV * 1000. ! PRINT '(A,f7.4)', 'The magnitude of Angular momentum is (au) :', MagnitudeJ ! PRINT '(A,f7.4)', 'The projection of Angular momentum on molecular axis is (au): ', ProjJ_K ! PRINT '(A,f7.4)', 'The projection of Angular momentum on z-axis is(au): ', ProjJ_M ! PRINT*,'' ! RETURN ! END SUBROUTINE ! ----------------------------------------- SUBROUTINE CheckRotationalEnergy( N, M, X, V, RotationalEnergy, MagnitudeJ, ProjJ_M, ProjJ_K ) IMPLICIT NONE INTEGER , INTENT(IN) :: N REAL, INTENT(OUT) :: MagnitudeJ, ProjJ_M, ProjJ_K 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( 4, 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 molecular symmetry axis,.. ProjJ_K = DOT_PRODUCT( FigureAxis , AngularMom) ! / SQRT( SUM ( FigureAxis**2.)) ! Not needed, since magnitude 1. ! .. 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 InitializeMoleculeInCOMCoordinates ! IMPLICIT NONE ! INTEGER :: i,l ! REAL, DIMENSION(3) :: COM ! ! IF(.NOT.EquilibriumCOMComputed) THEN ! CALL FindCOM( EquilibriumCoordinates, COM) ! DO i=1,5 ! DO l=1,3 ! EquilibriumCOMCoordinates(i,l) = EquilibriumCoordinates(i,l)-COM(l) ! ENDDO ! ENDDO ! EquilibriumCOMComputed = .TRUE. ! ENDIF ! ! Coordinates = EquilibriumCOMCoordinates ! ! END SUBROUTINE ! ----------------------------------------- ! SUBROUTINE FindCOM( Positions, COM) ! IMPLICIT NONE ! INTEGER :: i,l ! REAL, DIMENSION(5,3), INTENT(IN) :: Positions ! REAL, DIMENSION(3) , INTENT(OUT) :: COM ! ! COM = 0. ! DO i=1,3 ! DO l=1,5 ! COM( i) = COM( i) + Masses( l) * Positions( l, i) / TotMass ! ENDDO ! ENDDO ! ! 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(4,3), INTENT(IN) :: Positions REAL, DIMENSION(3,3), INTENT(OUT) :: Inertia Inertia= 0. DO i1 = 1,3 DO i2 = 1,3 DO i3 = 1,4 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 ! END MODULE