program gdiis
      implicit real*8 (a-h,o-z)
      logical fail
      parameter(iop=10,nqlim=180,nmax=60)
      dimension x(180),g(180),lopt(iop),wk(iop)
      character*4 wk,wo
      common /work/ hlast(180,180),b(32400),hh(32400),fi(180),ll(180)
     1,mmm(180),qi(180),c(60)
      common/work1/ xx(180,60), ee(180,60), gg(180,60),qq(180,60)
      common /optim/h(16290),e(180),yk(180),pk(180),rkl, e11e(3),
     *tk(180),w(180),d1(180,3), don(60,3), at(3,60),
     *ix(60),ipa(180,2),id(1770),id1,ie(60,180),iee(180),itab(180),
     *nk,nat, mk(60),iblanc,qp(3,5400), ak(3,3,60)
      common/optio/ option, opnclo, cndo, indo,mind,closed,open
     *, stndr,optim,idummy,ecrit,rho,ih,iup,idiag,kell
      common/sumary/en(20),gsum(20),isum,ismax
      common/ite/iter,intnl,fail,iix
      common /files/ inp,iout,ix9,ipun,nshout
      common/info/xc(60,3),natoms,ica,mult,an(60),nnn,xm(180),
     1jfix(180),nfix,symb(60)
      integer option,opnclo,cndo,closed,open,stndr,optim,an,symb
      integer gdiisf
      parameter(tol=1.d-8,limit=20)
      data wk/'mole','stop','gdii','tol ','bmat','m&s ','k   ','fixc',
     x 'fmat','card'/
      data zero,one/0.0d0,1.0d0/
c
c     optimize geometry by diis algorithm
c
c     on input
c
c     x  array of variables (geometry parameters ie. cartesians
c        or internal coordinates, the latter is preferred),dimens.=n
c
c     f   function value  (i.e. energy)
c
c     g  array of gradients(dim.=n)
c
c     eps  convergence limit .
c        exit if predicted changes in variables become less than
c        eps, (a suggested value is 1.d-4)
c
c     limit  allowed number of iterational steps (limit.lt.20)
c
c     ih=0  start by unit matrix for the hessian (default value)
c        1  start by a diagonal hessian
c        2  start by a hessian readed  from cards
c        the hessian must be stored in array h( of common/optim/ )
c        before call gdiis. if (ih.le.0) inversion of h not carried out.
c
c     iup = 0   update of h-1 not requiered (default value)
c     iup = 1   update of h-1 according murtagh&sargent
c
c
c      open files
c
      inp=5
      iout=6
      open (unit=inp,file='gdiisin',status='old')
      open (unit=iout,file='gdiisls')
c       file shout short output
      nshout=13
      open(unit=nshout,file='gdiism',status='unknown',
     &access='sequential',form='formatted')
      rewind nshout
c        file geomfl geometry file
      inp2=31
      open(unit=inp2,file='geom',status='old')
c        file cartgrd cartesian gradient
      inp1=32
      open(unit=inp1,file='cartgrd',status='old')
c        file gdiisf  gdiis information
      gdiisf=33
      open(unit=gdiisf,file='gdiisfl',status='unknown',
     $     form='unformatted')
c        file scrfl   for internal coordinate system
      iix=34
      open(unit=iix,file='scrfl',status='unknown')
      rewind gdiisf
      read(gdiisf,end=3100) m
 3000 format(i8)
      if(m.gt.1)goto 3900
 3100 m=1
      call setup(ih,lopt,iop,nfix,jfix,kell,eps,natp,m,nqlim,inp
     x ,iout,xita,wk,iix)
c
c     compute gradient at start point
c
      call objet(xc,en(1),g,natoms,inp2,iout,natp,nmax,inp1,nshout,
     x niter,symb,xm)
      ij=0
      idsp=0
      do 5 k=1,natoms
      do 5 i=1,3
      ij=ij+1
      x(ij)=xc(k,i)*0.52917706d0
      xx(ij,1)=x(ij)
    5 g(ij)=g(ij)*8.238856d0
      call conver(x,g,natoms,qi,fi,nq,xm,idsp,ncard)
      do 15 i=1,nq
      qq(i,1)=qi(i)
   15 gg(i,1)=fi(i)
      call dzero(nq*nq,hh,1)
      if(ih-1) 20,50,80
   20 ji=0
      do 40 i=1,nq
      ii=(i-1)*nq+i
   40 hh(ii)=1.0d0
      goto 100
   50 ji=0
      read(inp2,501)(e(i),i=1,nq)
  501 format(8f10.7)
      do 70 i=1,nq
      ii=(i-1)*nq+i
   70 hh(ii)=one/e(i)
      goto 100
   80 continue
      do 85 i=1,nq
      read(inp2,501)(h(j),j=1,i)
      do 85 j=1,i
      ij=(i-1)*nq+j
      ji=(j-1)*nq+i
      hh(ij)=h(j)
      hh(ji)=h(j)
      if(nfix.eq.0.or.i.eq.j)goto85
      do86k=1,nfix
      ii=jfix(k)
      if(i.ne.ii.and.j.ne.ii)goto86
      hh(ij)=zero
      hh(ji)=zero
  86  continue
  85  continue
      call osinv(hh,nq,det,tol,ll,mmm)
      if (det.gt.zero) goto 100
      write (iout,90)m,det
      write (nshout,90)m,det
   90 format (1x,7hnstep =, i3,26hhessian is singular  det =, e13.6 )
      stop ' hessian is singular'
  100 continue
      if(nfix.eq.0)goto105
      ij=0
      do 106 i=1,nq
      do 106 j=1,nq
      ij=ij+1
      do 106 k=1,nfix
      ii=jfix(k)
      if(i.eq.ii.or.j.eq.ii)hh(ij)=zero
  106 continue
  105 continue
      call relax(hh,fi,qi,nq,e)
      do 120 i=1,nq
      ee(i,1)=e(i)
  120 continue
      call displ(x,qi,natoms,nq,xm,symb,ncard,e)
      ij=0
      do 110 k=1,natoms
      do 110 i=1,3
      ij=ij+1
  110 x(ij)=x(ij)/0.52917706d0
      niter=niter+1
      m=2
      if(kell.le.0)kell=nq+1
      if(kell.gt.limit)kell=limit
      call dump(m,kell,nq,natoms,niter,xita,x,natp,lopt,gdiisf,
     x hh,inp2,iout,xm,nshout)
      stop
 3900 continue
      call setup(ih,lopt,iop,nfix,jfix,kell,eps,natp,m,nqlim,inp
     x ,iout,xita,wk,iix)
      call rest(m,kell,lang,nq,xita,lopt,gdiisf,hh,natp)
      call objet(xc,en(lang),g,natoms,inp2,iout,natp,nmax,inp1,nshout,
     x niter,symb,xm)
      ij=0
      do 115 k=1,natoms
      do 115 i=1,3
      ij=ij+1
      x(ij)=xc(k,i)*0.52917706d0
      xx(ij,lang)=x(ij)
  115 g(ij)=g(ij)*8.238856d0
      call conver(x,g,natoms,qi,fi,nq,xm,idsp,ncard)
      do 121 i=1,nq
      qq(i,lang) = qi(i)
      gg(i,lang) = fi(i)
  121 continue
      call relax(hh,fi,qi,nq,e)
      do 122 i=1,nq
      ee(i,lang) = e(i)
  122 continue
      if(m.le.2) goto 600
