PROGRAM GMOPT
c     optimized geometry from MOPAC output
      implicit real*8(a-h,o-z)
      implicit integer*4 (i-n)
      parameter (nat0=1000,nag0=10)
      DIMENSION r(3,nat0),iz(nat0)
      CHARACTER*50 FIL
      CHARACTER*80 FN,ES
      LOGICAL auto
      CHARACTER*1 oo
      real*8 XX(3),YY(3),A(3),P(3),E0(3),SSE0E0,SSPP,AD,
     1SSAE0,SSPE0,SSXXXX,SSYYYY,SSXXYY,APX,APY,APZ,SIGN,
     2CB,CA,BANGLE,ANGLE,PI,SSE0AP
      dimension itypep(nag0),L0p(nag0),L1p(nag0),L2p(nag0),L3p(nag0)
      logical lex
      CHARACTER*10 fnp(nag0)
      CHARACTER*2 atsy(89),symbol
      data atsy/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne',
     3'Na','Mg','Al','Si',' P',' S','Cl','Ar',
     4' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn',
     4          'Ga','Ge','As','Se','Br','Kr',
     5'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd',
     5          'In','Sn','Sb','Te',' I','Xe',
     6'Cs','Ba','La',
     6               'Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho',
     6               'Er','Tm','Yb','Lu',
     6'Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg',
     6          'Tl','Pb','Bi','Po','At','Rn',
     7'Fr','Ra','Ac'/
     
      nat=0
      iz(1)=0
      ES=' '
      r(1,1)=0.0d0
      inquire(file='AUTO',exist=auto)
      oo='n'
      if(auto)then
       FIL='M.OUT'
      else
       if(iargc().eq.1)then
        call getarg(1,FIL)
       else
        WRITE(*,*)' Full filename of the MOPAC output:'
        READ(*,'(A)')FIL
       endif
      endif

      inquire(file='LIST.PAR',exist=lex)
      if(lex)then
       write(6,*)'LIST.PAR found'
       open(45,file='LIST.PAR')
       open(51,file='LIST.TXT')
       write(6,*)' Appending LIST.TXT'
78781  read(51,*,end=858,err=858)
       goto 78781 
858    iout=45
       read(45,*)nag
       if(nag.gt.nag0)then
        write(6,*)'Too many parameters'
        stop
       endif
       do 21 iag=1,nag
       iout=iout+1
       read(45,4500)fnp(iag)
       write(6,4500)fnp(iag)
4500   format(a10)
       read(45,*)itypep(iag),L0p(iag),L1p(iag),L2p(iag),L3p(iag)
21     write(6,6005)itypep(iag),L0p(iag),L1p(iag),L2p(iag),L3p(iag)
6005   format(5i6)
       close(45)
      endif


      OPEN(4,FILE='FILE.X')
      write(6,*)' Appending FILE.X'
7878  read(4,*,end=888,err=888)
      goto 7878 
888   NG=0

      iy=0
      OPEN(2,FILE=FIL)
      iopt=0

1     READ(2,2000,END=1000)FN
2000  FORMAT(A80)
      if(FN(11:22).eq.'TOTAL ENERGY')then
       read(FN(36:52),*,err=2010)energy
       ES=FN
       ieer=0
2010   continue
      endif

      if(FN(11:33).eq.'FINAL HEAT OF FORMATION')iy=1

      IF(FN(12:29).EQ.'Empirical Formula:')then
       READ(FN(50:55),*)nat
       if(nat.gt.nat0)then
        write(6,*)'too many atoms'
        goto 1000
       endif
      endif

      IF(FN(1:29).EQ.'   ATOM   CHEMICAL          X'.and.iy.eq.1)then
       NG=NG+1
       DO 4 I=1,2
4      READ(2,*)
       do 5 l=1,nat
       READ(2,2000)FN
       symbol=FN(13:14)
       READ(FN(16:34),*)r(1,l)
       READ(FN(38:50),*)r(2,l)
       READ(FN(54:66),*)r(3,l)
       iz(l)=0
       do 6  j=1,89
6      if(atsy(j).eq.symbol)iz(l)=j
       if(iz(l).eq.0)write(6,*)l,symbol,' unknown symbol'
5      continue

       iopt=iopt+1
       write(4,4400)FIL(1:20)//ES(1:60)
       if(.not.auto)write(6,4400)FIL(1:20)//ES(36:60)
       write(6,*)energy
