C    *******************************************************************        
C    ** FICHE F.15                                                    **        
C    ** THREE METHODS FOR UNIFORM RANDOM ROTATION OF LINEAR MOLECULE. **        
C    **                                                               **        
C    ** ROUTINES SUPPLIED:                                            **        
C    **                                                               **        
C    ** SUBROUTINE FLIP1 ( EXIOLD, EYIOLD, EZIOLD, DPHIMX, DCOSMX,    **        
C    **                    EXINEW, EYINEW, EZINEW         )           **        
C    ** SUBROUTINE FLIP2 ( EXIOLD, EYIOLD, EZIOLD, DGAMAX,            **        
C    **                    EXINEW, EYINEW, EZINEW         )           **        
C    ** SUBROUTINE FLIP3 ( EXIOLD, EYIOLD, EZIOLD, DOTMIN,            **        
C    **                    EXINEW, EYINEW, EZINEW         )           **        
C    **                                                               **        
C    ** ROUTINE REQUIRED:                                             **        
C    **                                                               **        
C    ** REAL FUNCTION RANF ( DUMMY )                                  **        
C    **    RETURNS A UNIFORM RANDOM VARIATE ON THE RANGE ZERO TO ONE  **        
C    **    INCLUSIVE (SEE F.11).                                      **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** REAL      EXIOLD,EYIOLD,EZIOLD  OLD AXIAL VECTOR FOR MOL. I   **        
C    ** REAL      EYINEW,EYINEW,EZINEW  TRIAL AXIAL VECTOR FOR MOL. I **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
        SUBROUTINE FLIP1 ( EXIOLD, EYIOLD, EZIOLD, DPHIMX, DCOSMX,              
     :                     EXINEW, EYINEW, EZINEW         )                     
                                                                                
C    *******************************************************************        
C    ** MAKES A RANDOM CHANGE IN THE POLAR ANGLES PHI AND THETA.      **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** REAL      DCOSMX  MAXIMUM CHANGE IN COS(THETA)                **        
C    ** REAL      DPHIMX  MAXIMUM CHANGE IN PHI                       **        
C    ** REAL      PHIOLD  PHI IN THE OLD STATE                        **        
C    ** REAL      PHINEW  PHI IN THE NEW TRIAL STATE                  **        
C    ** REAL      COSOLD  COS(THETA) IN THE OLD STATE                 **        
C    ** REAL      COSNEW  COS(THETA) IN THE NEW TRIAL STATE           **        
C    **                                                               **        
C    ** USAGE:                                                        **        
C    **                                                               **        
C    ** FLIP1 MAKES A RANDOM CHANGE IN PHI AND COS(THETA).            **        
C    ** THE MAXIMUM ALLOWED CHANGES IN THESE VARIABLES ARE CONTROLLED **        
C    ** BY THE PARAMETERS DPHIMX AND DCOSMX RESPECTIVELY.             **        
C    ** PHI AND THETA ARE THE EULER ANGLES DESCRIBING THE ORIENTATION **        
C    ** OF THE AXIAL VECTOR. THIS METHOD CAN BE READILY EXTENDED TO   **        
C    ** POLYATOMICS BY CHANGING THE THIRD EULER ANGLE , PSI, IN THE   **        
C    ** RANGE ZERO TO TWOPI. IT WOULD BE FASTER TO PASS THE VARIABLES **        
C    ** PHIOLD AND COSOLD DIRECTLY TO FLIP1 IF THEY ARE AVAILABLE IN  **        
C    ** THE MAIN PROGRAM. SIMILARLY PHINEW AND COSNEW COULD BE        **        
C    ** PASSED DIRECTLY BACK THROUGH THE SUBROUTINE HEADER            **        
C    *******************************************************************        
                                                                                
        REAL        EXIOLD, EYIOLD, EZIOLD, EXINEW, EYINEW, EZINEW              
        REAL        DPHIMX, DCOSMX                                              
                                                                                
        REAL        COSNEW, COSOLD, PHINEW, PHIOLD, SINNEW                      
        REAL        TWOPI, PI                                                   
        REAL        RANF, DUMMY                                                 
                                                                                
        PARAMETER ( TWOPI = 6.2831853 )                                         
                                                                                
