program xgrad
      implicit none
      integer*4 n,nat,i,maxstep
      real*8 glim,gn,sft,pi,xmax,sp,gno,gm,gmo,rm,dm,slim
      real*8,allocatable::r(:),g(:),go(:)
      integer*4,allocatable::z(:)
      logical lgo,ls
      pi=4.0d0*datan(1.0d0)
      call readopt(glim,maxstep,xmax,slim)
      call initoutput(maxstep)
      inquire(file='GR.SCR',exist=lgo)
      inquire(file='S.SCR',exist=ls)
      if(ls)then
       open(17,file='S.SCR')
       read(17,*)sft
       close(17)
       write(3,603)'Current sft:',sft
 603   format(1x,a12,f12.6)
      else
       sft=1.0d0
       write(3,603)'Initial sft:',sft
      endif
      allocate(z(1),r(1))
      call rg('G98.OUT',0,nat,z,r,.true.)
      deallocate(z,r)
      n=3*nat
      allocate(z(nat),r(n),g(n),go(n))
      call rg('G98.OUT',1,nat,z,r,.true.)
      call wx(r,nat,z)

      call rg('G98.OUT',0,nat,z,r,.true.)
      call rgg(n,g,'G98.OUT',gn,gm)
      write(3,103)gn,gm,glim
103   format(' Norm of gradient:',f15.7, ' Largest element:',f15.7,/,
     1       ' Limit norm grad :',f15.7)
      if(gn.lt.glim)call report('Converged')
      go=0
      if(lgo)call rgg(n,go,'GR.SCR',gno,gmo)
      call wg(n,z,g,'GR.SCR')
      if(lgo)then
       sp=0.0d0
       do 3 i=1,n
3      sp=sp+g(i)*go(i)
       write(3,301)acos(sp/gn/gno)*180.0d0/pi
301    format(' Angle between old and new gradients:',f8.2,' deg.')
       if(sp.lt.0)then
        sft=sft/2.0d0
       else
        if(sft.lt.0.5d0)sft=sft*1.25d0
       endif
      endif
      rm=0.0d0
      do 2 i=1,n
      dm=sft*g(i)*xmax/dabs(gm)
      if(dabs(dm).gt.dabs(rm))rm=dm
2     r(i)=r(i)+dm
      write(3,302)rm,slim
302   format(' Largest shift:',f12.6,' A, limit:',f12.6)
      if(dabs(rm).lt.slim)call report('Stopped for too small shift')
      open(17,file='S.SCR')
      write(17,*)sft
      close(17)
      call ginp(r,nat,z,0)
      end


      subroutine readopt(glim,maxstep,xmax,slim)
      implicit none
      integer*4 maxstep
      real*8 glim,xmax,slim
      logical lex1,lex2
      character*80 ts
c     glim ... limit for gradient norm
c     slim ... limit for coordinate shift
c     xmax ... maximal coordinate shift
c     maxstep ... maximal number of steps
      glim=1.0d-4
      slim=1.0d-4
      xmax=0.1d0
      maxstep=-1
      inquire(file='Q.OPT',exist=lex1)
      inquire(file='X.OPT',exist=lex2)
      if(lex1.or.lex2)then
       if(lex1)open(10,file='Q.OPT')
       if(lex2)open(10,file='X.OPT')
1      read(10,1000,end=999,err=999)ts
1000   format(a80)
       if(ts(1:4).eq.'glim')read(10,*)glim
       if(ts(1:4).eq.'slim')read(10,*)slim
       if(ts(1:4).eq.'xmax')read(10,*)xmax
       if(ts(1:4).eq.'qmax')read(10,*)xmax
       if(ts(1:7).eq.'MAXSTEP')read(10,*)maxstep
       if(ts(1:7).eq.'maxstep')read(10,*)maxstep
       goto 1
999    close(10)
      endif
      return
      end

      subroutine rg(fo,ic,nat,z,r,lzmat)