c
c   update hessian if required
c
      if (iup.eq.0) goto 400
      do 300 i=1,nq
  300 yk(i)=g(i)-yk(i)
 1000 continue
c     murtagh&sargent
      xita = xita/dfloat(m)
      do 310 i=1,nq
  310 e(i) = -xita*e(i) - updot(h,yk,i,nq)
      sg = dot(e,yk,nq)
      shs = dot(e,e,nq)
      shg = dot(e,g,nq)
      roshs = 1.0d-3*shs
      if (dabs(sg) - roshs) 320,320,330
  320 sg = dsign(roshs,sg)
  330 if (shg/sg+1.0d-6) 340,340,360
  340 ji=0
      do 350 i=1,nq
      do 350 j=1,nq
      ji = ji + 1
  350 hh(ji) = hh(ji) + e(i)*e(j)/sg
  360 continue
c   update error vectors
      do 390 j=1,m
      do 370 i=1,nq
  370 fi(i)=gg(i,j)
      call relax(hh,fi,qi,nq,e)
      do 380 i=1,nq
  380 ee(i,j)=e(i)
  390 continue
  400 continue
c
c   exit if max(abs(e(i)).lt.eps.or.m.gt.limit
c
      s=zero
      do 440 i=1,nq
      if (dabs(e(i)).gt.s) s=dabs(e(i))
  440 continue
      if(s.lt.eps) goto 500
      if(m+1.gt.limit) goto 700
      goto 600
c...  exit
  500 continue
      write(iout,1011) m
      write(nshout,1011) m
 1011 format(1x,20(4h****)/1x,23hgeometry optimized in   ,i4,6h-steps,/)
c
c     change r1g in jcl to 1 and terminate iterations
c
c      call cmscmd('globalv set     r1g     1               ',irt)
      stop ' convergent in gdiis'
  700 continue
      stop ' limit is excited '
  600 continue
c
c  start diis
c
      m1=nq+1
      i1=1
      mm=lang
      iii=0
c&&
c     print*,' kell=',kell,' m1=',m1,' lang=',lang,' i1=',i1,' mm',mm
c     do 133 j=i1,mm
c     write(6,132)' ee',(ee(i,j),i=1,nq)
c     write(6,132)' gg',(gg(i,j),i=1,nq)
c     write(6,132)' qq',(qq(i,j),i=1,nq)
c 132 format(a,/,(10f10.6))
c 133 continue
c&&
      do 175 i=i1,mm
      iii=iii+1
      do 170 j=1,nq
      ji=(iii-1)*m1+j
      iiei=ee(j,i)*1.d5
      b(ji)=iiei*1.d-5
  170 e(j)=zero
      ij=iii*m1
      if(i1.eq.1)ij=i*m1
  175 b(ij)=one
      e(m1)=one
c
c   pseudoinversion with linlsq subroutine
c
      eps1=1.d-12
      call linlsq(b,c,e,m1,mm,eps1,ier,hlast,yk,h,pk,qi,ll,w)
      if(ier.eq.0)goto 190
      write(iout,9876)ier
9876  format(/,5x,'ier erteke a linlsq-ban',i5)
      stop ' after linlsq'
c
c   diis coefficients and the lagrangian
c
  190 continue
      s=zero
      do 191 i=1,mm
  191 s=s+c(i)
      do 192 i=1,mm
  192 c(i)=c(i)/s
      write (iout,201)
      write (nshout,201)
  201 format (/, 2x,'the c values ')
      write(iout,203) (c(i),i=1,mm)
      write(nshout,203) (c(i),i=1,mm)
  203 format(1(1x,10f10.6))
c
c   relaxation
c
      do 200 i=1,nq
      qi(i)= zero
      yk(i)=fi(i)
  200 fi(i)=zero
      ii1=0
      do 220 i=i1,mm
      ii1=ii1+1
      do 210 j=1,nq
      fi(j)=fi(j)+idint(c(ii1)*gg(j,i)*1.d6)*1.d-6
  210 qi(j)=qi(j)+idint(c(ii1)*qq(j,i)*1.d6)*1.d-6
  220 continue
      call relax(hh,fi,qi,nq,e)
      do 225 i=1,nq
      e(i)=qi(i)-qq(i,mm)
  225 continue
      call displ(x,qi,natoms,nq,xm,symb,ncard,e)
      ij=0
      do 230 k=1,natoms
      do 230 i=1,3
      ij=ij+1
  230 x(ij)=x(ij)/0.52917706d0
      niter=niter+1
      m=m+1
      call dump(m,kell,nq,natoms,niter,xita,x,natp,lopt,gdiisf,
     x hh,inp2,iout,xm,nshout)
      stop
      end
      subroutine linlsq(a,x,b,m,n,eta,ier,qr,alpha,e,y,r,piv,z)
      implicit real*8(a-h,o-z)
      integer piv
      dimension a(m,n),x(n),b(m)
      dimension qr(m,n),alpha(n),e(n),y(n),r(m),piv(n),z(n)
      data zero/0.d0/,var/0.0625d0/
      ier=0
      lstep=0
      do 10 j=1,n
      do 10 i=1,m
  10  qr(i,j)=a(i,j)
      call decomp(m,n,qr,alpha,e,y,piv,ier)
      if (ier.eq.0) goto 15
      write(6,14)ier
  14  format(/,5x,'matrix of coeff. is singular',5x,'ier =',i3)
      stop
  15  eta2=eta*eta
      do 20 i=1,m
  20  r(i)=b(i)
      call solve(m,n,qr,alpha,piv,r,y,z)
      do 40 i=1,m
      ss=-b(i)
      do 30 j=1,n
      ai=a(i,j)
      bi=y(j)
  30  ss=ss+ai*bi
  40  r(i)=-ss
      call solve(m,n,qr,alpha,piv,r,e,z)
      y0=zero
      e1=zero
      lstep=0
      do 50 i=1,n
      y0=y0+y(i)*y(i)
  50  e1=e1+e(i)*e(i)
      if(e1-y0*var) 60,60,55
  55  ier=2
      write(6,16)ier
  16  format(/5x,'no improvement in linlsq',5x,'ier =',i3)
      return
c     no attempt to improve the solution is made unless the norm of
c     the first correction is significantly smaller than the norm of
c     the initial solution
c
  60  continue
      lstep=lstep+1
      if(lstep.gt.10)goto110
c     improve solution
      do 70 i=1,n
  70  y(i)=y(i)+e(i)
      if(e1-eta2*y0) 110,80,80
  80  do 95 i=1,m
      ss=-b(i)
      do 90 j=1,n
      ai=a(i,j)
      bi=y(j)
  90  ss=ss+ai*bi
  95  r(i)=-ss
      call solve(m,n,qr,alpha,piv,r,e,z)
      e0=e1
      e1=zero
      do 100 i=1,n
 100  e1=e1+e(i)*e(i)
      if(e1-e0*var) 60,60,110
 110  continue
c
c     terminate if either the correction was of little significance,
c     or the norm of correction failed to decrease sufficiently as
c     compared to the prvious correction
c
      do 120 i=1,n
 120  x(i)=y(i)
      if(lstep.ne.1) write(6,130)lstep
 130  format(41x,'a javitashoz felhasznalt lepesszam:',i3)
      return
      end
      subroutine solve(m,n,qr,alpha,piv,r,y,z)
      implicit real*8(a-h,o-z)
      integer piv
      dimension qr(m,n),r(m),alpha(n),y(n),piv(n),z(n)
      data zero/0.d0/
      do 30 j=1,n
      ss=zero
      do 10 i=j,m
      ai=qr(i,j)
      bi=r(i)
  10  ss=ss+ai*bi
      gam=ss/(alpha(j)*qr(j,j))
      do 20 i=j,m
  20  r(i)=r(i)+gam*qr(i,j)
  30  continue
      z(n)=r(n)/alpha(n)
      n1=n-1
      do 50 ii=1,n1
      i=n1-ii+1
      i1=i+1
      ss=-r(i)
      do 40 j=i1,n
      ai=qr(i,j)
      bi=z(j)
  40  ss=ss+ai*bi
  50  z(i)=-ss/alpha(i)
      do 60 i=1,n
      ip=piv(i)
  60  y(ip)=z(i)
      return
      end
      subroutine decomp(m,n,qr,alpha,sum,y,piv,ier)
      implicit real*8(a-h,o-z)
      integer piv
      dimension qr(m,n),alpha(n),piv(n),sum(n),y(n)
      data zero/0.0d0/
      do 20 j=1,n
      ss=zero
      do 10 i=1,m
      ai=qr(i,j)
  10  ss=ss+ai*ai
      sum(j)=ss
  20  piv(j)=j
      do 190 k=1,n
      sigma=sum(k)
      jbar=k
      if(k-n) 30,50,50
  30  k1=k+1
      do 45 j=k1,n
      if(sigma-sum(j)) 40,45,45
  40  sigma=sum(j)
      jbar=j
  45  continue
  50  if(jbar-k) 60,80,60
  60  i=piv(k)
      piv(k)=piv(jbar)
      piv(jbar)=i
      sum(jbar)=sum(k)
      sum(k)=sigma
      do 70 i=1,m
      sigma=qr(i,k)
      qr(i,k)=qr(i,jbar)
  70  qr(i,jbar)=sigma
  80  continue
      ss=zero
      do 90 i=k,m
      ai=qr(i,k)
  90  ss=ss+ai*ai
      sigma=ss
      if(sigma) 110,100,110
 100  ier=1
c     a is singular
      return
 110  qrkk=qr(k,k)
      ss=dsqrt(sigma)
      if(qrkk) 120,130,130
 120  alphak=ss
      alpha(k)=alphak
      goto 140
 130  alphak=-ss
      alpha(k)=alphak
 140  beta=1.d0/(sigma-qrkk*alphak)
      qr(k,k)=qrkk-alphak
      k1=k+1
      if(k1-n) 145,145,165
 145  do 160 j=k1,n
      ss=zero
      do 150 i=k,m
      ai=qr(i,k)
      bi=qr(i,j)
 150  ss=ss+ai*bi
 160  y(j)=beta*ss
 165  if(k1-n) 166,166,185
 166  do 180 j=k1,n
      do 170 i=k,m
 170  qr(i,j)=qr(i,j)-qr(i,k)*y(j)
 180  sum(j)=sum(j)-qr(k,j)*qr(k,j)
 185  continue
 190  continue
      return
      end
      subroutine norm(u)
      implicit real*8(a-h,o-z)
      dimension u(3)
      x=1.0000000000/ dsqrt(scalar(u,u))
      do 1 i=1,3
      u(i)=u(i)*x
    1 continue
      return
      end
      function s2(x)
      implicit real*8(a-h,o-z)
      s2=dsqrt(1.0-x**2)
      return
      end
      function scalar(u,v)
      implicit real*8(a-h,o-z)
      dimension u(3),v(3)
      scalar=0.0
      do 1 i=1,3
      scalar=scalar+u(i)*v(i)
    1 continue
      return
      end
      subroutine normal(u,v,w)
      implicit real*8(a-h,o-z)
      dimension u(3),v(3),w(3)
      w(1)=u(2)*v(3)-u(3)*v(2)
      w(2)=u(3)*v(1)-u(1)*v(3)
      w(3)=u(1)*v(2)-u(2)*v(1)
      call norm(w)
      return
      end
      subroutine machb(nek,bmat,ncmax,nqmax,xa,ya,za,qq,n,inp,iout,nq,
     1 wri,qonly,ncard)
      implicit real*8(a-h,o-z)
      logical wri,qonly
c        nek=3*number of atoms input parameter
c        bmat is the transpose of the b matrix, dimensioned im the
c        call ing program as b(ncmax,nqmax) - output
c        ncmax, nqmax - input
c        xa,ya,za are the x,y,z coordinates of the atoms in angstrom
c        qq contains the values of the internal coordinates (output)
c        inp and iout are the input and output files
c        nq is the number of internal coordinates, wri is .true. if the
c        definition of internal coordinates is to be o printed
c        qonly is true if only coordinates (no b matrix) is  are to be c
c        calculated
c                                   input data
c        each elementary valence coordinate on a separate card
c        col. 1 'k' or' ' (blank). if 'k' a new coordinate begins, if
c        blank then the composite internal coordimnate begun earlier is
c         continued. any other character terminates the input
c        cols. 2-9 scale factor for the total coordinate (only if there
c        is k in column 1. blank or zero is interpreted as 1.0
c        cols. 21-24 coordinate type stre,bend,out ,tors,lin1,lin2
c                    invr,tor2
c        cols. 31-40,41-50,51-60,61-70 participating atoms a,b,c,d
c         format 4f10.x
c         a and b are given for 'stre' order arbitrary
c        a,b,c for bend - a and b are end atoms, c is the apex atom a-c-
c        atom a out of the bcd plane - c is the central atom - coordina
c        te positive if a is displaced toward the vector product db*dc
c        torsion a-b-c-d, positive as in the wilson-decius-cross book
c        note that the value of the coordinate varies between -pi/2 to
c        3pi/2   n o t  between -pi/2 to +pi/2
c        lin1 l  collinear bending a-b-c distorted in the planr of abd
c        positive if a moves toward d
c         lin2 linear bending a-c-b distorted perpendicular to the plane
c        abd - positive if a moves toward the vector cross product cd*ca
c        the linear bendings are a-c-b, i. e. the apex atom is third
      dimension qq(1)
      dimension tipus(8)
      dimension wort(3)
      dimension a(4),ia(4),u(3),v(3),w(3),z(3),x(3),uu(3),vv(3),ww(3),zz
     1 (3),uv(12)
      equivalence(ka,ia(1)),(kb,ia(2)),(kc,ia(3)),(kd,ia(4))
      equivalence (uv(1),uu(1)),(uv(4),vv(1)),(uv(7),ww(1)),(uv(10),zz(1
     1 ))
      dimension bmat(ncmax,nqmax),xa(1),ya(1),za(1)
      integer tipus,typ,we,wort,blank,tlast
      data tipus/4hstre,4hinvr,4hbend,4hout ,4htors,4hlin1,4hlin2,4htor2
     1/
      data wort/1hk,1h ,1h*/
      data anull/1.0d0/
      data blank/4h     /
      ncard=0
      o=1.0d0
      nab=nek
      nab=nab/3.0d0+0.1d0
      i=0
      c1=0.0d0
      ccx=0.0d0
      tlast=blank
      if(wri) write(iout,1901)
 1901 format(/,1x,34hdefinition of internal coordinates   ,/)
  311 read(inp,1101) we
 1101 format(a1)
      backspace inp
      do 3301 k=1,2
      if(we.eq.wort(k)) goto 3302
 3301 continue
      c1=dsqrt(1.0d0/c1)*ccx
      qq(i)=qq(i)*c1
      if(qonly) goto 310
      do 6601 k=1,nek
      bmat(k,i)=bmat(k,i)*c1
 6601 continue
      go to 310
 3302 continue
      ncard=ncard+1
      read(inp,1100) we,cc,c,typ,a
 1100 format(a1,f9.5,f10.4,a4,6x,5f10.4)
      if(typ.eq.blank) typ=tlast
      tlast=typ
      if(cc.eq.0.0d0) cc=1.0d0
      if(we.eq.wort(1)) ccc=ccx
      if(we.eq.wort(1)) ccx=cc
      if(c.eq.0.0d0) c=1.0d0
      if(we.eq.wort(2)) c1=c1+c**2
      if(we.ne.wort(1)) goto 9109
      if(i.eq.0) goto 9107
      if(wri) write(iout,5300)
 5300 format(1x)
      c1=dsqrt(1.0d0/c1)*ccc
      if(qonly) goto 6610
      do 6600 k=1,nek
      bmat(k,i)=bmat(k,i)*c1
 6600 continue
 6610 qq(i)=qq(i)*c1
 9107 i=i+1
      qq(i)=0.0d0
      c1=c**2
      if(qonly) goto 9109
      do 107 j=1,nek
      bmat(j,i)=0.0d0
  107 continue
 9109 do 108 k=1,4
      ia(k)=a(k)+0.1d0
  108 continue
      do 109 k=1,8
      if(typ.eq.tipus(k)) goto 15
  109 continue
      error=7
      write(iout,1902) i
 1902 format(/,1x,38hundefined int.coordinate type at no.  ,i3,/,10(4h**
     1**))
      goto 1111
   15 if(wri) write(iout,7000) i,typ,ia,c,ccx
 7000 format(1x,i3,1h.,a8,4i10,f12.5,f12.6)
      if(ka.lt.1.or.ka.gt.nab.or.kb.lt.1  .or.kb.gt.nab)goto 1310
      if(k.gt.2.and.(   kc.lt.1.or.kc.gt.nab))goto 1310
      if(k.gt.3.and.(kd.lt.1.or.kd.gt.nab))goto 1310
      goto(1,7,2,3,4,5,6,8),k
c..... stretch
    1 call vektor(uu,r1,ka,kb,xa,ya,za)
      uu(1)=uu(1)*anull
      uu(2)=uu(2)*anull
      uu(3)=uu(3)*anull
      vv(1)=-uu(1)
      vv(2)=-uu(2)
      vv(3)=-uu(3)
      ia(3)=0
      ia(4)=0
      qq(i)=qq(i)+r1*c
      goto 1500
c...  inverse
    7 continue
      call vektor(uu,r1,ka,kb,xa,ya,za)
      rm1=1.d0/r1
      rm2=rm1*rm1
      uu(1)=-rm2*uu(1)*anull
      uu(2)=-rm2*uu(2)*anull
      uu(3)=-rm2*uu(3)*anull
      vv(1)=-uu(1)
      vv(2)=-uu(2)
      vv(3)=-uu(3)
      ia(3)=0
      ia(4)=0
      qq(i)=qq(i)+rm1*c
      goto 1500
c.....bending
    2 call vektor(u,r1,ka,kc,xa,ya,za)
      call vektor(v,r2,kb,kc,xa,ya,za)
      co=scalar(u,v)
      si=s2(co)
      do 420 l=1,3
      uu(l)=(co*u(l)-v(l))/(si*r1)
      vv(l)=(co*v(l)-u(l))/(si*r2)
      ww(l)=-uu(l)-vv(l)
  420 continue
      ia(4)=0
      qq(i)=qq(i)+c*darcos(co)
      goto 1500
c.....out of plane
    3 call vektor(u,r1,ka,kd,xa,ya,za)
      call vektor(v,r2,kb,kd,xa,ya,za)
      call vektor(w,r3,kc,kd,xa,ya,za)
      call normal(v,w,z)
      steta=scalar(u,z)
      cteta=s2(steta)
      cfi1=scalar(v,w)
      sfi1=s2(cfi1)
      cfi2=scalar(w,u)
      cfi3=scalar(v,u)
      den=cteta*sfi1**2
      st2=(cfi1*cfi2-cfi3)/(r2*den)
      st3=(cfi1*cfi3-cfi2)/(r3*den)
      do 430 l=1,3
      vv(l)=z(l)*st2
      ww(l)=z(l)*st3
  430 continue
      call normal(z,u,x)
      call normal(u,x,z)
      do 431 l=1,3
      uu(l)=z(l)/r1
      zz(l)=-uu(l)-vv(l)-ww(l)
  431 continue
      cx=-c
      if(steta.lt.0.0) cx=c
      qq(i)=qq(i)-cx*darcos(cteta)
      goto 1500
c..... torsion
    4 call vektor(u,r1,ka,kb,xa,ya,za)
      call vektor(v,r2,kc,kb,xa,ya,za)
      call vektor(w,r3,kc,kd,xa,ya,za)
      call normal(u,v,z)
      call normal(w,v,x)
      co=scalar(u,v)
      co2=scalar(v,w)
      si=s2(co)
      si2=s2(co2)
      do 440 l=1,3
      uu(l)=z(l)/(r1*si)
      zz(l)=x(l)/(r3*si2)
      vv(l)=(r1*co/r2-1.00000000d0)*uu(l)-r3*co2/r2*zz(l)
      ww(l)=-uu(l)-vv(l)-zz(l)
  440 continue
      co=scalar(z,x)
      u(1)=z(2)*x(3)-z(3)*x(2)
      u(2)=z(3)*x(1)-z(1)*x(3)
      u(3)=z(1)*x(2)-z(2)*x(1)
      co2=scalar(u,v)
      s=darcos(-co)
      if(co2.lt.0.0) s=-s
      if(s.gt.1.5707963272d0) s=s-6.2831853072d0
      qq(i)=qq(i)-c*s
c.... remember that the range of this coordinate is -pi/2 to 3*pi/2
c.... in order to shift the discontinuity off the planar position
      goto 1500
c..... torsion2  (only other range as in tors)
    8 call vektor(u,r1,ka,kb,xa,ya,za)
      call vektor(v,r2,kc,kb,xa,ya,za)
      call vektor(w,r3,kc,kd,xa,ya,za)
      call normal(u,v,z)
      call normal(w,v,x)
      co=scalar(u,v)
      co2=scalar(v,w)
      si=s2(co)
      si2=s2(co2)
      do 880 l=1,3
      uu(l)=z(l)/(r1*si)
      zz(l)=x(l)/(r3*si2)
      vv(l)=(r1*co/r2-1.00000000d0)*uu(l)-r3*co2/r2*zz(l)
      ww(l)=-uu(l)-vv(l)-zz(l)
  880 continue
      co=scalar(z,x)
      u(1)=z(2)*x(3)-z(3)*x(2)
      u(2)=z(3)*x(1)-z(1)*x(3)
      u(3)=z(1)*x(2)-z(2)*x(1)
      co2=scalar(u,v)
      s=darcos(-co)
      if(co2.lt.0.0) s=-s
      if(s.gt.           0d0) s=s-6.2831853072d0
      qq(i)=qq(i)-c*s
c.... remember that the range of this coordinate is 0 to pi/2
c.... use it for perpendicular geometries
      goto 1500
c......linear complanar bending
    5 call vektor(u,r1,ka,kc,xa,ya,za)
      call vektor(v,r2,kd,kc,xa,ya,za)
      call vektor(x,r2,kb,kc,xa,ya,za)
      co=scalar(v,u)
      co2=scalar(x,v)
      qq(i)=qq(i)+c*(3.1415926536-darcos(co)-darcos(co2))
      call normal(v,u,w)
      call normal(u,w,z)
      call normal(x,v,w)
      call normal(w,x,u)
c..... coordinate positiv if atom a moves towards atom d
      do 450 l=1,3
      uu(l)=z(l)/r1
      vv(l)=u(l)/r2
      ww(l)=-uu(l)-vv(l)
  450 continue
      ia(4)=0
      goto 1500
c.....linear perpendicular bending
    6 call vektor(u,r1,ka,kc,xa,ya,za)
      call vektor(v,r2,kd,kc,xa,ya,za)
      call vektor(z,r2,kb,kc,xa,ya,za)
      call normal(v,u,w)
      call normal(z,v,x)
      do 460 l=1,3
      uu(l)=w(l)/r1
      vv(l)=x(l)/r2
      ww(l)=-uu(l)-vv(l)
  460 continue
      ia(4)=0
      co=scalar(u,w)
      co2=scalar(z,w)
      qq(i)=qq(i)+c*(3.1415926536-darcos(co)-darcos(co2))
 1500 if(qonly) goto 311
      do 16 j=1,4
      m=ia(j)
      if(m.le.0) go to 16
      m=m-1
      j1=3*(j-1)
      do 170 l=1,3
      m1=3*m+l
      l1=j1+l
      bmat(m1,i)=uv(l1)*c+bmat(m1,i)
 170  continue
   16 continue
      go to 311
 1310 error=6
      write(iout,1903) i
 1903 format(/,1x,41hatoms erronously defined,coordinate no.  ,i3,/,1x,
     110(4h****))
  310 nq=i
 1111 return
      end
      function darcos(x)
      implicit real*8(a-h,o-z)
      if(x.ge.1.0d0) goto 100
      if(x.le.-1.0d0) goto 200
      x1=dsqrt(1.0d0-x**2)
      if(dabs(x).lt.1.0d-11) goto 300
      s=datan(x1/x)
      if(x.lt.0.0d0) s=s+3.1415926536d0
      darcos=s
      return
  100 darcos=0.0d0
      return
  200 darcos=3.1415926536d0
      return
  300 darcos=1.5707963269d0
      return
      end
      subroutine vektor(u,r,i,j,xa,ya,za)
      implicit real*8(a-h,o-z)
      dimension u(3),xa(1),ya(1),za(1)
      u(1)=xa(i)-xa(j)
      u(2)=ya(i)-ya(j)
      u(3)=za(i)-za(j)
      r= dsqrt(scalar(u,u))
      call norm(u)
      return
      end
      subroutine osinv (a,n,d,tol,l,m)
      implicit real*8(a-h,o-z)
c
c     parameters:  a - input matrix , destroyed in computation and repla
c                      by resultant inverse (must be a general matrix)
c                  n - order of matrix a
c                  d - resultant determinant
c            l and m - work vectors of lenght n
c                tol - if pivot element is less than this parameter the
c                      matrix is taken for singular (usually = 1.0e-8)
c     a determinant of zero indicates that the matrix is singular
c
      dimension a(1), m(1), l(1)
      d=1.0d0
      nk=-n
      do 180 k=1,n
         nk=nk+n
         l(k)=k
         m(k)=k
         kk=nk+k
         biga=a(kk)
         do 20 j=k,n
            iz=n*(j-1)
         do 20 i=k,n
            ij=iz+i
c
c     10 follows
c
            if (dabs(biga)-dabs(a(ij))) 10,20,20
   10       biga=a(ij)
            l(k)=i
            m(k)=j
   20    continue
         j=l(k)
         if (j-k) 50,50,30
   30    ki=k-n
         do 40 i=1,n
            ki=ki+n
            holo=-a(ki)
            ji=ki-k+j
            a(ki)=a(ji)
   40    a(ji)=holo
   50    i=m(k)
         if (i-k) 80,80,60
   60    jp=n*(i-1)
         do 70 j=1,n
            jk=nk+j
            ji=jp+j
            holo=-a(jk)
            a(jk)=a(ji)
   70    a(ji)=holo
   80    if (dabs(biga)-tol) 90,100,100
   90    d=0.0d0
         return
  100    do 120 i=1,n
            if (i-k) 110,120,110
  110       ik=nk+i
            a(ik)=a(ik)/(-biga)
  120    continue
         do 150 i=1,n
            ik=nk+i
            ij=i-n
         do 150 j=1,n
            ij=ij+n
            if (i-k) 130,150,130
  130       if (j-k) 140,150,140
  140       kj=ij-i+k
            a(ij)=a(ik)*a(kj)+a(ij)
  150    continue
         kj=k-n
         do 170 j=1,n
            kj=kj+n
            if (j-k) 160,170,160
  160       a(kj)=a(kj)/biga
  170    continue
         d=d*biga
         a(kk)=1.0d0/biga
  180 continue
      k=n
  190 k=k-1
      if (k) 260,260,200
  200 i=l(k)
      if (i-k) 230,230,210
  210 jq=n*(k-1)
      jr=n*(i-1)
      do 220 j=1,n
         jk=jq+j
         holo=a(jk)
         ji=jr+j
         a(jk)=-a(ji)
  220 a(ji)=holo
  230 j=m(k)
      if (j-k) 190,190,240
  240 ki=k-n
      do 250 i=1,n
         ki=ki+n
         holo=a(ki)
         ji=ki+j-k
         a(ki)=-a(ji)
  250 a(ji)=holo
      go to 190
  260 return
      end
      subroutine relax(c,fi,qq,nq,qq1)
      implicit real*8(a-h,o-z)
      dimension c(1),fi(1),qq(1),qq1(1)
      common /files/ inp,iout,ix,ipun,nshout
      data ajoule/.22936758d0/,zero/0.0d0/,half/0.5d0/
      etot=0.d0
      ij=0
      edecr=zero
      do 480 i=1,nq
      s=zero
      do 470 j=1,nq
      ij=ij+1
      s=s+c(ij)*fi(j)
  470  continue
      qq1(i)=-s
      edecr=edecr-half*s*fi(i)
  480  continue
      edec=edecr*ajoule
      et=etot+edec
      write (iout,490) edecr,edec,et
  490 format(//,1x,26hprojected energy lowering=,f13.8,4x,7hajoule,4x,
     *f13.8,8h hartree,4x,15h total energy= ,f13.8,8h hartree)
      write (iout,1200)
      write (nshout,1200)
      do 500 i=1,nq
      qq(i)=qq(i)+qq1(i)
      write (iout,1160) i,qq1(i),qq(i)
      write (nshout,1160) i,qq1(i),qq(i)
  500 continue
 1160 format(2x,i2,4x,2f12.7)
 1200 format(/1x,44hgeometry change and new internal coordinates,/)
      return
      end
      subroutine conver(x,f,natoms,qq,fi,nq,xm,idsp,ncard)
      implicit real*8(a-h,o-z)
      logical fail
      dimension x(1),f(1),qq(1),fi(1),xm(1)
      common/work/bb(180,180),xa(60),ya(60),zza(60),ccc(16290),
     1cc(180),l(180),m(180)
      common /files/ inp,iout,ikk,ipun,nshout
      common /ite/ iter,intnl,fail,iix
      ix=iix
      do 10 i=1,natoms
      ii=3*i-2
      xa(i)=x(ii)
      ya(i)=x(ii+1)
  10  zza(i)=x(ii+2)
 1140 format(/,5x,16h determinant  = ,d16.5)
 1150 format(/1x,34hinternal coordinates and gradients,/)
 1160 format(2x,i2,4x,2f12.7)
 1170 format(8f10.7)
      rewind ix
      nek=3*natoms
      ncmax=180
      nqmax=180
      n=natoms
      nq=0
      call machb(nek,bb,ncmax,nqmax,xa,ya,zza,qq,n,ix,iout,nq,.false.,
     1 .false.,ncard)
      i1=0
      do 2240 i=1,nq
      do 2240 j=1,nq
      i1=i1+1
      s=0.d0
      do 2230 k=1,nek
 2230 s=s+bb(k,i)*bb(k,j)*xm(k)
 2240 ccc(i1)=s
      tol=1.d-8
      call osinv(ccc,nq,d,tol,l,m)
      if(d.lt.1.d-10) write (iout,1140)d
      if(idsp.ne.0)return
      write (iout,1150)
      write (nshout,1150)
      do 2260 i=1,nq
      ttt=0.0d0
      do 2250 j=1,nek
 2250 ttt=ttt+bb(j,i)*f(j)*xm(j)
 2260 cc(i)=ttt
      ij=0
      do 320 i=1,nq
      ttt=0.d0
      do 2280 j=1,nq
      ij=ij+1
 2280 ttt=ttt+ccc(ij)*cc(j)
      fi(i)=ttt
      write (iout,1160) i,qq(i),fi(i)
      write (nshout,1160) i,qq(i),fi(i)
      write (ipun,1170) qq(i),fi(i)
  320 continue
      return
      end
      subroutine displ(x1xx,qd,natoms,nq,xm,symb,ncard,ccin)
      implicit real*8(a-h,o-z)
      logical ltex,fail
      dimension x1xx(1),qd(1),xm(1),symb(1),ccin(180)
      common /info/ xcc(180),nat,ika,mult,ia(60),nnn
      common /work/ bb(180,180),xa(60),ya(60),zza(60),ccc(16290),
     1cc(180),f(180),fi(180),xy(180),qq1(180),qq(180),xin(180),
     2ndtt(12)
      common /files/ inp,iout,ikk,ipun,nshout
      common /ite/ iter,intnl,fail,iix
      integer symb
      ltex=.true.
      ix=iix
      idsp=1
      call conver(x1xx,f,natoms,qq,fi,nq,xm,idsp,ncard)
      do 200 i=1,nq
      xin(i)=i
  200 qq1(i)=ccin(i)
      nek=3*natoms
      ncmax=180
      nqmax=180
      n=natoms
      tol=1.d-8
      rewind ix
      write (iout,1250)
      write (iout,1260)  (qq(j),j=1,nq)
      write (iout,1270)  (xin(j),ccin(j),j=1,nq)
      write (nshout,1270)  (xin(j),ccin(j),j=1,nq)
      irep=0
      do 640 i=1,n
      xy(i)=xa(i)
       j=i+n
      xy(j)=ya(i)
      j=j+n
      xy(j)=zza(i)
  640  continue
  650  continue
      irep=irep+1
      ij=0
      do 670 i=1,nq
      s=0.0d0
      do 660 j=1,nq
      ij=ij+1
      s=s+ccc(ij)*qq1(j)
  660  continue
      fi(i)=s
  670  continue
      do 690 i=1,nek
      s=0.0d0
      do 680 j=1,nq
      s=s+bb(i,j)*fi(j)
  680  continue
      qq1(i)=s*xm(i)
  690  continue
      i3=0
      do 700 i=1,n
      i3=i3+3
      xy(i)=xy(i)+qq1(i3-2)
      j=i+n
      xy(j)=xy(j)+qq1(i3-1)
      j=j+n
      xy(j)=xy(j)+qq1(i3)
  700  continue
      i2=2*n+1
      do 710 i=1,ncard
      backspace ix
  710 continue
      call machb (nek,bb,ncmax,nqmax,xy(1),xy(n+1),xy(i2),qq1,n,ix,
     1iout,nq,.false.,.true.,ncard)
      write (iout,1280) irep
      write (iout,1260)  (qq1(i),i=1,nq)
      qmax=0.0d0
      do 720 i=1,nq
      qd(i)=qq1(i)
      qq1(i)=qq(i)+ccin(i)-qq1(i)
      if(qmax.lt.dabs(qq1(i))) qmax=dabs(qq1(i))
  720  continue
      if (irep.lt.15.and.qmax.gt.1.d-7) go to 650
      if (qmax.gt.1.d-6) write (iout,1330)
      if (qmax.gt.1.d-6) write (nshout,1330)
      write (iout,1290)
      write (ipun,1300)  (xin(i),ccin(i),i=1,4)
      i3=0
      do 740 i=1,n
      i3=i3+3
      ym=1.0d0/xm(i3)
      xiat=ia(i)
      j=i+n
      j1=j+n
       x=xy(i)-xa(i)
      y=xy(j)-ya(i)
      z=xy(j1)-zza(i)
      write (iout,1310)  i,xiat,x,y,z,xy(i),xy(j),xy(j1),ym
      xa(i)=xy(i)
      ya(i)=xy(j)
      zza(i)=xy(j1)
      if(ltex) write (ipun,1320) symb(i),xiat,xy(i),xy(j),xy(j1),ym
  740 continue
      do 760 i=1,n
      j=i+n
      j1=j+n
      i1=3*i-2
      i2=i1+1
      i3=i2+1
      x1xx(i1)=xy(i)
      x1xx(i2)=xy(j)
      x1xx(i3)=xy(j1)
  760 continue
 1250 format (/1x,30horiginal internal coordinates ,/)
 1260 format (2x,5f12.7,3x,5f12.7)
 1270 format (/,1x,21hinternal displacement,4(f4.0,f12.6))
 1280 format (1x,5hstep=,i4)
 1290 format (/,1x,53hcartesian displacements and new cartesian coordina
     1tes,/)
 1300 format (7hdispl. ,4(f3.0,1h,,f6.4),3x,3(a4,a6))
 1310 format (1x,i2,2h n,f3.0,3f10.6,1x,3f11.7,1x,f6.3)
 1320 format (2hn=,a4,4x,5f10.6)
 1330 format (/,1x,48h*****caution,within 15 steps no convergence*****/)
      return
      end
      subroutine setup(ih,lopt,iop,nfix,ifix,kell,tol,natp,m,nqlim,inp
     x ,iout,xita,wk,iix)
      implicit real*8 (a-h,o-z)
      dimension wk(iop),lopt(iop),ifix(nqlim)
      character*4 wk,wo
      character*8 a(10),astp
      parameter(one=1.0d0,ten=10.d0)
      data astp/'stop'/
c
c      set some deault parameter
c
      m=1
      tol=1.d-4
      ih=0
      nfix=0
c
c      read options
c
  10  read(inp,500) wo,x
      do 20 i=1,iop
      if(wk(i).eq.wo) goto 30
  20  continue
      write(iout,511) wo
  30  continue
      if(i.eq.7) then
      backspace inp
  50  read(inp,502) a
      write(iix,502) a
      if(a(1).eq.astp) goto 40
      goto 50
      endif
      lopt(i)=1
      write(iout,510) wo,x
      if(i.eq.3) kell=x+0.01d0
      if(i.eq.4) then
      iex=x+0.01d0
      tol=ten**(-iex)
      endif
      if(i.eq.5) ih=2
      if(i.eq.8) then
      nfix=nfix+1
      ifix(nfix)=x+0.01d0
      endif
      if(i.eq.9) then
      ih=x+0.01d0
      if(ih.ge.2)ih=2
      if(ih.eq.1)ih=1
      endif
      if(i.eq.10) natp=x+0.01d0
      if(i.eq.6)xita=one
      goto 10
 40   continue
  500 format (a4,6x,f10.5)
  510 format(1x,a4,2x,13hoption is on ,2f10.5)
  511 format(1x,'not such option',2x,a4)
  502 format(10a)
      return
      end
      subroutine objet(xc,en,g,natoms,inp2,iout,natp,nmax,inp1,nshout,
     x niter,symb,xm)
      implicit real*8 (a-h,o-z)
      common/vari/type(50),numb1(50),numb2(50),chg(50),iatno1(50)
      integer type,numb1,numb2,symb
      dimension xc(nmax,3),g(*),symb(*),xm(*),ia(50)
      data  one/1.0d0/
      rewind inp2
      read(inp2,900)niter
      write(iout,910)niter
      write (iout,530)
      na=0
      do 70 i=1,natp
      read(inp2,920)type(i),numb1(i),numb2(i),chg(i),iatno1(i)
      numb=numb2(i)
      do 71 j=1,numb
      na=na+1
      i3=na*3
      symb(na)=type(i)
      ia(na)=chg(i)+0.1
      read(inp2,930)xc(na,1),xc(na,2),xc(na,3),xm(i3)
         if (xm(i3).eq.0.0) xm(i3)=1.0
      write (iout,550) na,ia(na),xc(na,1),xc(na,2),xc(na,3),xm(i3)
         xm(i3)=one/xm(i3)
         xm(i3-1)=xm(i3)
         xm(i3-2)=xm(i3)
   71 continue
   70 continue
      natoms=na
      if(na.gt.0.and.na.le.nmax)go to 72
      write(iout,520)na
      stop 'objet'
   72 continue
      nek=3*na
      read (inp1,560) (g(i),i=1,nek)
c     do 73 i=1,nek
c     g(i)=-g(i)
c  73 continue
      write (iout,570)
      write (iout,580) (g(i),i=1,nek)
      write (nshout,570)
      write (nshout,580) (g(i),i=1,nek)
      xfqx=g(1)
      do 80 i=2,nek
   80 xfqx=xfqx+g(i)
      if (dabs(xfqx).gt.1.0d-5) write (iout,590) xfqx
      if (dabs(xfqx).gt.1.0d-5) write (nshout,590) xfqx
      if (dabs(xfqx).gt.1.0d-5) stop 'objet'
  520 format (1x,16htoo many nuclei ,i5)
  530 format (20x,'nuclear coordinates in a.u.'/)
  550 format (5x,i2,3x,i2,3f12.7,3x,f12.6)
  560 format (3d15.7)
  570 format (//,1x,'gradients'/)
  580 format (1x,3f12.6)
  590 format (/,1x,26hforces do not vanish, sum=,f15.7,/)
  900 format(i4)
  910 format('1',/t21,'program "gdiis 1.0"',
     +  //t10,'algorithm to optimizing molecular geometry',
     +  /t10,'the version based on subroutin in program "geomo"',
     +  /t10,'from p. csaszar and a.g. csaszar, elte budapest',
     + //t10,'reference:',t25,'p. csaszar and p. pulay,',
     +  ' j. mol. struct. 114, 31 (1984).',
     +  //t10,'programmed by: peter g. szalay',
     +  /t25,'institute for theoretical chemistry',
     +  /t25,'of univesity vienna',
     +  /t25,'waehringer str 17,  a-1090 wien ',
     +  /t25,'austria',
     +  /t25,'phone (222)43-61-41/70',
     +  /t25,'bitnet: a8441gad at awiuni11'
     +  //t10,'version date:',t25,'06-aug-87',
     +  //////,'          iteration number',i4,///)
  920 format(a3,2i3,f3.0,i3)
  930 format(4f14.8)
      return
      end
      subroutine dump(m,kell,nq,natoms,niter,xita,x,natp,lopt,gdiisf,
     x hh,inp2,iout,xm,nshout)
      implicit real*8 (a-h,o-z)
      common/work1/xx(180,60),ee(180,60),gg(180,60),qq(180,60)
      common/vari/type(50),numb1(50),numb2(50),chg(50),iatno1(50)
      integer type,numb1,numb2
      parameter(iop=10)
      dimension x(*),hh(*) ,lopt(iop),xm(*)
      integer gdiisf
      rewind gdiisf
      if((lopt(5).gt.0).or.(kell.le.1)) goto 201
      write(gdiisf) m
      lang=kell-1
      if(kell.ge.m) lang=m-1
      write(gdiisf) kell,lang,nq,natp,lopt
c&&
c     print*,' m',m,' kell',kell,' lang',lang
c&&
      write(gdiisf)((ee(i,j),j=1,lang),i=1,nq)
      write(gdiisf)((qq(i,j),j=1,lang),i=1,nq)
      write(gdiisf)((gg(i,j),j=1,lang),i=1,nq)
      nqq=nq**2
      write(gdiisf)(hh(i),i=1,nqq)
      write(gdiisf) xita
 201  rewind inp2
      write(inp2,900) niter
      write(iout,772)
      write(nshout,772)
      jj=0
      do 491 i=1,natp
      write(inp2,920)type(i),numb1(i),numb2(i),chg(i),iatno1(i)
      numb=numb2(i)
      xiat=chg(i)
      do 492 j=1,numb
      jj=jj+1
      jj1=(jj-1)*3+1
      jj2=jj1+1
      jj3=jj2+1
      xa=x(jj1)
      ya=x(jj2)
      za=x(jj3)
      xxmk=1.d0/xm(3*jj)
      write(inp2,940)xa,ya,za,xxmk
      write(iout,792)jj,xiat,xa,ya,za,xxmk
      write(nshout,792)jj,xiat,xa,ya,za,xxmk
  492 continue
  491 continue
      if((lopt(5).gt.0).or.(kell.le.1))then
      jj=1
      do 485 i=1,nq
      ij1=jj+i-1
      write(inp2,720)(hh(j),j=jj,ij1)
      jj=jj+nq
  485 continue
  101 format(10i8)
  102 format(10a8)
  720 format (8f10.5)
  772 format(/' cartesian coordinates in a.u.'/)
  792 format(1x,i3,4h   n,f3.0,4f12.7)
  900 format(i4)
  920 format(a3,2i3,f3.0,i3)
  940 format(4f14.5,24x)
      endif
      return
      end
      subroutine rest(m,kell,lang,nq,xita,lopt,gdiisf,hh,natp)
      implicit real*8 (a-h,o-z)
      common/work1/xx(180,60),ee(180,60),gg(180,60),qq(180,60)
c     common/vari/type(50),numb1(50),numb2(50),chg(50),iatno1(50)
c     integer type,numb1,numb2
      parameter(iop=10)
      dimension hh(*) ,lopt(iop)
      integer gdiisf
      rewind gdiisf
      read(gdiisf) m
      if(m.le.1) stop ' m in rest kleiner als 2'
      read(gdiisf) kell,lang,nq,natp,lopt
      print*,' m=',m,'kell=',kell,'lang=',lang
      read(gdiisf)((ee(i,j),j=1,lang),i=1,nq)
      read(gdiisf)((qq(i,j),j=1,lang),i=1,nq)
      read(gdiisf)((gg(i,j),j=1,lang),i=1,nq)
      nqq=nq**2
      read(gdiisf)(hh(i),i=1,nqq)
      read(gdiisf) xita
      lang=lang+1
  101 format(10i8)
  102 format(10a8)
      return
      end
      subroutine dzero(n,a,ia)
c
c  zero the elements of the vector a().
c   n  = number of elements in the vector.
c   a  = working precision vector.
c   ia = spacing between vector elements.
c
c  written by ron shepard.
c
      double precision a(ia,*)
      double precision zero
      parameter(zero=0d0)
c
      if(n.le.0)return
      do 10 i=1,n
10    a(1,i)=zero
      return
      end
      real*8 function dot()
      dot = 0.d0
      write (6,'('' ****** error''/'' ****** dot called'')')
      stop 'gdiis called dot'
      end
      real*8 function updot()
      updot = 0.d0
      write (6,'('' ****** error''/'' ****** updot called'')')
      stop 'gdiis called updot'
      end