      SUBROUTINE READMO
      implicit double precision (A-H, O-Z)
      include 'SIZES'
      include 'MOP6FILE'
C
C MODULE TO READ IN GEOMETRY FILE, OUTPUT IT TO THE USER,
C AND CHECK THE DATA TO SEE IF IT IS REASONABLE.
C EXIT IF NECESSARY.
C
C
C
C  ON EXIT NATOMS    = NUMBER OF ATOMS PLUS DUMMY ATOMS (IF ANY).
C          KEYWRD    = KEYWORDS TO CONTROL CALCULATION
C          KOMENT    = COMMENT CARD
C          TITLE     = TITLE CARD
C          LABELS    = ARRAY OF ATOMIC LABELS INCLUDING DUMMY ATOMS.
C          GEO       = ARRAY OF INTERNAL COORDINATES.
C          LOPT      = FLAGS FOR OPTIMIZATION OF MOLECULE
C          NA        = ARRAY OF LABELS OF ATOMS, BOND LENGTHS.
C          NB        = ARRAY OF LABELS OF ATOMS, BOND ANGLES.
C          NC        = ARRAY OF LABELS OF ATOMS, DIHEDRAL ANGLES.
C          LATOM     = LABEL OF ATOM OF REACTION COORDINATE.
C          LPARAM    = RC: 1 FOR LENGTH, 2 FOR ANGLE, AND 3 FOR DIHEDRAL
C          REACT(200)= REACTION COORDINATE PARAMETERS
C          LOC(1,I)  = LABEL OF ATOM TO BE OPTIMIZED.
C          LOC(2,I)  = 1 FOR LENGTH, 2 FOR ANGLE, AND 3 FOR DIHEDRAL.
C          NVAR      = NUMBER OF PARAMETERS TO BE OPTIMIZED.
C          XPARAM    = STARTING VALUE OF PARAMETERS TO BE OPTIMIZED.
C
************************************************************************
C *** INPUT THE TRIAL GEOMETRY  \IE.  KGEOM=0\
C   LABEL(I) = THE ATOMIC NUMBER OF ATOM\I\.
C            = 99, THEN THE I-TH ATOM IS A DUMMY ATOM USED ONLY TO
C              SIMPLIFY THE DEFINITION OF THE MOLECULAR GEOMETRY.
C   GEO(1,I) = THE INTERNUCLEAR SEPARATION \IN ANGSTROMS\ BETWEEN ATOMS
C              NA(I) AND (I).
C   GEO(2,I) = THE ANGLE NB(I):NA(I):(I) INPUT IN DEGREES; STORED IN
C              RADIANS.
C   GEO(3,I) = THE ANGLE BETWEEN THE VECTORS NC(I):NB(I) AND NA(I):(I)
C              INPUT IN DEGREES - STORED IN RADIANS.
C  LOPT(J,I) = -1 IF GEO(J,I) IS THE REACTION COORDINATE.
C            = +1 IF GEO(J,I) IS A PARAMETER TO BE OPTIMIZED
C            =  0 OTHERWISE.
C *** NOTE:    MUCH OF THIS DATA IS NOT includeD FOR THE FIRST 3 ATOMS.
C     ATOM1  INPUT LABELS(1) ONLY.
C     ATOM2  INPUT LABELS(2) AND GEO(1,2) SEPARATION BETWEEN ATOMS 1+2
C     ATOM3  INPUT LABELS(3), GEO(1,3)    SEPARATION BETWEEN ATOMS 2+3
C              AND GEO(2,3)              ANGLE ATOM1 : ATOM2 : ATOM3
C
************************************************************************
C
      dimension LOPT(3,NUMATM)
      character KEYWRD*241, KOMENT*81, TITLE*81, LINE*80, BANNER*80
      character KEYS(80)*1, SPACE*1, SPACE2*2, CH*1, CH2*2
      character ELEMNT*2, IDATE*28, GETNAM*80
      common /KEYWRD/ KEYWRD
      common /TITLES/ KOMENT,TITLE
      common /GEOVAR/ NVAR, LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
      common /PATH  / LATOM,LPARAM,REACT(200)
      common /MESH  / LATOM1, LPARA1, LATOM2, LPARA2
      common /ELEMTS/ ELEMNT(107)
      common /OKMANY/ ISOK
      common /ISTOPE/ AMS(107)
      common /GEOM  / GEO(3,NUMATM)
      common /NUMCAL/ NUMCAL
      common /GEOKST/ NATOMS,LABELS(NUMATM),
     1NA(NUMATM),NB(NUMATM),NC(NUMATM)
      common /GEOSYM/ NDEP, LOCPAR(MAXPAR), IDEPFN(MAXPAR),
     1                      LOCDEP(MAXPAR)
      logical INT, AIGEO, ISOK
      save SPACE, SPACE2, IREACT, INT
      dimension COORD(3,NUMATM),VALUE(40)
      equivalence (KEYS(1),KEYWRD)
      data SPACE, SPACE2/' ','  '/
      CONVTR=2.D0*ASIN(1.D0)/180.D0
      AIGEO=.FALSE.
   10 continue