c     reads geometry from gaussian output
      implicit none
      integer*4 ic,nat,z(*),I,l,ix
      real*8 r(*),q
      character*(*) fo
      character*80 s80
      logical lzmat
      open(2,file=fo)
1     read(2,2000,end=99,err=99)s80
2000  format(a80)
      IF((lzmat.and.(s80(19:39).EQ.'Z-Matrix orientation:'.OR.
     1               s80(26:46).EQ.'Z-Matrix orientation:'.OR.
     1               s80(20:37).EQ.'Input orientation:'.OR.
     1               s80(27:44).EQ.'Input orientation:'))
     1    .OR.
     1    ((.not.lzmat).and.
     2              (s80(20:40).EQ.'Standard orientation:'.OR.
     2               s80(26:46).EQ.'Standard orientation:')))THEN
       DO 2004 I=1,4
2004   READ(2,*)
       if(ic.eq.0)then
        l=0
2005    READ(2,2000)s80
        IF(s80(2:4).NE.'---')THEN
         l=l+1
         BACKSPACE 2
         READ(2,*)q,q
         IF(q.EQ.-1)l=l-1
         GOTO 2005
        ENDIF
        nat=l
       else
        do 2 l=1,nat
   2    read(2,*)z(l),z(l),r(1+3*(l-1)),(r(ix+3*(l-1)),ix=1,3)
        write(3,300)nat
300     format(i6,' atoms read in ',$)
        do 3 i=1,len(fo)
   3    write(3,301)fo(i:i)
301     format(a1,$)
        write(3,*)
       endif
       close(2)
       return
      endif
      goto 1
99    nat=0
      close(2)
      call report('geometry not found')
      end

      subroutine rgg(n,g,fn,gn,gm)
c     reads gradient from fn
      implicit none
      integer*4 n,i,j,nat
      real*8 g(n),gn,gm
      character*80 s80
      character*(*) fn
      nat=n/3
      open(17,file=fn,status='old')
101   read(17,80,end=77,err=77)s80
 80   format(a80)
      if(s80(38:59).eq.'Forces (Hartrees/Bohr)')then
       read(17,*)
       read(17,*)
       do 5 i=1,nat
5      read(17,*)(g(1+3*(i-1)),j=1,3),(g(j+3*(i-1)),j=2,3)
       gn=0.0d0
       gm=0.0d0
       do 6 i=1,n
       if(dabs(g(i)).gt.dabs(gm))gm=g(i)
6      gn=gn+g(i)**2
       gn=dsqrt(gn)
       close(17)
       return
      endif
      goto 101