4400   format(A80)
       write(4,*)nat
       do 34 l=1,nat
34     if(iz(l).gt.0)write(4,4000)iz(l),(r(i,l),i=1,3),(0,i=1,7),0.0d0
4000   format(I3,3F12.6,7(1x,i1),f4.1)
       if(lex)then
        do 221 iag=1,nag
        itype=itypep(iag)
        L0=L0p(iag)
        L1=L1p(iag)
        L2=L2p(iag)
        L3=L3p(iag)
        iout=51
        PI=3.141592653589D0
        A(1)=r(1,L2)-r(1,L3)
        A(2)=r(2,L2)-r(2,L3)
        A(3)=r(3,L2)-r(3,L3)
        P(1)=r(1,L0)-r(1,L1)
        P(2)=r(2,L0)-r(2,L1)
        P(3)=r(3,L0)-r(3,L1)
        E0(1)=r(1,L1)-r(1,L2)
        E0(2)=r(2,L1)-r(2,L2)
        E0(3)=r(3,L1)-r(3,L2)
        SSE0E0=E0(1)*E0(1)+E0(2)*E0(2)+E0(3)*E0(3)
        DO 8 KK=1,3
8       E0(KK)=E0(KK)/DSQRT(SSE0E0)
        SSPP=P(1)*P(1)+P(2)*P(2)+P(3)*P(3)
        AD=DSQRT(SSPP)
        SSAE0=A(1)*E0(1)+A(2)*E0(2)+A(3)*E0(3)
        SSPE0=P(1)*E0(1)+P(2)*E0(2)+P(3)*E0(3)
        DO 9 KK=1,3
        XX(KK)=A(KK)-E0(KK)*SSAE0
9       YY(KK)=P(KK)-E0(KK)*SSPE0
        SSXXXX=XX(1)*XX(1)+XX(2)*XX(2)+XX(3)*XX(3)
        SSYYYY=YY(1)*YY(1)+YY(2)*YY(2)+YY(3)*YY(3)
        SSXXYY=XX(1)*YY(1)+XX(2)*YY(2)+XX(3)*YY(3)
        APX=XX(2)*YY(3)-XX(3)*YY(2)
        APY=XX(3)*YY(1)-XX(1)*YY(3)
        APZ=XX(1)*YY(2)-XX(2)*YY(1)
        SSE0AP=APX*E0(1)+APY*E0(2)+APZ*E0(3)
        SIGN=1.0D0
        IF (SSE0AP.NE.0.0D0)SIGN=-SSE0AP/ABS(SSE0AP)
        CB=0.0D0
        IF (AD.NE.0.0D0)CB=-SSPE0/AD
        CA=0.0D0
        IF (SSXXXX.NE.0.0D0.AND.SSYYYY.NE.0.0D0)
     1  CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY)
        BANGLE=0.0D0
        ANGLE=0.0D0
        IF (CB.EQ.-1.0D0)BANGLE=PI
        IF (ABS(CB).LT.1.0D0)BANGLE=ACOS(CB)
        IF (CA.LE.-1.0D0) ANGLE=PI
        IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA)
        ANGLE=180.0D0* ANGLE/PI
        BANGLE=180.0D0*BANGLE/PI
        ANGLE=SIGN*ANGLE
        if(ANGLE.lt.-180.0d0)ANGLE=360.0d0+ANGLE
c       torsion, distance, band angle
        if(itype.eq.2)WRITE(iout,322)AD
        if(itype.eq.3)WRITE(iout,322)BANGLE
221     if(itype.eq.4)WRITE(iout,322)ANGLE
322     format(f15.6,$)
        if(ieer.eq.0)then
         WRITE(iout,3221)energy
3221     format(f20.8)
        else
         WRITE(iout,3222)ES(20:80)
3222     format(a20)
        endif
        goto 1000
       endif
       goto 1000

      ENDIF

      GOTO 1

1000  CLOSE(2)
      if(lex)close(51)
      if(.not.auto)then
       if(ng.eq.0)then
        write(*,*)'Geometry not found'
       else
        WRITE(*,*)' File FILE.X written'
        write(*,*)ng,' geometries found'
        write(*,*)iopt,' optimized'
        write(*,*)nat,' atoms'
       endif
      endif

      STOP
      END