      SUBROUTINE ORIENT( GEOM1, IE1, GEOM2, IE2, DISPL, ANGLES)
*
* COMPUTE ANGLE AND DISPLACEMENT RELATION BETWEEN GEOM1 AND GEOM2
*  THE RESULT WILL BE DISPL AND ANGLES WHICH WILL CONVERT
*  GEOM1 INTO GEOM2
*
* THE ASSUMPTION IS THAT GEOM1 IS NOT ROTATED
*   GEOM1( *, 1):  0.000    0.000   0.000
*   GEOM1( *, 2):  X1       0.000   0.000
*   GEOM1( *, 3):  X2       Y2      0.000
*   GEOM1( *, 4):  X3       Y3      Z3
*
* The method here is to step-by-step find the displacements and
* rotation angles to return GEOM2 back to the starting orientation
* represented by GEOM1.
* Throughout the part the rotation routines will be called in such
* a manner as to pretend the entire GEOM2 array (here used as GEOML
* a logcal copy) is a sub-group of a larger geometry.  In this way,
* the TOTROT information will not be accumulated.
* At the end, the process will be reversed and the rotation
* information will be accumulatd via TOTROT.
*
      IMPLICIT REAL( A-H, O-Z)
      INCLUDE 'SIZES'
      DIMENSION IE1( NUMATM), IE2( NUMATM)
      DIMENSION GEOM1( 3,NUMATM), GEOM2( 3,NUMATM), GEOML( 3, NUMATM)
      DIMENSION DISPL( 3), ANGLES( 3), IGROUP( NUMATM)
      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
*
      IF ( DEBUG) THEN
        WRITE (*,*) ' ORIENT:   geom1                         geom2'
        DO 9 I=1, 5
           WRITE (*, '( I3, I4, 3F10.5, 4X, I4, 3F10.5)') I, IE1(I), 
     .    (GEOM1( K, I),K=1,3), IE2( I), (GEOM2( K, I),K=1,3)
   9    CONTINUE
      ENDIF

      NTOTAL = 0
      DO 1 I=1, NUMATM
        IF ( IE1(I).NE.0 .OR. IE2(I).NE.0) THEN
           IGROUP( I) = I
           NTOTAL = NTOTAL + 1
        ELSE
           GOTO 2
        ENDIF
 1    CONTINUE
*
  2   CONTINUE
      I1 = 1
      I2 = 1
 10   CONTINUE
      IF( IE1( I1) .GE. 99 ) THEN
        I1 = I1 + 1
        GOTO 10
      ENDIF

      IF( IE2( I2) .GE. 99 ) THEN
        I2 = I2 + 1
        GOTO 10
      ENDIF

      IF ( IE1( I1) .NE. IE2( I2) ) THEN
         CALL DEBUGR( 'ORIENT: FIRST ATOMS DO NOT MATCH.')
         RETURN
      ENDIF

      DO 100 J= 1, 3
        DISPL( J) = GEOM1( J, I1) - GEOM2( J, I2)
 100  CONTINUE

      DO 105 K= 1, NTOTAL
      DO 105 J= 1, 3
         GEOML( J, K) = GEOM2( J, K) + DISPL( J)
 105  CONTINUE

      J1 = I1+1
      J2 = I2+1

 110  CONTINUE
      IF( IE1( J1) .GE. 99 ) THEN
        J1 = J1 + 1
        GOTO 110
      ENDIF

      IF( IE2( J2) .GE. 99 ) THEN
        J2 = J2 + 1
        GOTO 110
      ENDIF

      IF ( IE1( J1) .NE. IE2( J2) ) THEN
         CALL DEBUGR( 'ORIENT: SECOND ATOMS DO NOT MATCH.')
         RETURN
      ENDIF

      PI = 2.0D0 * ASIN(1.0D0)
      CONV = PI / 180.0D0

