      SUBROUTINE GETGEO(IREAD,LABELS,GEO,LOPT,NA,NB,NC,AMS,NATOMS,INT)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension GEO(3,*),NA(*),NB(*),NC(*),AMS(*), LOPT(3,*)
     1,LABELS(*)
      logical INT
************************************************************************
*
*   GETGEO READS IN THE GEOMETRY. THE ELEMENT IS SPECIFIED BY IT'S
*          CHEMICAL SYMBOL, OR, OPTIONALLY, BY IT'S ATOMIC NUMBER.
*
*  ON INPUT   IREAD  = CHANNEL NUMBER FOR READ, NORMALLY 5
*             AMS    = DEFAULT ATOMIC MASSES.
*
* ON OUTPUT LABELS = ATOMIC NUMBERS OF ALL ATOMS, INCLUDING DUMMIES.
*           GEO    = INTERNAL COORDINATES, IN ANGSTROMS, AND DEGREES.
*           LOPT   = INTEGER ARRAY, A '1' MEANS OPTIMIZE THIS PARAMETER,
*                    '0' MEANS do NOT OPTIMIZE, AND A '-1' LABELS THE
*                    REACTION COORDINATE.
*           NA     = INTEGER ARRAY OF ATOMS (SEE DATA INPUT)
*           NB     = INTEGER ARRAY OF ATOMS (SEE DATA INPUT)
*           NC     = INTEGER ARRAY OF ATOMS (SEE DATA INPUT)
*           ATMASS = ATOMIC MASSES OF ATOMS.
************************************************************************
      common /PATH  / IDUM(2),REACT(3,66), DUMM1,DUMM2
      common /SIMBOL/ SIMBOL(MAXPAR)
      common /ATMASS/ ATMASS(NUMATM)
      common /ATOMTX/ LTXT, TXTATM(NUMATM)
      common /KEYWRD/ KEYWRD
      dimension ISTART(40), XYZ(3,NUMATM), VALUE(4)
      logical LEADSP, IRCDRC
      character KEYWRD*241, TXTATM*8, SIMBOL*10, LTXT*1
      character ELEMNT(107)*2, LINE*80, SPACE*1, NINE*1,ZERO*1,
     1TAB*1, COMMA*1, STRING*80, ELE*2, TURN*1
      save ELEMNT, COMMA, SPACE, NINE, ZERO
      data (ELEMNT(I),I=1,107)/'H','HE',
     1 'LI','BE','B','C','N','O','F','NE',
     2 'NA','MG','AL','SI','P','S','CL','AR',
     3 'K','CA','SC','TI','V','CR','MN','FE','CO','NI','CU',
     4 'ZN','GA','GE','AS','SE','BR','KR',
     5 'RB','SR','Y','ZR','NB','MO','TC','RU','RH','PD','AG',
     6 'CD','IN','SN','SB','TE','I','XE',
     7 'CS','BA','LA','CE','PR','ND','PM','SM','EU','GD','TB','DY',
     8 'HO','ER','TM','YB','LU','HF','TA','W','RE','OS','IR','PT',
     9 'AU','HG','TL','PB','BI','PO','AT','RN',
     1 'FR','RA','AC','TH','PA','U','NP','PU','AM','CM','BK','CF','XX',
     2 'FM','MD','CB','++','+','--','-','TV'/
      data COMMA,SPACE,NINE,ZERO/',',' ','9','0'/
      TAB=CHAR(9)
      IRCDRC=(INDEX(KEYWRD,'IRC')+INDEX(KEYWRD,'DRC') .NE.0)
      ILOWA = ICHAR('a')
      ILOWZ = ICHAR('z')
      ICAPA = ICHAR('A')
      ICAPZ = ICHAR('Z')
      MAXTXT=0
      NATOMS=0
      NUMAT=0
      ISERR=0
      do 10 I=1,MAXPAR
   10 SIMBOL(I)= '---'
   20 read(IREAD,'(A)',end=130,ERR=230)LINE
      if (LINE.EQ.' ') go to 130
      NATOMS=NATOMS+1