C
      call GETTXT
      if (INDEX(KEYWRD,'ECHO').NE.0) then
         rewind MIOA
         ISOK=.FALSE.
         do 50 I=1,1000
            read(MIOA,'(A)',end=60)KEYWRD
            do 20 J=80,2,-1
   20       if (KEYWRD(J:J).NE.' ') go to 30
            J=1
   30       do 40 K=1,J
   40       if (ICHAR(KEYWRD(K:K)).LT.32) KEYWRD(K:K)='*'
            write(MIOB,'(1X,A)')  KEYWRD(1:J)
   50    continue
   60    continue
         rewind MIOA
         call GETTXT
      end if
      if (INDEX(KEYWRD,'ECHO').NE.0) write(MIOB,'(''1'')')
      if (KEYWRD(1:1) .NE. SPACE) then
         CH=KEYWRD(1:1)
         KEYWRD(1:1)=SPACE
         do 70 I=2,239
            CH2=KEYWRD(I:I)
            KEYWRD(I:I)=CH
            CH=CH2
            if (KEYWRD(I+1:I+2) .EQ. SPACE2) then
               KEYWRD(I+1:I+1)=CH
               go to 80
            end if
   70    continue
         CH2=KEYWRD(240:240)
         KEYWRD(240:240)=CH
         KEYWRD(241:241)=CH2
   80    continue
      end if
      if (KOMENT(1:1) .NE. SPACE) then
         CH=KOMENT(1:1)
         KOMENT(1:1)=SPACE
         do 90 I=2,79
            CH2=KOMENT(I:I)
            KOMENT(I:I)=CH
            CH=CH2
            if (KOMENT(I+1:I+2) .EQ. SPACE2) then
               KOMENT(I+1:I+1)=CH
               go to 100
            end if
   90    continue
         CH2=KOMENT(80:80)
         KOMENT(80:80)=CH
         KOMENT(81:81)=CH2
  100    continue
      end if
      if (TITLE(1:1) .NE. SPACE) then
         CH=TITLE(1:1)
         TITLE(1:1)=SPACE
         do 110 I=2,79
            CH2=TITLE(I:I)
            TITLE(I:I)=CH
            CH=CH2
            if (TITLE(I+1:I+2) .EQ. SPACE2) then
               TITLE(I+1:I+1)=CH
               go to 120
            end if
  110    continue
         CH2=TITLE(80:80)
         TITLE(80:80)=CH
         TITLE(81:81)=CH2
  120    continue
      end if
      do 121 I=1,200
  121 REACT(I)=0.D0
      LATOM=0
      LPARAM=0
      if (INDEX(KEYWRD,'OLDGEO').EQ.0) then
         NVAR=0
         NDEP=0
         if (AIGEO.OR.INDEX(KEYWRD,'AIGIN').NE.0) then
            call GETGEG(MIOA,LABELS,GEO,LOPT,NA,NB,NC,AMS,NATOMS,INT)
         else
            call GETGEO(MIOA,LABELS,GEO,LOPT,NA,NB,NC,AMS,NATOMS,INT)
            if (NATOMS.LT.0) then
               rewind MIOA
               if (NUMCAL.NE.1) then
                  write(MIOB,'(//,A)')
     %             '   GAUSSIAN INPUT REQUIRES STAND-ALONE JOB'
                  write(MIOB,'(/,A)') '   OR KEYWORD "AIGIN"'