77    close(17)
      call report('gradient not found in '//fn)
      end

      subroutine wg(n,z,g,fn)
      implicit none
      integer*4 n,i,j,nat,z(*)
      real*8 g(*)
      character*(*) fn
      open(17,file=fn)
      nat=n/3
      write(17,100)
100   format(1x,67(1h-),/,
     1' Center     Atomic',19x,'Forces (Hartrees/Bohr)',/,
     11x,2(6hNumber,5x),9x,1hX,14x,1hY,14x,1hZ,/,67(1h-))
      do 5 i=1,nat
5     write(17,101)i,z(i),(g(j+3*(i-1)),j=1,3)
101   format(i7,i9,7x,3f15.9)
      write(17,102)
102   format(1x,67(1h-))
      close(17)
      return
      end

      subroutine report(r)
      character*(*) r
      write(3,*)r
      close(3)
      stop
      end

      subroutine msg(r)
      character*(*) r
      write(3,*)r
      return
      end
      
      subroutine initoutput(maxstep)
      character*80 ts
      ist=0
      open(3,file='QGRAD.OUT')
2     read(3,3333,end=200,err=200)ts
3333  format(a80)
      if(ts(15:24).eq.'-   Step -')ist=ist+1
      goto 2
200   backspace 3
      write(3,3000)ist
3000  format(/,'   ------------   Step --------------',/,
     1         '   ------------ ',i6,' --------------',/)
      if(maxstep.ge.0.and.ist.ge.maxstep)call report(' Too many steps.')
      return
      end

      subroutine wx(x,nat,ity)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION x(*),ity(*)
c
      write(3,300)
300   format(/,' Cartesian coordinates:',
     1/,25x,'Z-Matrix orientation:',
     2/,1x,69(1h-),
     3/,' Center     Atomic     Atomic',14x,'Coordinates (Angstroms)',
     4/,' Number     Number      Type',14x,'X           Y           Z',
     5/,1x,69(1h-))
      do 1 i=1,nat
1     write(3,3001)i,ity(i),0,(x(3*(i-1)+j),j=1,3)
3001  format(i5,i11,i14,4x,3f12.6)
      write(3,301)
301   format(1x,69(1h-))
      return
      end

      subroutine ginp(x,nat,ity,nq)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      PARAMETER (MENDELEV=89)
      LOGICAL*4 ldalton,notemptyline,lturbo,lgauss,
     1ltinker,lcct
      CHARACTER*80 filename,ts
      CHARACTER*5 s5
      DIMENSION x(*),ity(*),ib(7)
      real*8,allocatable::cq(:),xq(:)
      CHARACTER*2 sy(MENDELEV)
      data sy/' 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'/
c
      allocate(cq(nq),xq(3*nq))
      filename='G98.INP'
      ldalton=.false.
      lcct=.false.
      ltinker=.false.
      lturbo=.false.
      lgauss=.true.
      open(10,file='Q.OPT')
1     read(10,1000,end=999,err=999)ts
1000  format(a80)
      if(ts(1:14).eq.'gaussian input')read(10,1000)filename
      if(ts(1:4).eq.'dalt')read(10,*)ldalton
      if(ts(1:4).eq.'lcct')read(10,*)lcct
      if(ts(1:4).eq.'tink')read(10,*)ltinker
      if(ts(1:4).eq.'turb')read(10,*)lturbo
      goto 1
999   close(10)
c
      if(lcct.or.ldalton.or.lturbo.or.ltinker)lgauss=.false.

      if(ltinker)filename='tinker.xyz'
      if(ldalton)filename='d.mol'
      if(lcct)filename='FILE.X'
      if(lturbo)filename='coord'

      open(10,file=filename)
      open(11,file='INP.NEW')

c     ccccccccccc cct:   cccccccccccccccccccccccccccccccccccc
      if(lcct)then
       read(10,1000)ts
       write(11,1000)ts
       read(10,*)nat
       write(11,*)nat
       do 6 ia=1,nat
       read(10,*)qt,qt,qt,qt,(ib(ii),ii=1,7),qt
6      write(11,1200)ity(ia),(x(3*(ia-1)+ix),ix=1,3),(ib(ii),ii=1,7),qt
1200   format(i5,3f12.6,7i5,f8.4)
      endif
c     cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

C     DDDDDDDDDDD Dalton: DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD
      if(ldalton)then
       read(10,1000,end=899,err=899)ts
       write(11,1000)ts
       iab=1
       if(ts(1:9).eq.'ATOMBASIS')iab=0
       if(iab.eq.1)then
        read(10,1000,end=899,err=899)ts
        write(11,1000)ts
       endif
       do 901 i=1,3
       read(10,1000,end=899,err=899)ts
901    write(11,1000)ts
       ia=0
902    read(10,1000,end=899,err=899)ts
       idp=0
       do 903 i=1,80
903    if(ts(i:i).eq.'.')idp=idp+1
       if(idp.eq.3)then
        ia=ia+1
        write(11,7701)ts(1:3),(x(3*(ia-1)+ix)/0.529177d0,ix=1,3)
7701    format(a3,3f20.15)
        if(ia.lt.nat)goto 902
       else
        write(11,1000)ts
        goto 902
       endif
      endif
C     DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD

C     IIIIIIIIII   Tinker: IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
      if(ltinker)then
       read(10,1000,end=899,err=899)ts
       write(11,1000)ts
       do 905 ia=1,nat
       read(10,1000,end=899,err=899)ts
905    write(11,7702)ts(1:11),(x(3*(ia-1)+ix),ix=1,3),
     1 ts(48:80)
7702   format(a11,3f12.6,a33)
      endif
C     IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

C     TTTTTTTTTTTT Turbomole: TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
      if(lturbo)then
       write(11,1001)
       read(10,*)
1001   format('$coord')
       do 9037 ia=1,nat
       read(10,7741)s5
7741   format(66x,a5)
9037   write(11,7731)(x(3*(ia-1)+ix)/0.529177d0,ix=1,3),s5
7731   format(f20.14,2f22.14,2x,a5)
       write(11,1002)
       read(10,*)
       read(10,*)
1002   format('$user-defined bonds',/,'$end')
      endif
C     TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT

      if(lgauss)then
c     GGGGGGGGGG Gaussian : GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
c
c      rewrite header
       lchk=0
       loniom=0
2      read(10,1000,end=899,err=899)ts
       call chkchek(ts,lchk)
       call chkoniom(ts,loniom)
       if(notemptyline(ts))then
        write(11,1000)ts
        goto 2
       else
        if(lchk.eq.0)write(11,1111)
1111    format('guess=checkpoint')
        write(11,1000)ts
       endif 

c
c      name
21     read(10,1000,end=899,err=899)ts
       write(11,1000)ts
       if(notemptyline(ts))goto 21
c
c      charge and multiplicity
       read(10,1000,end=899,err=899)ts
       write(11,1000)ts
c

c      Cartesian coordinates - replace by new guess
       do 3 ia=1,nat
       if(loniom.eq.1)then
        read(10,1000)ts
        write(11,7700)ts(1:19),(x(3*(ia-1)+ix),ix=1,3),ts(56:80)
7700    format(a19,3f12.6,a25)
       else
        read(10,1000)ts
        write(11,1100)sy(ity(ia)),(x(3*(ia-1)+ix),ix=1,3)
1100    format(a2,2x,3f15.8)
       endif
3      continue
c
c      fix point charges to keep constant potentials at nuclei:
       if(nq.gt.0)then 
        read(10,1000,end=899,err=899)ts
        write(11,1000)ts
        do 4 i=1,nq
4       read(10,*)(xq(3*(i-1)+ix),ix=1,3),cq(i)
        do 5 i=1,nq
5       write(11,119)(xq(3*(i-1)+ix),ix=1,3),cq(i)
119     format(4f12.6)
       endif
c
c      optional parameters until end of file:
31     read(10,1000,end=899,err=899)ts
       write(11,1000)ts
       goto 31
      endif
c     GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
c      
899   close(10)
      close(11)
      call msg('New input written in INP.NEW')
      return
      end

      subroutine chkchek(ts,lchk)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      CHARACTER*(*) ts
      CHARACTER*16 s
c
      l=len(ts)
      do 1 i=1,l-15
      do 2 ii=1,16
      ia=ichar(ts(i+ii-1:i+ii-1))
      if(ia.ge.65.and.ia.le.90)ia=ia+32
2     s(ii:ii)=char(ia)
1     if(s.eq.'guess=checkpoint')lchk=1
      return
      end

      subroutine chkoniom(ts,loniom)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      CHARACTER*(*) ts
      CHARACTER*16 s
c
      l=len(ts)
      do 1 i=1,l-4
      do 2 ii=1,5
      ia=ichar(ts(i+ii-1:i+ii-1))
      if(ia.ge.65.and.ia.le.90)ia=ia+32
2     s(ii:ii)=char(ia)
1     if(s.eq.'oniom')loniom=1
      return
      end

      function notemptyline(ts)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      LOGICAL*4 notemptyline
      CHARACTER*(*) ts
c
      l=len(ts)
      notemptyline=.false.
      do 1 i=1,l
1     if(ts(i:i).ne.' ')notemptyline=.true.
      return
      end