C
C   SEE IF TEXT IS ASSOCIATED WITH THIS ELEMENT
C
      I=INDEX(LINE,'(')
      if (I.NE.0) then
C
C  YES, ELEMENT IS LABELLED.
C
         K=INDEX(LINE,')')
         TXTATM(NATOMS)=LINE(I:K)
         MAXTXT=MAX(MAXTXT,K-I+1)
         STRING=LINE(1:I-1)//LINE(K+1:)
         LINE=STRING
      else
         TXTATM(NATOMS)=' '
      end if
*   CLEAN THE INPUT DATA
************************************************************************
      do 30 I=1,80
         ILINE=ICHAR(LINE(I:I))
         if (ILINE.GE.ILOWA.AND.ILINE.LE.ILOWZ) then
            LINE(I:I)=CHAR(ILINE+ICAPA-ILOWA)
         end if
   30 continue
************************************************************************
      ICOMMA=ICHAR(COMMA)
      ITAB=ICHAR(TAB)
      do 40 I=1,80
         KHAR=ICHAR(LINE(I:I))
         if (KHAR.EQ.ICOMMA.OR.KHAR.EQ.ITAB) LINE(I:I)=SPACE
   40 continue
*
*   INITIALIZE ISTART TO INTERPRET BLANKS AS ZERO'S
      do 50 I=1,10
   50 ISTART(I)=80
*
* FIND INITIAL DIGIT OF ALL NUMBERS, CHECK FOR LEADING SPACES FOLLOWED
*     BY A character AND STORE IN ISTART
      LEADSP=.TRUE.
      NVALUE=0
      do 60 I=1,80
         if (LEADSP.AND.LINE(I:I).NE.SPACE) then
            NVALUE=NVALUE+1
            ISTART(NVALUE)=I
         end if
         LEADSP=(LINE(I:I).EQ.SPACE)
   60 continue
*
* ESTABLISH THE ELEMENT'S NAME AND ISOTOPE, CHECK FOR ERRORS OR E.O.DATA
*
      WEIGHT=0.D0
      STRING=LINE(ISTART(1):ISTART(2)-1)
      if ( STRING(1:1) .GE. ZERO .AND. STRING(1:1) .LE. NINE) then
*  ATOMIC NUMBER USED: NO ISOTOPE ALLOWED
         LABEL = READA(STRING,1)
         if (LABEL.EQ.0) go to 120
         if (LABEL.LT.0.OR.LABEL.GT.107) then
            write(MIOB,'(''  ILLEGAL ATOMIC NUMBER'')')
            go to 240
         end if
         go to 80
      end if
*  ATOMIC SYMBOL USED
      real=ABS(READA(STRING,1))
      if (REAL.LT.1.D-15) then
*   NO ISOTOPE
         ELE=STRING(1:2)
      else
         WEIGHT=REAL
         if ( STRING(2:2) .GE. ZERO .AND. STRING(2:2) .LE. NINE) then
            ELE=STRING(1:1)
         else
            ELE=STRING(1:2)
         end if
      end if
*   CHECK FOR ERROR IN ATOMIC SYMBOL
      if (ELE(1:1).EQ.'-'.AND.ELE(2:2).NE.'-')  ELE(2:2)=' '
      do 70 I=1,107
         if (ELE.EQ.ELEMNT(I)) then
            LABEL=I
            go to 80
         end if
   70 continue
      if (ELE(1:1).EQ.'X') then
         LABEL=99
         go to 80
      end if
      write(MIOB,'(''  UNRECOGNIZED ELEMENT NAME: ('',A,'')'')') ELE
      go to 240