c
                  call filclo  !  J/M
                  stop ' in READMO - gaussian input'
               end if
               AIGEO=.TRUE.
               go to 10
            end if
         end if
c
         if (NATOMS.EQ.0) then
           call filclo  !  J/M
           stop ' in READMO - there are NO atoms!'
         end if
c
      else
         DEGREE=90.D0/ASIN(1.D0)
         do 130 I=1,NATOMS
            do 130 J=2,3
  130    GEO(J,I)=GEO(J,I)*DEGREE
      end if
      if (INDEX(KEYWRD,'FORCE').NE.0 .AND. LABELS(NATOMS).EQ.107) then
      do 131 I=1,NA(NATOMS)
      if (LABELS(I).EQ.99) then
      write(MIOB,'(A)')  ' NO DUMMY ATOMS ALLOWED BEFORE TRANSLATION'
      write(MIOB,'(A)')  ' ATOM IN A FORCE CALCULATION'
c
      call filclo  !  J/M
      stop ' in READMO -  XX atom before translation one in FORCE!'
c
      end if
  131 continue
      end if
C
C
C OUTPUT FILE TO UNIT 6
C
C    WRITE HEADER
      IDATE=' '
      call DATE(IDATE)
c     write(*,*) ' after DATE, IDATE=',IDATE
      write(MIOB,'(1X,15(''*****''),''****'')')
C
C     CHANGE THE FOLLOWING LINE TO SUIT LOCAL ENVIRONMENT, IF DESIRED
C
      BANNER=
     % ' **              Jagiellonian University                 '
      BANNER(80:80)='*'
      write(MIOB,'(A)')  BANNER
      BANNER=
     % ' ** Department of Computational Methods in Chemistry     '
      BANNER(80:80)='*'
      write(MIOB,'(A)')  BANNER
      BANNER=
     % ' ** adapted for OS/2 system by J.M. A.M.,    July 1993   '
      BANNER(80:80)='*'
      write(MIOB,'(A)')  BANNER