* FIRST ROTATE ABOUT Z AXIS TO ACHIEVE VERTICAL COINCIDENCE

      D12 = ( GEOM1( 1, I1) - GEOM1( 1, J1))**2 +
     .      ( GEOM1( 2, I1) - GEOM1( 2, J1))**2 
      D13 = ( GEOM1( 1, I1) - (GEOML( 1, J2)))**2 +
     .      ( GEOM1( 2, I1) - (GEOML( 2, J2)))**2 
      D23 = ( GEOM1( 1, J1) - (GEOML( 1, J2)))**2 +
     .      ( GEOM1( 2, J1) - (GEOML( 2, J2)))**2 
      XY = SQRT( D12*D13)
      IF( XY.LT.1.0D-2) THEN
        CALL DEBUGR( 'ORIENT: ATOMS HAVE ZERO SPACING.' )
        RETURN
      ENDIF

      TEMP = 0.5D0 * (D12+D13-D23) / XY
      IF(TEMP.GT.1.0D0) THEN
        TEMP=1.0D0
      ELSEIF(TEMP.LT.-1.0D0)THEN
        TEMP=-1.0D0
      ENDIF
C??      ANGLES( 1) = -ACOS( TEMP) / CONV
      ANGLES( 1) = ACOS( TEMP) / CONV

* USE CARTESIAN ROTATION ROUTINE WITH NGROUP SET
      CALL RCART( NTOTAL, GEOML, GEOML, I2, 0.0D0, 0.0D0, ANGLES(1),
     .    NTOTAL, IGROUP)


* NOW WORK-OUT THE ROTATION ABOUT Y-AXIS TO MAKE THEM MATCH
      D23 = (GEOM1( 3, I2) - GEOML( 3, J2))**2

      TEMP = 0.5D0 * (D12+D13-D23) / XY
      IF(TEMP.GT.1.0D0) THEN
        TEMP=1.0D0
      ELSEIF(TEMP.LT.-1.0D0)THEN
        TEMP=-1.0D0
      ENDIF
      ANGLES( 2) = -ACOS( TEMP) / CONV

      CALL RCART( NTOTAL, GEOML, GEOML, I2, 0.0D0, ANGLES( 2), 0.0D0,
     .    NTOTAL, IGROUP)

 200  CONTINUE
      K1 = J1 + 1
      K2 = J2 + 1

 210  CONTINUE
      IF( IE1( K1) .GE. 99 ) THEN
        K1 = K1 + 1
        GOTO 210
      ENDIF

      IF( IE2( K2) .GE. 99 ) THEN
        K2 = K2 + 1
        GOTO 210
      ENDIF

      IF ( IE1( K1) .NE. IE2( K2) ) THEN
         CALL DEBUGR( 'ORIENT: SECOND ATOMS DO NOT MATCH.')
         RETURN
      ENDIF
* THE FOLLOWING CODE IS MODIFIED FROM ROUTINE DIHED
*  ATOM I IS GEOML( X, K2)
*       J    GEOM1( X, J1)
*       K    GEOM1( X, I1)
*       L    GEOM1( X, K1)

      XI1 = GEOML( 1, K2) - GEOM1( 1, I1)
      XJ1 = GEOM1( 1, J1) - GEOM1( 1, I1)
      XL1 = GEOM1( 1, K1) - GEOM1( 1, I1)

      YI1 = GEOML( 2, K2) - GEOM1( 2, I1)
      YJ1 = GEOM1( 2, J1) - GEOM1( 2, I1)
      YL1 = GEOM1( 2, K1) - GEOM1( 2, I1)

      ZI1 = GEOML( 3, K2) - GEOM1( 3, I1)
      ZJ1 = GEOM1( 3, J1) - GEOM1( 3, I1)
      ZL1 = GEOM1( 3, K1) - GEOM1( 3, I1)

      DIST = SQRT(XJ1**2+YJ1**2+ZJ1**2)
      COSA = ZJ1 / DIST
      IF(COSA.GT.1.0D0) COSA = 1.0D0
      IF(COSA.LT.-1.0D0) COSA = -1.0D0
      DDD=1.0D0-COSA**2
      IF(DDD.LE.0.0) GOTO 315
      YXDIST=DIST* SQRT(DDD)
      IF(YXDIST.GT.1.0D-9) GOTO 310
  315 CONTINUE
      XI2=XI1
      XL2=XL1
      YI2=YI1
      YL2=YL1
      COSTH=COSA
      SINTH=0.D0
      GOTO 311
  310 COSPH=YJ1/YXDIST
      SINPH=XJ1/YXDIST
      XI2=XI1*COSPH-YI1*SINPH
      XJ2=XJ1*COSPH-YJ1*SINPH
      XL2=XL1*COSPH-YL1*SINPH
      YI2=XI1*SINPH+YI1*COSPH
      YJ2=XJ1*SINPH+YJ1*COSPH
      YL2=XL1*SINPH+YL1*COSPH