*
* ALL O.K.
*
   80 if (LABEL.NE.99) NUMAT=NUMAT+1
      if (WEIGHT.NE.0.D0) then
         write(MIOB,'('' FOR ATOM'',I4,''  ISOTOPIC MASS:''
     %    ,F15.5)') NATOMS, WEIGHT
         ATMASS(NUMAT)=WEIGHT
      else
         if (LABEL .NE. 99)  ATMASS(NUMAT)=AMS(LABEL)
      end if
      if (NATOMS.GT.NUMATM) then
         write(MIOB,
     %         '(//10X,''****  MAX. NUMBER OF ATOMS ALLOWED:'',I4)')
     %         NUMATM
         call filclo ! J/M
         stop ' in GETGEO - too many atoms'
      end if
      LABELS(NATOMS)   =LABEL
      GEO(1,NATOMS)    = READA(LINE,ISTART(2))
      GEO(2,NATOMS)    = READA(LINE,ISTART(4))
      GEO(3,NATOMS)    = READA(LINE,ISTART(6))
      if (IRCDRC) then
         TURN=LINE(ISTART(3):ISTART(3))
         if (TURN.EQ.'T') then
            LOPT(1,NATOMS)=1
            if (NATOMS.EQ.1)
     %        write(MIOB,'(A)') ' IN DRC MONITOR POTENTIAL ENERGY'//
     %                          ' TURNING POINTS'
         else
            LOPT(1,NATOMS)=0
         end if
         TURN=LINE(ISTART(5):ISTART(5))
         if (TURN.EQ.'T') then
            LOPT(2,NATOMS)=1
         else
            LOPT(2,NATOMS)=0
         end if
         TURN=LINE(ISTART(7):ISTART(7))
         if (TURN.EQ.'T') then
            LOPT(3,NATOMS)=1
         else
            LOPT(3,NATOMS)=0
         end if
      else
         LOPT(1,NATOMS)   = READA(LINE,ISTART(3))
         LOPT(2,NATOMS)   = READA(LINE,ISTART(5))
         LOPT(3,NATOMS)   = READA(LINE,ISTART(7))
         do 90 I=3,7,2
            if (ICHAR(LINE(ISTART(I):ISTART(I))).GE.ICAPA .AND.
     1         ICHAR(LINE(ISTART(I):ISTART(I))).LE.ICAPZ)  ISERR=1
   90    continue
      end if
      NA(NATOMS)       = READA(LINE,ISTART(8))
      NB(NATOMS)       = READA(LINE,ISTART(9))
      NC(NATOMS)       = READA(LINE,ISTART(10))
C
C  SPECIAL CASE OF USERS FORGETTING TO ADD DIHEDRAL DATA FOR ATOM 3
C
      if (NATOMS.EQ.3) then
         if (LOPT(3,3).EQ.2) then
            NA(3)=1
            NB(3)=2
            GEO(3,3)=0.D0
            LOPT(3,3)=0
         else if (LOPT(3,3).EQ.1.AND.ABS(GEO(3,3)-2.D0).LT.1.D-4) then
            NA(3)=2
            NB(3)=1
            GEO(3,3)=0.D0
            LOPT(3,3)=0
         end if
      end if
      if (LOPT(1,NATOMS).GT.1.OR.LOPT(2,NATOMS).GT.1 .OR.
     1   LOPT(3,NATOMS).GT.1) ISERR=1
      if (ISERR.EQ.1) then
C
C  MUST BE GAUSSIAN GEOMETRY INPUT
C
         do 110 I=2,NATOMS
            do 110 K=1,3
               J=GEO(K,I)+0.4D0
               if (ABS(GEO(K,I)-J).GT.1.D-5) then
C
C   GEOMETRY CANNOT BE GAUSSIAN
C
                  write(MIOB,'(A)')
     %             ' GEOMETRY IS FAULTY.  GEOMETRY READ IN IS'
                  CONST=3.141592653598D0/180.D0
                  do 100 L=1,NATOMS
                     GEO(2,L)=GEO(2,L)*CONST
  100             GEO(3,L)=GEO(3,L)*CONST
                  call GEOUT(MIOB)
c
                  call filclo ! J/M
                  stop ' in GETGEO - geometry is faulty'
c
               end if
  110    continue
         NATOMS=-1
         return
      end if
      go to 20
