      SUBROUTINE GMETRY(GEO,COORD)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1NA(NUMATM),NB(NUMATM),NC(NUMATM)
     2       /EULER / TVEC(3,3), ID
      COMMON /REACTN/ STEP, GEOA(3,NUMATM), GEOVEC(3,NUMATM),COLCST
      COMMON /GEOOK/ IGEOOK
      COMMON /NUMCAL/ NUMCAL
      DIMENSION GEO(3,*),COORD(3,*)
      CHARACTER *15 NDIMEN(4)
      LOGICAL GEOOK
      SAVE ICALCN, NDIMEN, GEOOK
      DATA ICALCN/0/
      DATA NDIMEN/' MOLECULE     ',' POLYMER       ',
     1'LAYER STRUCTURE',' SOLID         '/
C***********************************************************************
C
C    GMETRY  COMPUTES COORDINATES FROM BOND-ANGLES AND LENGTHS.
C *** IT IS ADAPTED FROM THE PROGRAM WRITTEN BY M.J.S. DEWAR.
C
C     THREE SEPARATE OPTIONS EXIST WITHIN GMETRY. THESE ARE:
C    (A) IF NA(1) IS EQUAL TO 99 (IMPOSSIBLE UNDER NORMAL CIRCUMSTANCES)
C        THEN GEO IS ASSUMED TO BE IN CARTESIAN RATHER THAN INTERNAL
C        COORDINATES, AND COORD IS THEN SET EQUAL TO GEO.
C    (B) IF STEP IS NON-ZERO (THIS IS THE CASE WHEN "SADDLE" IS USED)
C        THEN GEO IS FIRST MODIFIED BY SHIFTING THE INTERNAL COORDINATES
C        ALONG A RADIUS FROM GEOA TO PLACE GEO AT ADISTANCESTEPFROMGEOA.
C    (C) NORMAL CONVERSION FROM INTERNAL TO CARTESIAN COORDINATESISDONE.
C
C  ON INPUT:
C         GEO    = ARRAY OF INTERNAL COORDINATES.
C         NATOMS = NUMBER OF ATOMS, INCLUDING DUMMIES.
C         NA     = ARRAY OF ATOM LABELS FOR BOND LENGTHS.
C
C  ON OUTPUT:
C         COORD  = ARRAY OF CARTESIAN COORDINATES
C
C***********************************************************************
C                                     OPTION (B)
      GEOOK=(IGEOOK.EQ.99)
      IF(ABS(STEP) .GT. 1.D-4) THEN
         SUM=0.D0
         DO 10 J=1,3
C$DOIT VBEST
            DO 10 I=1,NATOMS
               GEOVEC(J,I)=GEO(J,I)-GEOA(J,I)
   10    SUM=SUM+GEOVEC(J,I)**2
         SUM=SQRT(SUM)
         ERROR=(SUM-STEP)/SUM
      ELSE
         ERROR=0.D0
      ENDIF
      DO 20 J=1,3
C$DOIT VBEST
         DO 20 I=1,NATOMS
   20 GEO(J,I)=GEO(J,I)-ERROR*GEOVEC(J,I)
C                                     OPTION (A)
      IF(NA(1).EQ.99) THEN
         DO 30 I=1,3
C$DOIT VBEST
            DO 30 J=1,NATOMS
   30    COORD(I,J)=GEO(I,J)
         GOTO 100
      ENDIF
C                                     OPTION (C)
      COORD(1,1)=0.0D00
      COORD(2,1)=0.0D00
      COORD(3,1)=0.0D00
      COORD(1,2)=GEO(1,2)
      COORD(2,2)=0.0D00
      COORD(3,2)=0.0D00
      IF(NATOMS.EQ.2) GOTO 100
      CCOS=COS(GEO(2,3))
      IF(NA(3).EQ.1)THEN
         COORD(1,3)=COORD(1,1)+GEO(1,3)*CCOS
      ELSE
         COORD(1,3)=COORD(1,2)-GEO(1,3)*CCOS
      ENDIF
      COORD(2,3)=GEO(1,3)*SIN(GEO(2,3))
      COORD(3,3)=0.0D00
      DO 90 I=4,NATOMS
         COSA=COS(GEO(2,I))
         MB=NB(I)
         MC=NA(I)
         XB=COORD(1,MB)-COORD(1,MC)
         YB=COORD(2,MB)-COORD(2,MC)
         ZB=COORD(3,MB)-COORD(3,MC)
         RBC=XB*XB+YB*YB+ZB*ZB
         IF(RBC.LT.1.D-16)THEN
C
C     TWO ATOMS ARE COINCIDENT.  A FATAL ERROR.
C
            WRITE(6,'(A,I4,A,I4,A)')' ATOMS',MB,' AND',MC,' ARE COINCIDE
     1NT'
            WRITE(6,'(A)')' THIS IS A FATAL ERROR, RUN STOPPED IN GMETRY
     1'
            STOP
         ELSE
            RBC=1.0D00/SQRT(RBC)
         ENDIF
         MA=NC(I)
         XA=COORD(1,MA)-COORD(1,MC)
         YA=COORD(2,MA)-COORD(2,MC)
         ZA=COORD(3,MA)-COORD(3,MC)
C
C     ROTATE ABOUT THE Z-AXIS TO MAKE YB=0, AND XB POSITIVE.  IF XYB IS
C     TOO SMALL, FIRST ROTATE THE Y-AXIS BY 90 DEGREES.
C
         XYB=SQRT(XB*XB+YB*YB)
         K=-1
         IF (XYB.GT.0.1D00) GO TO 40
         XPA=ZA
         ZA=-XA
         XA=XPA
         XPB=ZB
         ZB=-XB
         XB=XPB
         XYB=SQRT(XB*XB+YB*YB)
         K=+1