*  ROTATE KJ AROUND THE X AXIS SO KJ LIES ALONG THE Z AXIS
      COSTH=COSA
      SINTH=YJ2/DIST
  311 CONTINUE
      YI3=YI2*COSTH-ZI1*SINTH
      YL3=YL2*COSTH-ZL1*SINTH
      CALL DANG(XL2,YL3,XI2,YI3,ANGLES(3))
c?      ANGLES(3)=-ANGLES(3)
c?      IF (ANGLES(3) .LT. -180.0D0) ANGLES(3) = ANGLES(3) + 360.0D0
c?      IF (ANGLES(3) .GT.  180.0D0) ANGLES(3) = ANGLES(3) - 360.0D0

      CALL RPAIR( NTOTAL, GEOML, GEOML, I2, J2, ANGLES(3), 
     .    NTOTAL, IGROUP)
 
* Now that all components are known, do the same thing in reverse and
* allow the rotation routines to accumulate the info in TOTROT.

      DO 400 J= 1, 3
         ANGLES( J) = -ANGLES( J)
 400  CONTINUE

      DO 410 K= 1, NTOTAL
      DO 410 J= 1, 3
        GEOML( J, K) = GEOM1( J, K)
 410  CONTINUE

      CALL RPAIR( NTOTAL, GEOML, GEOML, I1, J1, ANGLES(3), 
     .  0, IGROUP)

      CALL RCART( NTOTAL, GEOML, GEOML, I1, 0.0D0, ANGLES(2), 0.0D0,
     .  0, IGROUP)

      CALL RCART( NTOTAL, GEOML, GEOML, I1, 0.0D0, 0.0D0, ANGLES(1),
     .  0, IGROUP)

      DO 440 K= 1, NTOTAL
      DO 440 J= 1, 3
        GEOML( J, K) = GEOML( J, K) - DISPL( J) 
 440  CONTINUE

c?      WRITE (*,*) ' ORIENT:   geom1                         geom2'
c?      M1 = 0
c?      M2 = 0
c?      DO 499 I=1, 5
c?         M1 = M1 + 1
c?         M2 = M2 + 1
c? 491     IF ( IE1( M1).EQ.99) THEN
c?            M1 = M1 + 1
c?            GOTO 491
c?         ENDIF
c? 492     IF ( IE2( M2).EQ.99) THEN
c?            M2 = M2 + 1
c?            GOTO 492
c?         ENDIF
c?         WRITE (*, '( I3, I4, 3F10.5, 4X, I4, 3F10.5)') 
c?     .    M1,IE1(M1),(GEOML(K,M1),K=1,3),IE2(M2),(GEOM2(K,M2),K=1,3)
c? 499  CONTINUE

 600  CONTINUE
      DO 610 I= 1, NUMATM
        DO 610 J= 1, 3
          GEOM1( J, I) = GEOML( J, I)
 610  CONTINUE

      RETURN
      END