C
C    THE BANNER DOES NOT APPEAR ANYWHERE else.
C
      write(MIOB,'(1X,79(''*''))')
      LINE='   MNDO'
      if (INDEX(KEYWRD,'MINDO') .NE. 0) LINE='MINDO/3'
      if (INDEX(KEYWRD,'AM1')   .NE. 0) LINE='    AM1'
      if (INDEX(KEYWRD,'PM3')   .NE. 0) LINE='    PM3'
      write(MIOB,
     % '(/29X,A,'' CALCULATION RESULTS'',28X,///1X,15(''*****''),
     1          ''****'' )')   LINE(:7)
      write(MIOB,'('' *'',10X,''MOPAC:  VERSION '',F5.2/
     1          '' *'',10X,''CALC''''D. '',A/'' *'')') VERSON, IDATE
C
C CONVERT ANGLES TO RADIANS
      do 140 J=2,3
C$DOIT VBEST
         do 140 I=1,NATOMS
            GEO(J,I) = GEO(J,I) * CONVTR
  140 continue
C
C CHECK DATA
C
      do 150 I=1,NATOMS
         if (LABELS(I) .LE. 0 ) then
            write(MIOB,'('' ATOMIC NUMBER OF '',I3,'' ?'')') LABELS(I)
            if (I.EQ.1) then
               write(MIOB,'(A)')  ' THIS WAS THE FIRST ATOM'
            else
               write(MIOB,'(A)')
     %           '    GEOMETRY UP TO, BUT NOT INCLUDING, THE'//
     %           ' FAULTY ATOM'
               NATOMS=I-1
               call GEOUT(MIOB)
            end if
c
            call filclo  !  J/M
            stop ' in READMO - check input!'
c
         end if
         if (  NA(I).GE.I.OR. NB(I).GE.I.OR. NC(I).GE.I
     1  .OR. (NA(I).EQ.NB(I))   .AND. I.GT.1
     2  .OR. (NA(I).EQ.NC(I).OR.NB(I).EQ.NC(I))  .AND. I.GT.2
     3  .OR.  NA(I)*NB(I)*NC(I).EQ.0  .AND. I.GT.3) then
            write(MIOB,'('' ATOM NUMBER '',I3,'' IS ILL-DEFINED'')') I
c
            if (I.EQ.1) then
              call filclo  !  J/M
              stop ' in READMO - problems with atom #1'
            end if
c
            write(MIOB,'(/,''  GEOMETRY READ IN'',/)')
            call GEOUT(MIOB)
c
            call filclo  !  J/M
            stop ' in READMO - problems with definitions of an atom'
c
         end if
  150 continue
C
C WRITE KEYWORDS BACK TO USER AS FEEDBACK
      call WRTKEY(KEYWRD)
      write(MIOB,'(1X,13(''*****''),''*'',I3.3,''BY'',I3.3)')
     %                                                    MAXHEV,MAXLIT
C
C FILL IN GEO MATRIX IF NEEDED
      if (INDEX(KEYWRD,'OLDGEO').EQ.0.AND.INDEX(KEYWRD,'SYM') .NE. 0
     1   .AND. NDEP.EQ.0) call GETSYM
      if (NDEP.NE.0)       call SYMTRY
C
C INITIALIZE FLAGS FOR OPTIMIZE AND PATH
      IFLAG = 0
      LATOM = 0
      NUMAT=0
      if (NVAR.NE.0) then
         NUMAT=NATOMS
      else
         do 180 I=1,NATOMS
            if (LABELS(I).NE.99.AND.LABELS(I).NE.107) NUMAT=NUMAT+1
            do 180 J=1,3
               if (LOPT(J,I) ) 160, 180, 170
C    FLAG FOR PATH
  160          CONVRT=1.D0
               if ( IFLAG .NE. 0 ) then
                  if (INDEX(KEYWRD,'STEP1').NE.0) then
                     LPARA1=LPARAM
                     LATOM1=LATOM
                     LPARA2=J
                     LATOM2=I
                     LATOM=0
                     IFLAG=0
                     go to 180
                  else
                     write(MIOB,
     %                 '('' ONLY ONE REACTION COORDINATE PERMITTED'')')
                     call filclo  !  J/M
                     stop ' in READMO - more than 1 rection coordinate'
c
                  end if
               end if
               LATOM  = I
               LPARAM = J
               if (J.GT.1) CONVRT=0.01745329252D00
               REACT(1)  = GEO(J,I)
               IREACT=1
               IFLAG = 1
               go to 180
C    FLAG FOR OPTIMIZE
  170          NVAR = NVAR + 1
               LOC(1,NVAR) = I
               LOC(2,NVAR) = J
               XPARAM(NVAR)   = GEO(J,I)
  180    continue
      end if
C READ IN PATH VALUES
      if (IFLAG.EQ.0) go to 220
      if (INDEX(KEYWRD,'NLLSQ').NE.0) then
         write(MIOB,'(A)') ' NLLSQ USED WITH REACTION PATH; '//
     %                     'THIS OPTION IS NOT ALLOWED'
c
         call filclo  ! J/M
         stop ' in READMO - NLLSQ & reaction path are forbidden'
c
      end if
      if (INDEX(KEYWRD,'SIGMA').NE.0) then
         write(MIOB,'(A)') ' SIGMA USED WITH REACTION PATH; '//
     %                     'THIS OPTION IS NOT ALLOWED'
c
         call filclo  ! J/M
         stop ' in READMO - SIGMA & reaction path are forbidden'
c
      end if
      if (INDEX(KEYWRD,'STEP')+INDEX(KEYWRD,'POINTS').NE.0) then
         STEP = READA(KEYWRD,INDEX(KEYWRD,'STEP=')+5)
         NPTS = READA(KEYWRD,INDEX(KEYWRD,'POINT=')+6)
         if (NPTS.GT.200) then
            write(MIOB,
     %        '(///,''    ONLY TWO HUNDRED POINTS ALLOWED IN REACT'',
     %              ''ION COORDINATE'')')
c
            call filclo  ! J/M
            stop ' in READMO - more than 200 points in  reaction path'
c
         end if
         if (LPARAM.EQ.1.AND.STEP.LE.0) then
            write(MIOB,
     %        '(///,''    STEP FOR BOND LENGTH SHOULD BE SET POSITIVE'',
     %              '' TO PREVENT TWO ATOMS COLLAPSE'')')
c
            call filclo  ! J/M
            stop ' in READMO - step < 0 for bond length is illegal'
c
         end if
         go to 220
      end if
  190 read(MIOA,'(A)',end=210) LINE
      call NUCHAR(LINE,VALUE,NREACT)
      if (NREACT.EQ.0) go to 210
      do 200 I=1,NREACT
         IJ=IREACT+I
         if (IJ.GT.200) then
            write(MIOB,
     %        '(///,''    ONLY TWO HUNDRED POINTS ALLOWED IN REACTION'',
     %              '' COORDINATE'')')
c
            call filclo  ! J/M
            stop ' in READMO - more than 200 points in reaction coord.'
c
         end if
         REACT(IJ)=VALUE(I)*CONVRT
         if (ABS(REACT(IJ)-REACT(IJ-1)).LT.1.D-5) then
            DUM1 = REACT(IJ)/CONVRT
            DUM2 = REACT(IJ-1)/CONVRT
            write(MIOB,'(///,'' TWO ADJACENT POINTS ARE IDENTICAL:  '',
     % F7.3,2X,F7.3,/,'' THIS IS NOT ALLOWED IN A PATH CALCULATION'')')
     %                                                        DUM1,DUM2
c
            call filclo  ! J/M
            stop ' in READMO - two adjacent points identical in path'
c
         end if
  200 continue
      IREACT=IREACT+NREACT
      go to 190
  210 continue
      DEGREE=1.D0
      if (LPARAM.GT.1) DEGREE=90.D0/ASIN(1.D0)
      if (IREACT.LE.1) then
         write(MIOB,'(//10X,'' NO POINTS SUPPLIED FOR REACTION PATH'')')
         write(MIOB,'(//10X,'' GEOMETRY AS READ IN IS AS FOLLOWS'')')
         call GEOUT(1)
c
         call filclo  ! J/M
         stop ' in READMO - no data for reaction path'
c
      else
         write(MIOB,'(//10X,'' POINTS ON REACTION COORDINATE'')')
         write(MIOB,'(10X,8F8.2)') (REACT(I)*DEGREE,I=1,IREACT)
      end if
      Iend=IREACT+1
      REACT(Iend)=-1.D12
C
C OUTPUT GEOMETRY AS FEEDBACK
C
  220 call WRTTXT(MIOB)
      if (INDEX(KEYWRD,'NOLOG').EQ.0) then
         open(unit=MIOE, form='FORMATTED', status='UNKNOWN',
     +       file=GETNAM('FOR011'))
         call WRTTXT(MIOE)
      end if
      call GEOUT(1)
      call GMETRY(GEO,COORD)
      if (INDEX(KEYWRD,'NOXYZ') .EQ. 0) then
         if (INDEX(KEYWRD,'0SCF').NE.0) then
C
C  WRITE OUT CARTESIAN COORDINATES FOR USE AS A DATA-SET
C
            write(MIOB,'(A)')
     %        '   GEOMETRY IN CARTESIAN COORDINATE FORMAT'
            call WRTTXT(MIOB)
            J=0
            do 230 I=1,NATOMS
               if (LABELS(I).NE.99) then
                  J=J+1
                  write(MIOB,'(2X,A,3(F19.13,I3))')
     %                        ELEMNT(LABELS(I)),(COORD(K,J),1,K=1,3)
               end if
  230       continue
         else
            write(MIOB,'(//10X,''CARTESIAN COORDINATES '',/)')
            write(MIOB,'(4X,''NO.'',7X,''ATOM'',9X,''X'',
     %                   9X,''Y'',9X,''Z'',/)')
            L=0
            do 240 I=1,NATOMS
               if (LABELS(I) .EQ. 99.OR.LABELS(I).EQ.107) go to 240
               L=L+1
               write(MIOB,'(I6,8X,A2,4X,3F10.4)')
     %                         L,ELEMNT(LABELS(I)),(COORD(J,L),J=1,3)
  240       continue
         end if
      end if
      if (   INDEX(KEYWRD,' XYZ') .NE. 0 ) then
         if ( NVAR .NE. 0 .AND.
     1 INT.AND.(NDEP .NE. 0 .OR.  NVAR.LT.3*NUMAT-6)) then
            if (NDEP.NE.0)
     1write(MIOB,'(//10X,'' INTERNAL COORDINATES READ IN, AND SYMMETRY''
     2 /10X,      '' SPECIFIED, BUT CALCULATION TO BE RUN IN CARTESIAN''
     3           ,'' COORDINATES'')')
            if (NVAR.LT.3*NUMAT-6)
     1        write(MIOB,'(//10X,'' INTERNAL COORDINATES READ IN, AND'',
     2  '' CALCULATION '',/10X,''TO BE RUN IN CARTESIAN COORDINATES, ''
     3       /10X,''BUT NOT ALL COORDINATES MARKED FOR OPTIMISATION'')')
            write(MIOB,
     %        '(//10X,'' THIS INVOLVES A LOGIcallLY ABSURD CHOICE''
     %       /10X,'' SO THE CALCULATION IS TERMINATED AT THIS POINT'')')
c
            call filclo  ! J/M
            stop ' in READMO - all cartesian coords must be optimized'
c
         end if
         SUMX=0.D0
         SUMY=0.D0
         SUMZ=0.D0
         do 250 J=1,NUMAT
            SUMX=SUMX+COORD(1,J)
            SUMY=SUMY+COORD(2,J)
  250    SUMZ=SUMZ+COORD(3,J)
         SUMX=SUMX/NUMAT
         SUMY=SUMY/NUMAT
         SUMZ=SUMZ/NUMAT
         do 260 J=1,NUMAT
            GEO(1,J)=COORD(1,J)-SUMX
            GEO(2,J)=COORD(2,J)-SUMY
  260    GEO(3,J)=COORD(3,J)-SUMZ
         NA(1)=99
         J=0
         NVAR=0
         do 280 I=1,NATOMS
            if (LABELS(I).NE.99) then
               J=J+1
               if (J.EQ.1) then
                  K=0
               else if (J.LT.4) then
                  K=MIN(3,I-1)
               else
                  K=3
               end if
               do 270 L=1,K
                  NVAR=NVAR+1
                  LOC(1,NVAR)=J
                  LOC(2,NVAR)=L
  270          XPARAM(NVAR)=GEO(L,J)
               LABELS(J)=LABELS(I)
            end if
  280    continue
         NATOMS=NUMAT
      else
         if (NVAR.EQ.0) return
         if ( .NOT. INT.AND.(NDEP .NE. 0 .OR.  NVAR.LT.3*NUMAT-6)) then
            if (NDEP.NE.0)
     1        write(MIOB,
     2          '(//10X,'' CARTESIAN COORDINATES READ IN, AND SYMMETRY''
     3       /10X,'' SPECIFIED, BUT CALCULATION TO BE RUN IN INTERNAL ''
     3           ,''COORDINATES'')')
            if (NVAR.LT.3*NUMAT-6)
     1       write(MIOB,'(//10X,'' CARTESIAN COORDINATES READ IN, AND'',
     2   '' CALCULATION '',/10X,''TO BE RUN IN INTERNAL COORDINATES, ''
     3       /10X,''BUT NOT ALL COORDINATES MARKED FOR OPTIMISATION'')')
            write(MIOB,
     1        '(//10X,''MOPAC, BY DEFAULT, USES INTERNAL COORDINATES''
     2    /10X,''TO SPECIFY CARTESIAN COORDINATES USE KEY-WORD :XYZ:'')
     3        ')
            write(MIOB,
     %      '(10X,''YOUR CURRENT CHOICE OF KEY-WORDS INVOLVES'',
     %            '' A LOGIcallLY''/10X,
     %            ''ABSURD CHOICE SO THE CALCULATION IS'',
     %            '' TERMINATED AT THIS POINT'')')
c
            call filclo  ! J/M
            stop ' in READMO - confusing choice of keywords'
c
         end if
      end if
      return
      end