C
C     ROTATE ABOUT THE Y-AXIS TO MAKE ZB VANISH
C
   40    COSTH=XB/XYB
         SINTH=YB/XYB
         XPA=XA*COSTH+YA*SINTH
         YPA=YA*COSTH-XA*SINTH
         SINPH=ZB*RBC
         COSPH=SQRT(ABS(1.D00-SINPH*SINPH))
         ZQA=ZA*COSPH-XPA*SINPH
C
C     ROTATE ABOUT THE X-AXIS TO MAKE ZA=0, AND YA POSITIVE.
C
         YZA=SQRT(YPA**2+ZQA**2)
         IF(YZA.LT.1.D-4)GOTO 60
         IF(YZA.LT.2.D-2 .AND. .NOT.GEOOK)THEN
            WRITE(6,'(//20X,'' CALCULATION ABANDONED AT THIS POINT'')')
            WRITE(6,'(//10X,'' THREE ATOMS BEING USED TO DEFINE THE'',/
     110X,'' COORDINATES OF A FOURTH ATOM, WHOSE BOND-ANGLE IS'')')
            WRITE(6,'(10X,'' NOT ZERO OR 180 DEGREEES, ARE '',
     1''IN AN ALMOST STRAIGHT'')')
            WRITE(6,'(10X,'' LINE.  THERE IS A HIGH PROBABILITY THAT THE
     1'',/10X,'' COORDINATES OF THE ATOM WILL BE INCORRECT.'')')
            WRITE(6,'(//20X,''THE FAULTY ATOM IS ATOM NUMBER'',I4)')I
            CALL GEOUT(1)
            WRITE(6,'(//20X,''CARTESIAN COORDINATES UP TO FAULTY ATOM'')
     1')
            WRITE(6,'(//5X,''I'',12X,''X'',12X,''Y'',12X,''Z'')')
            DO 50 J=1,I
   50       WRITE(6,'(I6,F16.5,2F13.5)')J,(COORD(K,J),K=1,3)
            WRITE(6,'(//6X,'' ATOMS'',I3,'','',I3,'', AND'',I3,
     1'' ARE WITHIN'',F7.4,'' ANGSTROMS OF A STRAIGHT LINE'')')
     2MC,MB,MA,YZA
            STOP
         ENDIF
         COSKH=YPA/YZA
         SINKH=ZQA/YZA
         GOTO 70
   60    CONTINUE
C
C   ANGLE TOO SMALL TO BE IMPORTANT
C
         COSKH=1.D0
         SINKH=0.D0
   70    CONTINUE
C
C     COORDINATES :-   A=(???,YZA,0),   B=(RBC,0,0),  C=(0,0,0)
C     NONE ARE NEGATIVE.
C     THE COORDINATES OF I ARE EVALUATED IN THE NEW FRAME.
C
         SINA=SIN(GEO(2,I))
         SIND=-SIN(GEO(3,I))
         COSD=COS(GEO(3,I))
         XD=GEO(1,I)*COSA
         YD=GEO(1,I)*SINA*COSD
         ZD=GEO(1,I)*SINA*SIND
C
C     TRANSFORM THE COORDINATES BACK TO THE ORIGINAL SYSTEM.
C
         YPD=YD*COSKH-ZD*SINKH
         ZPD=ZD*COSKH+YD*SINKH
         XPD=XD*COSPH-ZPD*SINPH
         ZQD=ZPD*COSPH+XD*SINPH
         XQD=XPD*COSTH-YPD*SINTH
         YQD=YPD*COSTH+XPD*SINTH
         IF (K.LT.1) GO TO 80
         XRD=-ZQD
         ZQD=XQD
         XQD=XRD
   80    COORD(1,I)=XQD+COORD(1,MC)
         COORD(2,I)=YQD+COORD(2,MC)
         COORD(3,I)=ZQD+COORD(3,MC)
   90 CONTINUE
C
C *** NOW REMOVE THE TRANSLATION VECTORS, IF ANY, FROM THE ARRAY COOR
C
  100 CONTINUE
      K=NATOMS
  110 IF(LABELS(K).NE.107) GOTO 120
      K=K-1
      GOTO 110
  120 K=K+1
      IF(K.GT.NATOMS) GOTO 170
C
C   SYSTEM IS A SOLID, OF DIMENSION NATOMS+1-K
C
      L=0
      DO 130 I=K,NATOMS
         L=L+1
         MC=NA(I)
         TVEC(1,L)=COORD(1,I)-COORD(1,MC)
         TVEC(2,L)=COORD(2,I)-COORD(2,MC)
         TVEC(3,L)=COORD(3,I)-COORD(3,MC)
  130 CONTINUE
      ID=L
      IF (ICALCN.NE.NUMCAL) THEN
         ICALCN=NUMCAL
         WRITE(6,140)NDIMEN(ID+1)
  140    FORMAT(/10X,'    THE SYSTEM IS A ',A15,/)
         IF(ID.EQ.0) GOTO 170
         WRITE(6,150)
         WRITE(6,160)(I,(TVEC(J,I),J=1,3),I=1,ID)
  150    FORMAT(/,'                UNIT CELL TRANSLATION VECTORS',/
     1/,'              X              Y              Z')
  160    FORMAT('    T',I1,' = ',F11.7,'    ',F11.7,'    ',F11.7)
      ENDIF
  170 CONTINUE
      J=0
      DO 190 I=1,NATOMS
         IF (LABELS(I).EQ.99.OR.LABELS(I).EQ.107) GO TO 190
         J=J+1
C$DOIT ASIS
         DO 180 K=1,3
  180    COORD(K,J)=COORD(K,I)
  190 CONTINUE
      RETURN
      END