*
* ALL DATA READ IN, CLEAN UP AND return
*
  120 NATOMS=NATOMS-1
  130 NA(2)=1
      LTXT=CHAR(MAXTXT)
      if (NATOMS.GT.3) then
         INT=(NA(4).NE.0)
      else
         if (GEO(2,3).LT.10.AND.NATOMS.EQ.3)
     1  write(MIOB,'(//10X,
     %     '' WARNING: INTERNAL COORDINATES ARE ASSUMED -'',/
     % 10X,'' FOR THREE-ATOM SYSTEMS '',//)')
         INT=.TRUE.
      end if
      if (INT)  GEO(2,2)=0
C
C     READ IN VELOCITY VECTOR, IF PRESENT
C
      if (INDEX(KEYWRD,'VELO').GT.0) then
         if (INT) then
            write(MIOB,
     %       '(A)')' COORDINATES MUST BE CARTESIAN WHEN VELOCITY'
     %     //' VECTOR IS USED.'
c
            call filclo ! J/M
            stop ' in GETGEO - cartesian coords must be used'
c
         end if
C#      write(MIOB,'(/10X,A)') 'INITIAL VELOCITY VECTOR FOR DRC'
         do 150 I=1,NATOMS
            read(MIOA,'(A)') LINE
            call NUCHAR(LINE,VALUE,NDMY)
            if (NDMY.NE.3) then
               write(MIOB,'(/10X,A)')
     %           '  THERE MUST BE EXACTLY THREE VELOCITY DATA PER LINE'
c
               call filclo ! J/M
               stop ' in GETGEO - three velocity data/line reqrd.'
c
            end if
            do 140 J=1,3
  140       REACT(J,I+2)=VALUE(J)
C#      write(MIOB,'(2X,A2,2X,3F13.5)')
C#     %      ELEMNT(LABELS(I)),(VALUE(J),J=1,3)
  150    continue
         do 160 I=1,3
            do 160 J=1,2
  160    REACT(I,J)=GEO(I,J+1)-GEO(I,1)
C
C  NOW TO ROTATE VELOCITY VECTOR TO SUIT INTERNAL COORDINATE DEFINITION
C
C
C   ROTATE AROUND THE 1-2 X-AXIS TO AS TO ELIMINATE REACT(3,2)
C   (PUT ATOM 2 IN X-Y PLANE)
         SA=REACT(3,1)/SQRT(REACT(2,1)**2+REACT(3,1)**2+1.D-20)
         CA=SIGN(SQRT(1.D0-SA**2),REACT(2,1))
C#      LABELS(NATOMS+1)=1
C#      LABELS(NATOMS+2)=1
C#      write(MIOB,*) ' FIRST ROTATION, ABOUT 1-2 X-AXIS'
         do 170 I=1,NATOMS+2
            TEMP1= REACT(2,I)*CA+REACT(3,I)*SA
            TEMP2=-REACT(2,I)*SA+REACT(3,I)*CA
            REACT(2,I)=TEMP1
            REACT(3,I)=TEMP2
C#      write(MIOB,'(2X,A2,2X,3F13.5)')
C#     %      ELEMNT(LABELS(I)),(REACT(J,I),J=1,3)
  170    continue
C   ROTATE AROUND THE 1-2 Z-AXIS TO AS TO ELIMINATE REACT(2,2)
C   (PUT ATOM 2 ON X AXIS)
         CA=REACT(1,1)/SQRT(REACT(2,1)**2+REACT(1,1)**2+1.D-20)
         SA=SIGN(SQRT(1.D0-CA**2),REACT(2,1))
C#      write(MIOB,*) ' SECOND ROTATION, ABOUT 1-2 Z-AXIS'
         do 180 I=1,NATOMS+2
            TEMP1= REACT(1,I)*CA+REACT(2,I)*SA
            TEMP2=-REACT(1,I)*SA+REACT(2,I)*CA
            REACT(1,I)=TEMP1
            REACT(2,I)=TEMP2
C#      write(MIOB,'(2X,A2,2X,3F13.5)')
C#     %      ELEMNT(LABELS(I)),(REACT(J,I),J=1,3)
  180    continue
C   ROTATE AROUND THE 2-3 X-AXIS TO AS TO ELIMINATE REACT(3,3)
C   (PUT ATOM 3 ON X-Y PLANE)
         SA=REACT(3,2)/SQRT(REACT(2,2)**2+REACT(3,2)**2+1.D-20)
         CA=SIGN(SQRT(1.D0-SA**2),REACT(2,2))
C#      write(MIOB,*) ' THIRD ROTATION, ABOUT 2-3 X-AXIS'
         do 190 I=1,NATOMS+2
            TEMP1= REACT(2,I)*CA+REACT(3,I)*SA
            TEMP2=-REACT(2,I)*SA+REACT(3,I)*CA
            REACT(2,I)=TEMP1
            REACT(3,I)=TEMP2
C#      write(MIOB,'(2X,A2,2X,3F13.5)')
C#     %      ELEMNT(LABELS(I)),(REACT(J,I),J=1,3)
  190    continue
C
C  STRIP OFF FIRST TWO COORDINATES; THESE WERE THE COORDINATE AXIS
C  DEFINITIONS
C
         do 200 I=1,NATOMS
            do 200 J=1,3
  200    REACT(J,I)=REACT(J,I+2)
      end if
      if (  .NOT. INT ) then
         do 210 I=1,NATOMS
            do 210 J=1,3
  210    XYZ(J,I)=GEO(J,I)
         DEGREE=90.D0/ASIN(1.D0)
         call XYZINT(XYZ,NATOMS,NA,NB,NC,DEGREE,GEO)
         if (INDEX(KEYWRD,' XYZ').EQ.0) then
C
C  UNCONDITIONALLY SET FLAGS FOR INTERNAL COORDINATES
C
            do 220 I=1,3
               do 220 J=I,3
  220       LOPT(J,I)=0
         end if
         if (ABS(GEO(2,3)-180.D0).LT.1.D-4 .OR.
     1      ABS(GEO(2,3)).LT.1.D-4)           then
            write(MIOB,'(A)')
     %       ' DUE TO PROGRAM BUG, THE FIRST THREE ATOMS MUST'//
     %       ' NOT LIE IN A STRAIGHT LINE.'
c
            call filclo ! J/M
            stop ' in GETGEO - first three atoms form a straight line'
c
         end if
      else if (.NOT.IRCDRC) then
         LOPT(2,2)=0
         if (LOPT(1,1)+LOPT(2,1)+LOPT(3,1)+LOPT(3,2)+
     1      LOPT(3,3) .GT. 0) then
c
            LOPT(1,1)=0
            LOPT(2,1)=0
            LOPT(3,1)=0
            LOPT(3,2)=0
            LOPT(3,3)=0
            write(MIOB,
     %       '(//10X,'' AN UNOPTIMIZABLE GEOMETRIC PARAMETER HAS''/
     %           10X,'' BEEN MARKED FOR OPTIMIZATION. THIS IS A'',
     %               '' NON-FATAL ERROR'')')
         end if
      end if
      if (NA(3).EQ.0) then
         NB(3)=1
         NA(3)=2
      end if
      return
* ERROR CONDITIONS
  230 if (IREAD.EQ.5) then
         write(MIOB,'( '' ERROR DURING READ AT ATOM NUMBER '', I3 )')
     %           NATOMS
      else
         NATOMS=0
         return
      end if
  240 J=NATOMS-1
      write(MIOB,'('' DATA CURRENTLY READ IN ARE '')')
      do 250 K=1,J
  250 write(MIOB,260) LABELS(K),(GEO(JJ,K),LOPT(JJ,K),JJ=1,3),
     1NA(K),NB(K),NC(K)
  260 format(I4,2X,3(F10.5,2X,I2,2X),3(I2,1X))
c
      call filclo ! J/M
      stop ' in GETGEO - '
c
      end