C    *******************************************************************        
                                                                                
C    ** CONVERT THE AXIAL VECTOR TO THE EULER ANGLES **                         
                                                                                
        COSOLD = EZIOLD                                                         
        PHIOLD = ATAN2 ( EYIOLD, EXIOLD )                                       
                                                                                
C    ** PERFORM THE DISPLACEMENTS **                                            
                                                                                
        PHINEW = PHIOLD + ( 2.0 * RANF ( DUMMY ) - 1.0 ) * DPHIMX               
        PHINEW = PHINEW - ANINT ( PHINEW / TWOPI ) * TWOPI                      
        COSNEW = COSOLD + ( 2.0 * RANF ( DUMMY ) - 1.0 ) * DCOSMX               
        COSNEW = COSNEW - ANINT ( COSNEW / 2.0 ) * 2.0                          
        SINNEW = SQRT ( 1.0 - COSNEW * COSNEW )                                 
                                                                                
C    ** CONVERT THE EULER ANGLES TO AXIAL VECTORS **                            
                                                                                
        EXINEW = COS ( PHINEW ) * SINNEW                                        
        EYINEW = SIN ( PHINEW ) * SINNEW                                        
        EZINEW = COSNEW                                                         
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE FLIP2 ( EXIOLD, EYIOLD, EZIOLD, DGAMAX,                      
     :                     EXINEW, EYINEW, EZINEW         )                     
                                                                                
C    *******************************************************************        
C    ** PERFORMS RANDOM ROTATION ABOUT SPACE-FIXED AXES.              **        
C    **                                                               **        
C    ** REFERENCE:                                                    **        
C    **                                                               **        
C    ** BARKER AND WATTS, CHEM PHYS LETTS 3, 144, 1969.               **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** REAL      DGAMAX  MAXIMUM ANGULAR DISPLACEMENT IN RADIANS     **        
C    ** INTEGER   IAXIS   SPACE FIXED AXIS FOR ROTATION               **        
C    **                   (1 = X, 2 = Y, 3 = Z)                       **        
C    **                                                               **        
C    ** USAGE:                                                        **        
C    **                                                               **        
C    ** FLIP2 CHOOSES ONE OF THE THREE SPACE-FIXED AXES AT RANDOM     **        
C    ** AND ROTATES THE MOLECULE AROUND THIS AXIS BY DGAMMA RADIANS.  **        
C    ** THE MAXIMUM ANGULAR DISPLACEMENT IS DGAMAX. THIS METHOD CAN   **        
C    ** READILY EXTENDED TO POLYATOMIC MOLECULES.                     **        
C    *******************************************************************        
                                                                                
        REAL    EXIOLD, EYIOLD, EZIOLD, EXINEW, EYINEW, EZINEW, DGAMAX          
                                                                                
        REAL    COSDG, SINDG, DGAMMA                                            
        REAL    RANF, DUMMY                                                     
        INTEGER IAXIS                                                           
                                                                                
C    *******************************************************************        
                                                                                
C    ** CHOOSE A SPACE FIXED AXIS AT RANDOM **                                  
                                                                                
        IAXIS = INT ( 3.0 * RANF ( DUMMY ) ) + 1                                
                                                                                
C    ** CHOOSE A RANDOM ROTATION **                                             
                                                                                
        DGAMMA = ( 2.0 * RANF ( DUMMY ) - 1.0 ) * DGAMAX                        
                                                                                
C    ** SET UP THE ROTATION MATRIX **                                           
                                                                                
        COSDG = COS ( DGAMMA )                                                  
        SINDG = SIN ( DGAMMA )                                                  
                                                                                
C    ** PERFORM ROTATIONS **                                                    
                                                                                
        IF ( IAXIS .EQ. 1 ) THEN                                                
                                                                                
           EXINEW = EXIOLD                                                      
           EYINEW = COSDG * EYIOLD + SINDG * EZIOLD                             
           EZINEW = COSDG * EZIOLD - SINDG * EYIOLD                             
                                                                                
        ELSE IF ( IAXIS .EQ. 2 ) THEN                                           
                                                                                
           EXINEW = COSDG * EXIOLD - SINDG * EZIOLD                             
           EYINEW = EYIOLD                                                      
           EZINEW = COSDG * EZIOLD + SINDG * EXIOLD                             
                                                                                
        ELSE                                                                    
                                                                                
           EXINEW = COSDG * EXIOLD + SINDG * EYIOLD                             
           EYINEW = COSDG * EYIOLD - SINDG * EXIOLD                             
           EZINEW = EZIOLD                                                      
                                                                                
        ENDIF                                                                   
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE FLIP3 ( EXIOLD, EYIOLD, EZIOLD, DOTMIN,                      
     :                     EXINEW, EYINEW, EZINEW         )                     
                                                                                
C    *******************************************************************        
C    ** CHOOSES A RANDOM DISPLACEMENT ON THE SURFACE OF A UNIT SPHERE.**        
C    **                                                               **        
C    ** REFERENCE:                                                    **        
C    **                                                               **        
C    ** MARSAGLIA, ANN MATHS STAT 43, 645, 1972.                      **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** REAL      DOT                   DOT PRODUCT OF OLD AND NEW    **        
C    **                                 AXIAL VECTORS                 **        
C    ** REAL      DOTMIN                PARAMETER TO ADJUST MAXIMUM   **        
C    **                                 DISPLACEMENT. DOTMIN SHOULD   **        
C    **                                 BE LESS THAN ONE              **        
C    **                                                               **        
C    ** USAGE:                                                        **        
C    **                                                               **        
C    ** FLIP3 USES A REJECTION TECHNIQUE TO CREATE A TRIAL            **        
C    ** ORIENTATION OF MOLECULE I SUBJECT TO THE CONSTRAINT THAT      **        
C    ** THE COSINE OF THE ANGLE BETWEEN THE OLD AND NEW AXIAL         **        
C    ** VECTORS IS GREATER THAN ( 1.0 - DOTMIN ).                     **        
C    *******************************************************************        
                                                                                
        REAL    EXIOLD, EYIOLD, EZIOLD, EXINEW, EYINEW, EZINEW, DOTMIN          
                                                                                
        REAL    DOT, XI1, XI2, XI, XISQ                                         
        REAL    RANF, DUMMY                                                     
                                                                                
C    *******************************************************************        
                                                                                
C    ** INITIALISE DOT **                                                       
                                                                                
        DOT  = 0.0                                                              
                                                                                
C    ** ITERATIVE LOOP **                                                       
                                                                                
1000    IF ( ( 1.0 - DOT ) .GE. DOTMIN ) THEN                                   
                                                                                
C       ** INITIALISE XISQ **                                                   
                                                                                
           XISQ = 1.0                                                           
                                                                                
C       ** INNER ITERATIVE LOOP **                                              
                                                                                
2000       IF ( XISQ .GE. 1.0 ) THEN                                            
                                                                                
              XI1  = RANF ( DUMMY ) * 2.0 - 1.0                                 
              XI2  = RANF ( DUMMY ) * 2.0 - 1.0                                 
              XISQ = XI1 * XI1 + XI2 * XI2                                      
                                                                                
              GOTO 2000                                                         
                                                                                
           ENDIF                                                                
                                                                                
           XI = SQRT ( 1.0 - XISQ )                                             
           EXINEW = 2.0 * XI1 * XI                                              
           EYINEW = 2.0 * XI2 * XI                                              
           EZINEW = 1.0 - 2.0 * XISQ                                            
           DOT    = EXINEW * EXIOLD + EYINEW * EYIOLD + EZINEW * EZIOLD         
                                                                                
           GOTO 1000                                                            
                                                                                
        ENDIF                                                                   
                                                                                
        RETURN                                                                  
        END                                                                     
