*====================================================================
* INPUT ROUTINES: ICOCOM  IPOCOM     IANCOM     ISACOM ICFCOM IZMCOM
*====================================================================
      subroutine icocom(str,ifc,nx,ifx,a,lfn,mxn,noa,nmz,imz)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifh,ifc
      character str*(*),prmp*79,ct*1
      dimension fnum(10),a(nx,3),imz(nx)
 100  it=1
      ct='x'
      prmp=
     +' Enter:(n) x  y  z      or: ct      (h = help)'
      prmp(75:79)='*ico*'
      write(6,*) prmp
 98   str='ico>'
      if(it.eq.2) str='icoC'
      if(it.eq.3) str='icoS'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         if(it.le.1) then
            write(6,1000)
            go to 100
         elseif(it.eq.2) then
            write(6,1002)
            go to 90
         elseif(it.eq.3) then
            write(6,1003)
            go to 90
         endif
      endif
      call chtflp(str(1:last),nun,fnum)
      if(nun.ne.4.and.nun.ne.3.and.nun.ne.0) then
 80      write(6,*) 
     +' Input wrong, input line must normally contain 3 or 4 numbers.'
         if(it.gt.1) go to 90
         go to 100
      endif
      if(nun.eq.0) then
         if(index(str,'s').ne.0) then
            it=3
            ct='s'
            go to 90
         elseif(index(str,'c').ne.0) then
            it=2
            ct='c'
            go to 90
         elseif(index(str,'x').ne.0) then
            it=1
            ct='x'
            go to 90
         else
            return
         endif
      endif
      if(nun.eq.4) n=int(fnum(1))
      if(nun.eq.3) n=lfn
      iix=0
      if(nun.eq.3) iix=1
      if(fls(nx,n)) go to 200
      if(nun.ge.3) then
         ca=fnum(2-iix)
         cb=fnum(3-iix)
         cc=fnum(4-iix)
         go to 97
      endif
      it=n
      if(it.gt.3.or.it.lt.1) return
      go to 90
   97 call nchek(str,ifc,nx,ifx,n,io)
      if(ifc) return
      if(io.eq.0) go to 99
   90 if(it.eq.1) prmp=' Enter:(n) x  y  z                (h = help)'
      if(it.eq.2) prmp=' Enter:(n) rh  phi  z             (h = help)'
      if(it.eq.3) prmp=' Enter:(n) r  tht  phi            (h = help)'
      prmp(74:79)=ct//'*ico*'
      write(6,*) prmp
      go to 98
   99 if(it.eq.2.or.it.eq.3) go to 20
   10 a(n,1)=ca
      a(n,2)=cb
      a(n,3)=cc
      ifx(n)=.true.
      write(6,5) n,ca,cb,cc
    5 format(2x,i3,3x,3(2x,f15.9))
      call naccnt(nx,ifx,lfn,mxn,noa)
      go to 90
   20 if(it.eq.3) go to 30
      call cyli(x,y,ca,cb,1)
      ca=x
      cb=y
      go to 10
   30 call spheri(x,y,z,ca,cc,cb,1)
      ca=x
      cb=y
      cc=z
      go to 10
 200  write(6,1010) nx
 1010 format(' Atom number n must be in range  1  to ',i3)
      go to 90
 1000          format(
     +'     Enter n, the number of the atom (if you wish)'/
     +' and CARTESIAN coordinates x y z  (Example of input'/
     +' line: 5  1.3408  -.31288  2.00764 ). Note that'/
     +' exponential notation cannot be used for input.'/
     +'     If n is not given the smallest free number is'/
     +' automatically choosen. If you give a number n'/
     +' which already has been used you are warned and'/
     +' must confirm it to replace old input.'/
     +'     To change to prompt for other coordinate type (ct)'/
     +' (CYLINDRICAL or SPHERICAL) Type the corresponding'/
     +' coordinate type (ct= c or s respectively).        '/
     +'     The number n and the stored cartesian coordi-'/
     +' nates (even if you are not using them) are echoed to '/
     +' the terminal after each input line.              '/
     +'     Back for new command, when finished, as usual'/
     +' with an empty line.')
 1002 format(
     +'   CYLINDRICAL COORDINATES are defined as follows:'/
     +' rh =  distance to z-axis (rh stands for the Greek letter rho).'/
     +' phi = angle between positive x-axis and the projection of the'/
     +'       point''s position vector in the xy-plane. phi is between'/
     +'       -180 and 180 degrees and has the same sign as the y-'/
     +'       coordinate of the projection.'/
     +' z =   the ordinary z-coordinate of the point.')
 1003 format(
     +'   SPHERICAL COORDINATES are defined as follows:'/
     +' r =   distance from point to origin'/
     +' tht = angle between positive z-axis and the point''s position'/
     +'       vector (between 0 and 180 degrees). (tht stands for the'/
     +'       Greek letter theta.)'/
     +' phi = angle between positive x-axis and the projection of the'/
     +'       point''s position vector in the xy-plane. phi is between'/
     +'       -180 and 180 degrees and has the same sign as the y-'/
     +'       coordinate of the projection.')
      end
*=====================================================================
      subroutine ipocom(str,ifc,nx,ifx,nq,a,mxn,ifz)
      implicit double precision(a-h,o-z)
      logical ifx(nx),ifz,ifzloc,ifh,ifc
      dimension nq(nx),a(nx,3),fnum(10),v1(3)
      character str*(*)
 97   write(6,*)
     +' Enter: name  (  size (  Z))                (h = help)     ',
     +'               *ipo*'
      str='ipo>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         if(last.eq.1) then
            write(6,1000)
            go to 97
         else
            write(6,1001)
            go to 97
         endif
      endif
      call chtflp(str(1:last),nuno,fnum)
      if(nuno.gt.3) go to 99
      nu=int(fnum(1))
      ndim=2
      if(index(str,'t ').ne.0) then
         nu=4
         ndim=3
      elseif(index(str,'o ').ne.0) then
         nu=6
         ndim=3
      elseif(index(str,'c ').ne.0) then
         nu=8
         ndim=3
      elseif(index(str,'i ').ne.0) then
         nu=12
         ndim=3
      elseif(index(str,'d ').ne.0) then
         nu=20
         ndim=3
      endif
      if(ndim.eq.2.and.nuno.eq.0) go to 99
      if(nu.lt.2) go to 99
      if(mxn+nu.gt.nx) then
         write(6,*) ' Not room enough for polygon or -hedron.'
         return
      endif
      r=1.D0
      it=0
      if(nuno.ge.1) then
         if(ndim.eq.2.and.nuno.ge.2) then
            r=fnum(2)
            it=2
         elseif(ndim.eq.3) then
            r=fnum(1)
            it=2
         endif
      endif
      if(r.lt.0.D0) then
         r=-r
         it=1
      elseif(r.eq.0.D0) then
         r=1.D0
         it=0
      endif
      nqz=1
      ifzloc=.false.
      if((ndim.eq.3.and.nuno.ge.2).or.(ndim.eq.2.and.nuno.ge.3)) then
         if(ndim.eq.2) nqz=int(fnum(3))
         if(ndim.eq.3) nqz=int(fnum(2))
         nqz=min(nqz,99)
         nqz=max(nqz,-9)
         ifz=.true.
         ifzloc=.true.
      endif
      mf=mxn+1
      ml=mxn+nu
      if(ndim.eq.3) go to 90
      pi=dacos(-1.D0)
      if(it.eq.2) r=r/(2.D0*dsin(pi/float(nu)))
      do 10 kk=mf,ml
         k=kk-mf
         arg=2.D0*float(k)*pi/float(nu)
         a(kk,1)=r*dcos(arg)
         a(kk,2)=r*dsin(arg)
         a(kk,3)=0.D0
         if(ifzloc) nq(kk)=nqz
 10      ifx(kk)=.true.
      go to 200
 90   continue
      if(nu.ne.4.and.nu.ne.6.and.nu.ne.8.and.nu.ne.12.and.nu.ne.20)
     +     go to 99
      if(nu.eq.6.or.nu.eq.12) go to 91
      mpl=min0(ml,mf+7)
      do 20 kk=mf,mpl
         ifx(kk)=.true.
         if(ifzloc) nq(kk)=nqz
         do 20 j=1,3
 20         a(kk,j)=1.D0
      neg=-1.D0
      a(mf+1,1)=neg
      a(mf+1,2)=neg
      a(mf+2,1)=neg
      a(mf+2,3)=neg
      a(mf+3,2)=neg
      a(mf+3,3)=neg
      if(nu.eq.4) go to 96
      do 22 i=1,3
 22      a(mf+3+i,i)=neg
      do 21 i=1,3
 21      a(mf+7,i)=neg
      if(nu.eq.8) go to 96
 91   mp=mf
      if(nu.eq.20) mp=mf+8
      do 30 kk=mp,ml
         ifx(kk)=.true.
         if(ifzloc) nq(kk)=nqz
         do 30 i=1,3
 30         a(kk,i)=0.D0
      tau=(1.D0+dsqrt(5.D0))/2.D0
      ut=1.D0
      if(nu.eq.20) ut=1.D0/tau
      do 31 ii=1,2
         i=(ii-1)*6
         a(mp+i,3)=ut
         a(mp+i+1,1)=ut
         a(mp+i+2,2)=ut
         do 33 j=1,3
 33         a(mp+i+2+j,j)=-ut
         if(nu.eq.6) go to 96
 31      continue
      do 32 ii=1,2
         i=(ii-1)*6
         a(mp+i,2)=tau
         a(mp+i+1,3)=tau
         a(mp+i+2,1)=tau
         a(mp+i+3,3)=tau
         a(mp+i+4,1)=tau
         a(mp+i+5,2)=tau
 32      tau=-tau
 96   continue
      edg=dsqrt(2.D0)
      if(nu.eq.4) edg=2.D0*edg
      if(nu.eq.8.or.nu.eq.12) edg=2.D0
      if(nu.eq.20) edg=2.D0*ut      
      fac=1.D0/edg
C
C Handle truncated polyhedra below      
C
      if(ndim.eq.3.and.index(str,' ').eq.3) then
         do 40 kk=mf,ml
            do 41 j=1,3
               a(kk,j)=fac*a(kk,j)
 41         continue
 40      continue
         if(str(1:1).eq.'t') then
            if(nu.eq.4.or.nu.eq.6.or.nu.eq.12) then
               xa=1.D0/3.D0
            elseif(nu.eq.8) then
               xa=1.D0/(2.D0+dsqrt(2.D0))
            elseif(nu.eq.20) then
               xa=2.D0/(5.D0+dsqrt(5.D0))
            endif
         elseif(str(1:1).eq.'d') then
            xa=0.5D0
         else
            go to 99
         endif
         nnew=0
         do 50 kk=mf,ml
            do 51 kl=kk+1,ml
               call bij(nx,a,kk,kl,v1,dd,ierr)
               if(dabs(1.D0-dd).gt.1.D-9) go to 51
               nnew=nnew+1
               do 52 j=1,3
                  mn=ml+nnew
                  if(mn.gt.nx) go to 98
                  a(mn,j)=a(kk,j)+xa*v1(j)
 52            continue
               if(xa.eq.0.5D0) go to 51
               nnew=nnew+1
               do 53 j=1,3
                  mn=ml+nnew
                  if(mn.gt.nx) go to 98
                  a(mn,j)=a(kl,j)-xa*v1(j)
 53            continue               
 51         continue
 50      continue         
      nu=nnew
      dmin=100.
      do 60 kk=mf,mf+nu-1
         do 61 j=1,3
            a(kk,j)=a(ml+(kk-mf)+1,j)
 61      continue
         if(kk.lt.mf+nu-1) then
            call bij(nx,a,mn,kk,v1,dd,ierr)
            if(dd.lt.dmin) dmin=dd
         endif
         ifx(kk)=.true.
         if(ifzloc) nq(kk)=nqz
 60   continue
      ml=mf+nu-1
      edg=1.D0-2*xa
      if(xa.eq.0.5D0) edg=dmin
      endif
C
C Truncation section of routine ends here
C
      if(it.ne.1.and.it.ne.2) go to 200
      dd=dsqrt(a(mf,1)**2+a(mf,2)**2+a(mf,3)**2)
      if(it.eq.1) fac=r/dd
      if(it.eq.2) fac=r/edg
      do 70 kk=mf,ml
         do 71 j=1,3
            a(kk,j) = fac*a(kk,j)
 71      continue
 70   continue
 200  call order(nx,ifx,nq,a,mf,ml,3,.false.)
      write(6,100) ndim,nu,mf,ml
 100  format(2x,i2,'-dimensional',i4,'-gon in range:',i3,' -',i3)
      return
 98   write(6,*) ' Not enough workspace to truncate.'
      return
 99   write(6,*) ' Input to polygon & polyhedron wrong.'
      go to 97
 1000    format(/
     +'   To generate a regular POLYGON (equilateral triangle, '/
     +' square, pentagon, or hexagon etc)'/
     +'   Enter: name = number of vertices wanted'/
     +' (i.e. name = 3, 4, 5, or 6 etc).'//
     +'   To generate one of the five regular (platonic) POLYHEDRA: '/
     +' Tetrahedron, Octahedron, Cube, Icosahedron, or Dodecahedron,'/
     +'   Enter: name = t ,  o ,  c ,  i , or  d  respectively'/
     +'   These  have   4,   6,   8,  12, and 20 vertices, resp.'//
     +'   If you wish to specify the length of an EDGE,'/
     +' put:   size= length of edge.'/
     +'   If you wish to specify the RADIUS, '/
     +' put:   size= -radius. '//
     +'   If size has been given you can also give the atomic number Z'/
     +' (charge) last in the input line.'/
     +'   (If size is not given a standard size is choosen; if atomic'/
     +' number is not given it will remain marked off or as it was.)'/
     +' The points will be centered on the origin and placed after the'/
     +' highest numbered entered atom.'/
     +'                                   (ht= help about truncation)')
 1001    format(/'    TRUNCATED REGULAR POLYHEDRA'/
     +'   To generate points at the vertices of one of the five semi-'/
     +' regular Archimedean polyhedra corresponding to the Truncated'/
     +' Tetrahedron, Octahedron, Cube, Icosahedron, or Dodecahedron,'/
     +'   Enter: name = tt ,  to ,  tc ,  ti , or  td  respectively'//
     +' These have 12, 24, 24, 60, and 60 vertices respectively.'/
     +' Note: the Truncated Icosahedron is the C60 football.'//
     +'   To generate the two new, quasi-regular, polyhedra that can'/
     +' be obtained by truncating the five Platonic polyhedra in such'/
     +' a way that each edge is replaced by a point (deep truncation),'/
     +'   Enter: name = dt ,  do ,  dc ,  di , or  dd  respectively'//
     +' Here dt          gives the  Octahedron again, but:'/
     +'      do and dc   give the   Cuboctahedron, and'/
     +'      di and dd   give the   Icosidodecahedron.'/
     +' These last two are the, so called, quasi-regular polyhedra.'/
     +' They have 12 and 30 vertices respectively.'/)
      end
*====================================================================
      subroutine iancom(str,ifc,nx,ifx,nq,noa,ifz)
      logical ifx(nx),fls,ifz,ifh,ifc
      character str*(*), atmn*2, atmnme*2,prmp*79
      integer nq(nx),nx,noa,ivec(10),nqq
 99   prmp=
     +' Enter: Z   or: n1 Z1 ( n2 Z2..( n5 Z5))   or: nf nl Z'//
     +'     (h = help)     *ian*'
      write(6,*) prmp
      str='ian>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh.and.last.eq.1) then
         write(6,1000)
         write(6,1001)
         go to 99
      endif
      call intstr(str(1:last),nun,ivec)
 500  continue
      if(nun.eq.0) then
         if(last.eq.2.and.str(1:2).eq.'rm') then
            do 1 i=1,nx
 1             nq(i)=1
            ifz=.false.
            write(6,*) ' Atomic numbers Z removed.'
            return
         else
            if(last.gt.2) go to 91
            atmn=str(1:last)
            nqz=nqq(atmn)
            if(nqz.gt.99) then
               write(6,*) ' Chemical symbol not recognized.'
               go to 99
            endif
            write(6,1002) atmn,nqz
 1002       format('   Atom: ',a2,'  has atomic number Z= ',i2)
            go to 99
         endif
      elseif(nun.eq.1) then
         nin=1
         ncount=1
         nst=1
 98      continue
         nct=1
         do 4 i=nst,nx
            if(noa.ne.0.and..not.ifx(i)) go to 4
            nqz=ivec(nct)
            if(nqz.gt.99.or.nqz.lt.-9) then
               write(6,*) ' Atomic number Z must be in range -9 to 99'
               nqz=min(nqz,99)
               nqz=max(nqz,-9)
            endif
            nq(i)=nqz
            write(6,1003) i,nqz,atmnme(nqz)
            ifz=.true.
            nct=nct+1
            ncount=ncount+1
            if(noa.ne.0.and.ncount.gt.noa) go to 93
            if(nct.gt.nin) go to 5
 4       continue
         if(i.gt.nx) go to 94
 5       nst=i+1
 97      prmp=
     +   ' Enter: Z1 ( Z2 ( Z3..( Z10)..))              (h = help)'
         prmp(74:79)='Z*ian*'
         write(6,*) prmp
         str='ianZ'
         call readin(ifc,str,last,ifh)
         if(ifc) return
         if(last.eq.0) go to 99
         if(ifh.and.last.eq.1) then
            write(6,1005)
            go to 97
         endif
         call intstr(str(1:last),nin,ivec)
         if(nin.eq.0) then
            if(last.gt.2) then
               write(6,*) ' Input to Z-only input wrong.'
               go to 97
            endif
            atmn=str(1:last)
            nqz=nqq(atmn)
            if(nqz.gt.99) then
               write(6,*) ' Chemical symbol not recognized.'
               go to 97
            endif
            write(6,1002) atmn,nqz
            go to 97
         endif
         go to 98
      elseif(nun.eq.3) then
         nf=ivec(1)
         nl=ivec(2)
         nqz=ivec(3)
         if(nqz.gt.99.or.nqz.lt.-9) go to 92
         if(fls(nx,nf).or.fls(nx,nl)) go to 94
         if(nf.gt.nl) go to 91
         do 10 n=nf,nl
            nq(n)=nqz
 10      continue
         write(6,1004) nf,nl,nqz,atmnme(nqz)
         ifz=.true.
         go to 99
      elseif(nun.eq.2.or.nun.gt.3) then
         if(mod(nun,2).ne.0) go to 91
         do 20 ii=1,nun/2
            n=ivec(2*ii-1)
            if(fls(nx,n)) go to 94
            nqz=ivec(2*ii)
            if(nqz.gt.99.or.nqz.lt.-9) go to 92
            nq(n)=nqz
            write(6,1003) n,nqz,atmnme(nqz)
 20      continue
         ifz=.true.
         go to 99
      endif
 1003 format(
     +' Atom no:',i3,'   given  Z=',i2,'  i.e. is a  ',a2,'-atom.')
 1004 format(
     +' Range ',i3,' -',i3,'  given  Z=',i2,'  i.e. are  ',a2,'-atoms.')
 1005 format(
     +'   If you make a MISTAKE during Z-only-input and have many'/
     +' charges to enter, don''t stop! (This is done by returning an'/
     +' empty line.) Note the number n and continue your input. Since'/
     +' the Z-only-input mode always restarts at the first atom you'/
     +' must otherwise do a lot of work again. When you are finished'/
     +' you can correct your mistakes using the  n, Z - input mode.'/
     +'   Now continue with your next Z-value.')
 91   write(6,*) ' Atomic number input wrong.'
      go to 99
 92   write(6,*) ' Atomic number Z must be in range -9 to 99'
      go to 99
 93   write(6,*) ' All entered atoms have now been given charge.'
      return
 94   write(6,1010) nx
 1010 format(' Atom number n must be in range  1  to ',i3)
      go to 99
 1000       format(
     +'   There are three different modes of input of atomic number Z:'/
     +'          (1)  INPUT OF Z only'/
     +'   If the first input line contains only one number the corre-'/
     +' sponding Z-value is given to the first atom (if you have not '/
     +' entered any atoms to atom 1, else to the lowest numbered en-'/
     +' tered atom). After this you will be prompted for Z-values only'/
     +' and they will be assigned to (entered) atoms in the order that'/
     +' they are given (now from 2 or second lowest and so on).'/
     +'          (2) INPUT OF n and Z'/
     +'   If you give an even number of numbers in the first line the'/
     +' first number in each pair is interpreted as the atom number n'/
     +' and the second as its Z-value. You can give up to 5 pairs (10'/
     +' numbers) per line.')
 1001       format(
     +'          (3) INPUT OF nf  nl  Z   (range and atomic number)'/
     +'   Three numbers in the input line will be taken to be  nf, nl'/
     +' and Z, i.e. the atoms in the corresponding range get charge Z.'/
     +'          SPECIAL OPTIONS'/
     +'   Type:rm    to remove all charges (all Z=1 and charge marked '/
     +' as off). Atomic number will then not be displayed.'/
     +'   Type:X =chem. symb. of atom  to get its Z displayed (e.g. Fe'/
     +' will display 26, note use of capital and lower case letters).')
      end
*====================================================================
      subroutine isacom(str,ifc,nx,ifx,nq,a,noa,mxn,ifz,nnh,nmz,imz,mz)
      double precision a(nx,3),x,y,z
      logical ifx(nx),ifz,ifzr,fex,open,ifc,ifh,empty
      character prmp*79,str*(*),numbst*3,line*10
      integer nx,nq(nx),nnh,noa,mxn, nmz,imz(nx), mz(nx,3),last
      integer k(15),in(5),jn(5),kn(5),izn(15), jmx, noar,nmzr,nnhr
      integer i,j,io,ncount,n,nsave,nsaves,ns,nun,ivec(10)
      inquire(unit=3, opened=open)
      if (open) go to 100
      inquire(file='ecldsv.dat',exist=fex)
      if(.not.fex) then
         write(6,*) ' There is no saved file (ecldsv.dat).'
         return
      endif
      open (3,file='ecldsv.dat', err=91, status= 'old')
 100  rewind 3
*
*  Read the file till end of file and check it:
      do 106 nsaves=1,900
*
         read(3,30,iostat=io) noar,ifzr,nnhr, nmzr
         if(io.ne.0) go to 90
*
         if(noar.gt.0) then
            do 102 i=1,noar
               read(3,40,iostat=io) n, x, y, z
               if(io.ne.0) go to 92
 102        continue
         endif
         if(nnhr.gt.0) then
            do 103 i=1,nnhr/15
               read(3,50,iostat=io) (k(j),izn(j), j=1,15)
               if(io.ne.0) go to 92
 103        continue
            jmx=nnhr-15*(nnhr/15)
            read(3,50,iostat=io) (k(j),izn(j), j=1,jmx)
            if(io.ne.0) go to 92
         endif
         if(nmzr.gt.0) then
            do 104 i=1,nmzr/5
               read(3,60,iostat=io)(k(j),in(j),jn(j),kn(j),izn(j),j=1,5)
               if(io.ne.0) go to 92
 104        continue
            jmx=nmzr-5*(nmzr/5)
            read(3,60,iostat=io)(k(j),in(j),jn(j),kn(j),izn(j),j=1,jmx)
            if(io.ne.0) go to 92
         endif
         read(3,70,iostat=io) line
         if(io.ne.0) go to 90
*     
 106  continue
*
      write(6,*) ' Too many saves on file (more than 900)!'
      return
*
 90   continue
      if(io.gt.0) then
         write(6,*) 
     +        ' Check of save file ended with IO status',io,' after',
     +        nsaves,' saves (or one less).'
         return
      endif
      nsave=nsaves-1
      empty=(noa.eq.0.and.(.not.ifz).and.nmz.eq.0) 
*
      if(.not.empty) then
         prmp=
     +' Enter: number of save  (last save is '//numbst(nsave)//')'//
     +'    (h = help)'
         prmp(75:79)='*isa*'
 80      write(6,*) prmp
         str='isa>'
         call readin(ifc,str,last,ifh)
         if(last.eq.0.or.ifc) return
         if(ifh) then
            write(6,2500)
            go to 80
         endif
         call intstr(str(1:last),nun,ivec)
         if(nun.le.0.or.nun.ge.2) then
            write(6,*) ' Input to  isa  wrong.'
            go to 80
         endif
         nsave=ivec(1)
         if(nsave.lt.1.or.nsave.ge.nsaves) then
            write(6,*) ' Input to  isa  wrong.'
            go to 80
         endif
      endif
*
      rewind 3
      do 6 ns=1,nsave
*
         ncount=0
         read(3,30,iostat=io) noar,ifzr,nnhr, nmzr
         if(io.ne.0) go to 92
         if(ns.eq.nsave) then
            noa=noar
            ifz=ifzr
            nnh=nnhr
            nmz=nmzr
            do 1 i=1,nx
               ifx(i)=.false.
               nq(i)=1
               imz(i)=0
 1          continue
         endif
*
         do 2 i=1,noar
            read(3,40,iostat=io) n, x, y, z
            if(io.ne.0) go to 92
            if(ns.lt.nsave) go to 2
            ncount=ncount+1
            a(n,1)=x
            a(n,2)=y
            a(n,3)=z
            ifx(n)=.true.
 2       continue
         if(nnhr.gt.0) then
            do 3 i=1,nnhr/15
               read(3,50,iostat=io) (k(j),izn(j), j=1,15)
               if(io.ne.0) go to 92
               if(ns.lt.nsave) go to 3
               do 23 j=1,15
                  nq(k(j)) = izn(j)
 23            continue
 3          continue
            jmx=nnhr-15*(nnhr/15)
            read(3,50,iostat=io) (k(j),izn(j), j=1,jmx)
            if(io.ne.0) go to 92
            if(ns.lt.nsave) go to 25
            do 24 j=1,jmx
               nq(k(j)) = izn(j)
 24         continue
 25         continue
         endif
         if(nmzr.eq.0) go to 5
         do 4 i=1,nmzr/5
            read(3,60,iostat=io) (k(j),in(j),jn(j),kn(j),izn(j),j=1,5)
            if(io.ne.0) go to 92
            if(ns.lt.nsave) go to 4
            do 33 j=1,5
               mz(k(j),1)=in(j)
               mz(k(j),2)=jn(j)
               mz(k(j),3)=kn(j)
               imz(k(j))=izn(j)
 33         continue
 4       continue
         jmx=nmzr-5*(nmzr/5)
         read(3,60,iostat=io)(k(j),in(j),jn(j),kn(j),izn(j),j=1,jmx)
         if(io.ne.0) go to 92
         if(ns.lt.nsave) go to 5
         do 34 j=1,jmx
            mz(k(j),1)=in(j)
            mz(k(j),2)=jn(j)
            mz(k(j),3)=kn(j)
            imz(k(j))=izn(j)
 34      continue
 5       continue
         read(3,70,iostat=io) line
*        
 6    continue
*
      write(6,10) ncount, nsave
 10   format(
     + tr1,i3,' atoms read from save number ',i3,' on file:ecldsv.dat')
      return
 91   write(6,*) ' File: ecldsv.dat  could not be opened on unit 3.'
      return
 92   continue
      if(io.lt.0) then
         write(6,*) ' Unexpected end-of-file on file:  ecldsv.dat .'
      else
         write(6,*) ' Error in access to file:  ecldsv.dat .'
         write(6,*) ' IO status is ', io,' nsaves',  nsaves
      endif
      return
 30   format(i4,l2,i4,i4)
 40   format(tr1,i3,3f25.18)
 50   format(tr3,15(i3,i2))
 60   format(5(tr1,4i3,i2))
 70   format(a10)
 2500 format('   INPUT OF SAVE  '/
     +'   Each time you have given the command  osa (or  s, when'/
     +' there are entered atoms in the program) the contents of'/
     +' the arrays of Evclid are written to the end of the file'/
     +' evcldsv.dat. If the program is completely empty the command'/
     +' s    means   isa  and the last  save is automatically read'/
     +' into the program. If you give the command   isa   when the'/
     +' program is not empty you get current prompt. This gives you'/
     +' the option to read in an earlier save than the last one.'/
     +' Just enter the number of the desired save. The number of the'/
     +' last one is shown in the prompt line.'//
     +'   WARNING: All current memory contents will be removed before'/
     +' the save is read in. If you might need it do  osa  first!'/)
      end
*=====================================================================
      subroutine icfcom(str,ifc,nx,ifx,nq,a,ifz,mxn)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifz,ifh,fex,ifc
      character str*(*), atnm*2,prmp*79,filnm*80,atmnme*2
      dimension a(nx,3),nq(nx),ivec(10),fnum(10)
      prmp=
     + ' Type: file-name         for coordinate input  (h = help)'
      prmp(74:79)='f*icf*'
 70   write(6,*) prmp
      str='file'
      call readin(ifc,str,last,ifh)
      if(last.eq.0) return
      if(ifh.and.last.eq.1) then
         write(6,2500)
         go to 70
      endif
      filnm=str(1:last)
      inquire(file=filnm,exist=fex)
      if(.not.fex) then
         write(6,*) ' File does not exist. Try another name.'
         go to 70
      endif
*
      nf=1
      nl=nx
      if(mxn.gt.0) then
 90      write(6,*)
     +        ' Enter:nf ( nl)       for atoms to be re',
     +        'ad from file     (h = help)       *icf*'
         str='icf>'
         call readin(ifc,str,last,ifh)
         if(last.eq.0.or.ifc) return
         if(ifh) then
            write(6,3000)
            go to 90
         endif
         call intstr(str(1:last),nint,ivec)
         if(nint.eq.0.or.nint.gt.2) then
            write(6,*) ' Range input wrong, try again.'
            go to 90
         endif
         nf=ivec(1)
         if(nint.eq.2) nl=ivec(2)
         if(fls(nx,nf).or.fls(nx,nl).or.nf.gt.nl) then
            write(6,*) ' Range input wrong, try again.'
            go to 90
         endif
         if(noar(nx,ifx,nf,nl).ne.0) then
            write(6,*)
     +        ' Warning!  Old contents may be over-written by file.'
            write(6,*)
     +        ' Type:y (for yes) to proceed, press return to cancel.'
            str='*y?>'
            call readin(ifc,str,last,ifh)
            if(index(str,'y').eq.0.or.last.eq.0.or.ifc) return
         endif
      endif
      open(1,file=filnm,err=93,status='old')
      rewind 1
      ncount=0
      do 2 n=nf,nl
 1       continue
         str='   '
         read(1,2000,iostat=io) str
 2000    format(a)
         call strip(str,last)
         if(io.ne.0.and.last.eq.0) go to 4
         call chtflp(str(1:last),nuno,fnum)
         if(nuno.lt.3) then
            write(6,*)' Line without coordinates encountered.'
            go to 1
         endif
         atnm='x '
         if(str(1:1).ge.'A'.and.str(1:1).le.'Z') atnm(1:1)=str(1:1)
         if(str(2:2).ge.'a'.and.str(2:2).le.'z') atnm(2:2)=str(2:2)
         if(atnm.eq.'x ') then
            do 3 i=-9,-1
               if(str(1:1)//' '.eq.atmnme(i)) atnm=atmnme(i)
 3          continue
         endif
         if(atnm.eq.'x ') then
            if(nuno.eq.3) nqn=1
            if(nuno.eq.4.or.nuno.gt.5) nqn=fnum(1)
            if(nuno.eq.5) nqn=fnum(2)
            if(nqn.lt.-9.or.nqn.gt.99) nqn=1
         else
            nqn=nqq(atnm)
            if(nqn.eq.100) nqn=1
         endif
         nq(n)=nqn
         if(nqn.ne.1) ifz=.true.
         icc=0
         if(nuno.gt.3) icc=nuno-3
         a(n,1)=fnum(1+icc)
         a(n,2)=fnum(2+icc)
         a(n,3)=fnum(3+icc)
         ifx(n)=.true.
         ncount=ncount+1
 2    continue
 4    continue
      write(6,1000) ncount,filnm
 1000 format(1x,i5,' atoms read on file:',tr1,a40)
      close (1)
      if(io.eq.0) return
      if(io.eq.-1.and.ncount.eq.0) 
     +     write(6,*) ' End of file found on coordinate file.'
      if(io.ne.-1) write(6,*) ' Error in access to coordinate file.'
      return
 93   write(6,*) ' File could not be opened on logical unit 1.'
      return
 3000    format(' The input is placed in the entered range nf,nl.'/
     +' If nl is not given the file is read until end of file or'/
     +' until there is no more room in the central atom store.'/
     +'   Do not place the input where you already have entered'/
     +' atoms since then you must confirm each time that you wish'/
     +' old input over written.')
 2500 format(
     +'   ON THE COORDINATE INPUT FILE '/
     +' The file is assumed to contain Cartesian coordinates of atoms.'/
     +' The format may vary but the last three numbers in each line'/
     +' must be the x, y, and z-coordinates respectively. Leftmost in'/
     +' the line information about atomic number is looked for. A'/
     +' chemical symbol Xx must appear first in the input line if it'/
     +' is to be correctly read. If no chemical symbol is found but'/
     +' there are more numbers than three the first number is taken to'/
     +' be the atomic number Z except for the case of five numbers in'/
     +' which case the second is assumed to be Z. Thus  icf  can read'/
     +' files produced with  ocf  (and many others, possibly after'/
     +' editing to one of the forms mentioned above).'/
     +'   Note that exponential notation cannot be read by icf.')
      end
*======================================================================
      subroutine izmcom(str,ifc,nx,ifx,nq,a,ifz,
     +     imz,mz,mxn,symbs,nxz,consts)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifz,ifh,fex,ifc,ifsymb
      character str*(*),prmp*79,filnm*80, line*100
      character symbs(nxz)*20,consts(nxz)*25,strr*60
      dimension flv(10), v1(3),v2(3),v3(3)
      dimension nq(nx),a(nx,3),ivec(10),imz(nx),mz(nx,3)
      common/pp/conv
      prmp=
     +' Type: file-name         for z-matrix input   (h = help)'
      prmp(74:79)='f*izm*'
 80   write(6,*) prmp
      str='file'
      call readin(ifc,str,last,ifh)
      if(last.eq.0) return
      if(ifh.and.last.eq.1) then
         write(6,2000)
         go to 80
      endif
      filnm=str(1:last)
      inquire(file=filnm,exist=fex)
      if(.not.fex) then
         write(6,*) ' File does not exist. Try another name.'
         go to 80
      endif
*      
      nf=1
      nl=nx
      if(mxn.gt.0) then
 90      write(6,*)
     +        ' Enter: nf (nl)    for input atoms from ',
     +        'z-matrix-file    (h = help)       *izm*'
         str='izm>'
         call readin(ifc,str,last,ifh)
         if(last.eq.0.or.ifc) return
         if(ifh) then
            write(6,4000)
            go to 90
         endif
         call intstr(str(1:last),nun,ivec)
         if(nun.eq.0.or.nun.gt.2) then
            write(6,*) ' Range input wrong, try again.'
            go to 90
         endif
         nf=ivec(1)
         if(nun.eq.2) nl=ivec(2)
         if(fls(nx,nf).or.fls(nx,nl).or.nf.gt.nl) then
            write(6,*) ' Range input wrong, try again.'
            go to 90
         endif
         if(noar(nx,ifx,nf,nl).ne.0) then
            write(6,*)
     +        ' Warning!  Old contents may be over-written by file.'
            write(6,*)
     +        ' Type:y (for yes) to proceed, press return to cancel.'
            str='*y?>'
            call readin(ifc,str,last,ifh)
            if(index(str,'y').eq.0.or.last.eq.0.or.ifc) return
         endif
      endif
      nadd=nf-1
      nct=0
      open(1,file=filnm,err=93,status='old')
      rewind 1
*
* Start loop which checks for possible symbolic z-matrix
*
      ifsymb=.false.
      do 1 n=1,nx+nxz+2
         str=' '
 999     read(1,1000,iostat=io) str
 1000    format(a80)
         call strip(str,last)
         if(last.eq.0.and.io.ne.0) go to 22
         if(n.eq.1.and.last.eq.0) go to 999
         if(n.gt.1.and.last.eq.0) then
            ifsymb=.true.
            go to 22
         endif
 1    continue
 22   continue
      natoms=n-1
      if(natoms.eq.0) then
         Write(6,*) ' No atoms could be read from this file.'
         close(1)
         go to 80
      endif
      nl=min0(nf+natoms-1,nl)
      nsymb=0
      if(ifsymb) then
 23      continue
         str='    '
         read(1,1000,iostat=io) str
         call strip(str,last)
         if(last.eq.0.and.io.ne.0) go to 50
         if(last.lt.3) then
            if(last.eq.0) go to 50
            go to 92
         endif
         ind=index(str,' ')
         if(ind.gt.21.or.nsymb.ge.nxz) go to 92
         nsymb=nsymb+1
         symbs(nsymb)=str(1:ind-1)
         strr=str(ind+1:last)
         call strip(strr,last)
         consts(nsymb)=strr
         if(io.ne.0) go to 50
         go to 23
      endif
 50   continue
      if(nsymb.eq.0) ifsymb=.false.
*
* Start loop to store cartesian coordinates
*
      rewind 1
      do 2 n=nf,nl
         lin=n-nf+1
         str='  '
 998     read(1,1000,iostat=io) str
         line = str
         call strip(line,last)
         if(last.eq.0.and.io.ne.0) go to 4
         if(n.eq.nf.and.last.eq.0) go to 998
         if(ifsymb) then
            do 24 i=1,nsymb
               call strip(symbs(i),les)
               call strip(consts(i),lec)
               call replac(line,symbs(i)(1:les),consts(i)(1:lec),nocc)
 24         continue
 25         call strip(line,last)
            ind=index(line(1:last),'     ')
            if(ind.gt.1) then
               line=line(1:ind-1)//line(ind+3:last)
               go to 25
            endif
            write(6,*) line(1:last)
         endif
         call chtflp(line(1:last),nun,flv)
         nq(n)=nqq(line(1:2))
         if(nq(n).ne.1) ifz=.true.
         if(nun.eq.0) then
            if(lin.ne.1) go to 5
            do 30 i=1,3
 30            a(n,i)=0.D0
            ifx(n)=.true.
            nct=nct+1
         elseif(nun.eq.2) then
            if(lin.ne.2) go to 5
            a(n,1)=flv(2)
            a(n,2)=0.D0
            a(n,3)=0.D0
            ifx(n)=.true.
            nct=nct+1
         elseif(nun.eq.4) then
            if(lin.ne.3) go to 5
            n1=flv(1)+nadd
            n2=flv(3)+nadd
            if(nadd.eq.0) then               
               imz(n)=1
               mz(n,1)=n1
               mz(n,2)=n2
            endif
            cc=dcos(conv*flv(4))
            ss=dsin(conv*flv(4))
            sig=1.D0
            if(n1.gt.n2) sig=-sig
            a(n,1)=a(n1,1)+sig*cc*flv(2)
            a(n,2)=ss*flv(2)
            a(n,3)=0.D0
            ifx(n)=.true.
            nct=nct+1
         else
            if(nun.lt.6.or.lin.lt.4.or.nun.gt.7) go to 5
            jj=0
            if(nun.eq.7) jj=flv(7)
            n1=flv(1)+nadd
            n2=flv(3)+nadd
            n3=flv(5)+nadd
            ierr=0
            call refsys(nx,a,n1,n2,n3,v1,v2,v3,sc,d1,d2,ierr)
            if(ierr.ne.0) go to 5
            cc=dcos(conv*flv(4))
            ss=dsin(conv*flv(4))
            ccp=dcos(conv*flv(6))
            ssp=dsin(conv*flv(6))
            imz(n)=1
            if(jj.ne.0) then
               if(iabs(jj).ne.1) go to 5
               imz(n)=-1
               psi=gprm(nx,a,n3,n1,n2,0)
               cps=dcos(conv*psi)
               sps=dsin(conv*psi)
               flv(6)=0.D0
               an=ss*sps
               if(dabs(an).lt.1.D-10) go to 70
               tlj=ccp-cc*cps
               flv(6)=tlj/an
               if(dabs(flv(6)).gt.1.D0) flv(6)=dint(flv(6))
               flv(6)=-jj*dacos(flv(6))/conv
 70            jj=0
               ccp=dcos(conv*flv(6))
               ssp=dsin(conv*flv(6))
            endif
            do 40 i=1,3
               mz(n,i)=flv(2*i-1)+nadd
 40            a(n,i)=a(n1,i)+flv(2)*(cc*v1(i)+ss*(ccp*v2(i)-ssp*v3(i)))
            ifx(n)=.true.
            nct=nct+1
         endif
 2    continue
      go to 6
 4    if(io.eq.-1) write(6,*) ' End of file found on z-matrix file.'
      if(io.ne.-1) write(6,*) ' Error in access to z-matrix file.'
      go to 6
 5    write(6,3000) lin
 3000 format(' Cannot interpret z-matrix on line:',i2)
 6    close (1)
      if(nct.eq.0) return
      call clean(nx,ifx,a,nf,nf+nct-1)
      write(6,5000) nct,filnm
 5000 format(tr2,i3,' atoms read from z-matrix file: ',a40)
      return
 93   write(6,*) ' File could not be opened on unit 1.'
      return
 92   write(6,3001) nsymb+1
 3001 format(
     +' Maximum symbol length or number exceeded on line ',i2,
     +'  of symbol list.')
      return
 2000 format(
     +'        ON Z-MATRIX FILE INPUT'/
     +'   A file produced with: ozm  can be read with this option'/
     +' as well as the z-matrix part of an input file to Gaussian 90,'/
     +' even a symbolic one. There should be no other lines in the'/
     +' file; the first line should contain a single chemical symbol'/
     +' and there should be a blank line before the symbol value'/
     +' assignments.'/
     +'   For an explanation of the z-matrix concept see: dzm and'/
     +' then: h . The Cartesian coordinates are calculated from'/
     +' the information in the z-matrix and then stored in the'/
     +' central atom store. The first atom is placed at the origin,'/
     +' the second on the positive x-axis, the third in the xy-plane'/
     +' with positive y-coordinate. The remaining atoms then have'/
     +' uniquely determined positions.The z-matrix connectivity, i.e.'/
     +' the pattern of reference atoms, is also stored; see  ezm.')
 4000 format( 
     +'        INPUT FROM Z-MATRIX FILE'/
     +'   The input is placed in the entered range   nf, nl.'/
     +' If  nl  is not given the file is read until end of file'/
     +' or until the central atom store is full. If old input is'/
     +' in danger of being over written you will be warned once.'/
     +'   When nl is greater than 1 the z-matrix connectivity'/
     +' for the third line atom of the file will not be stored'/
     +' since it is not consistent to have only two reference'/
     +' atoms for an atom which is not number 3.')
      end
*======================================================================
* DISPLAY  ROUTINES:    DCOCOM, DZMCOM, DGECOM
* OUTPUT ROUTINES:      OZMCOM, OCFCOM, OSACOM
*====================================================================
      subroutine dcocom(str,ifc,nx,ifx,nq,nnq,a,mxn,ifz,lnmx)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifz,ifh,ifc
      character str*(*),atmnme*2,iq*2,ct*2,prmp*79,symb*5
      integer nq(nx),nnq(nx),inum(10)
      dimension a(nx,3)
*
      prmp=
     +' Enter:(ct)   (nf (nl))                 '//
     +'        (h=help, q=skip)          *dco*'
 90   write(6,*) prmp
      str='dco>'
      call readin(ifc,str,last,ifh)
      if(ifc) return
      it=1
      nf=1
      nl=mxn
      nuno=0
      if(last.eq.0) go to 60
      if(ifh) then
         write(6,1000) 
         go to 90
      endif
      call intstr(str(1:last),nuno,inum)
      if(nuno.gt.2) go to 91
      ct=str(1:2)
      it=0
      if(ct.eq.'x '.or.ct.eq.'  ') it=1
      if(ct.eq.'c ') it=2
      if(ct.eq.'s ') it=3
      if(nuno.eq.0) then
         if(ct.eq.'q ') return
      elseif(nuno.eq.1) then
         if(it.eq.0) it=1
         nf=inum(1)
      elseif(nuno.eq.2) then
         if(it.eq.0) it=1
         nf=inum(1)
         nl=inum(2)
      endif
      if(it.eq.0) go to 91
      if(nf.eq.0) return
      if(fls(nx,nf).or.fls(nx,nl)) go to 91
      if(nf.gt.nl) go to 91
 60   continue
      if(ifz) then
         if(it.eq.1) write(6,3)
         if(it.eq.2) write(6,4)
         if(it.eq.3) write(6,5)
      else
         if(it.eq.1) write(6,30)
         if(it.eq.2) write(6,40)
         if(it.eq.3) write(6,50)
      endif
 3    format(' Symbol   Z            n ',6x,'x',tr15,'y',tr15,'z')
 4    format(' Symbol   Z            n ',5x,'rh',13x,'phi',tr15,'z')
 5    format(' Symbol   Z            n ',6x,'r',13x,'tht',13x,'phi')
 30   format('                n ',6x,'x',tr15,'y',tr15,'z')
 40   format('                n ',5x,'rh',13x,'phi',tr15,'z')
 50   format('                n ',6x,'r',13x,'tht',13x,'phi')
      ncou=1
      do 6 i=nf,nl
         if(.not.ifx(i)) go to 6
         x=a(i,1)
         y=a(i,2)
         z=a(i,3)
         if(it.eq.2) then
            call cyli(x,y,r,phi,-1)
            x=r
            y=phi
         elseif(it.eq.3) then
            call spheri(x,y,z,r,phi,the,-1)
            x=r
            y=the
            z=phi
         endif
*
*
         if(ifz) then
            iq=atmnme(nq(i))
            write(6,7) symb(iq,nnq(i)),nq(i),i,x,y,z
         else
            write(6,8) i,x,y,z
         endif
 7       format(tr2,a5,tr2,i2,tr10,i3,3(tr1,f15.8))
 8       format(tr14,i3,3(1x,f15.8))
         ncou=ncou+1
         if(i.eq.mxn) return
         if(ncou.eq.lnmx.and.nl.eq.mxn) then
            write(6,*)' Press return for more or Type:q to quit'
            read(5,2000,iostat=io) str
 2000       format(a)
            if(index(str,'q ').ne.0) return
            if(ifz) then
               if(it.eq.1) write(6,3)
               if(it.eq.2) write(6,4)
               if(it.eq.3) write(6,5)
            else
               if(it.eq.1) write(6,30)
               if(it.eq.2) write(6,40)
               if(it.eq.3) write(6,50)
            endif
            ncou=1
         endif
 6    continue
      return
 91   write(6,*) ' Input wrong to coordinate display; type:q to cancel.'
      go to 90
 1000    format('        COORDINATE DISPLAY '//
     +' For  Cartesian  coordinates, just push return; '/
     +' for Cylindrical coordinates, type:c (and then return); '/
     +' for  Spherical  coordinates, type:s (and then return). '//
     +'   If you wish to look at (display) a sub-range nf, nl'/
     +' then enter it. '//
     +'   To cancel (skip) display  type:q.'//
     +'         DEFINITIONS for Spherical and Cylindrical coordinates:'/
     +' rh =  distance to z-axis (rh stands for the Greek letter rho).'/
     +' phi = angle between positive x-axis and the projection of the'/
     +'       point''s position vector in the xy-plane. phi is between'/
     +'       -180 and 180 degrees and has the same sign as the y-'/
     +'       coordinate of the projection.'/
     +' r =   distance from point to origin'/
     +' tht = angle between positive z-axis and the point''s position'/
     +'       vector (between 0 and 180 degrees). (tht stands for the'/
     +'       Greek letter theta.)'/)
      end
*=====================================================================
      subroutine dzmcom(str,ifc,nx,ifx,nq,nnq,a,mxn,ifz,lnmx,nmz,imz,mz)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifz,alt,ord,cal,ifh,ifc,mix,fxd,arenot
      character str*(*),atmnme*2,isym*2,prmp*79,symb*5,sgn*2
      dimension inum(10),v1(3),v2(3),v3(3),nq(nx),nnq(nx),a(nx,3)
      integer imz(nx),mz(nx,3)
*
* First find the numbers n1, n2, n3 of the first three entered atoms
      call trefst(nx,ifx,n1,n2,n3)
      nil=0
 90   prmp=
     +' Enter: ( nf ( nl ( i j k)))   (a) (o) ('//
     +'#)    (h=help, q=skip)            *dzm*'
      write(6,*) prmp
      str='dzm>'
      call readin(ifc,str,last,ifh)
      if(ifc.or.str(last:last).eq.'q') return
      nun=0
      if(last.eq.0) go to 92
      if(ifh) then
         if(str(1:last).ne.'hh') then
            write(6,1000)
            write(6,1001)
         else
            write(6,1002)
         endif
         go to 90
      endif
      call intstr(str(1:last),nun,inum)
 92   if(n3.eq.0) then
         write(6,*) ' Z-matrix requires three reference atoms.'
         return
      endif
      if(nun.eq.3.or.nun.eq.4.or.nun.gt.5) go to 91
      nf=1
      if(nun.ge.1) nf=inum(1)
      nl=mxn
      if(nun.ge.2) nl=inum(2)
      if(fls(nx,nf).or.fls(nx,nl).or.nl.lt.nf) then
         write(6,*) ' Range input (nf,nl) wrong.'
         go to 90
      endif
      fxd=.false.
      alt=.false.
      ord=.false.
      cal=.false.
      if(index(str,'a').ne.0) alt=.true.
      if(index(str,'o').ne.0) ord=.true.
      if(index(str,'#').ne.0) cal=.true.
      if(alt.and.ord) go to 91
      if(nun.ge.5) then
         n1=inum(3)
         n2=inum(4)
         n3=inum(5)
         fxd=.true.
         if(arenot(nx,ifx,n1,n2,n3)) then
            write(6,*) ' Entered reference atom numbers not accepted.'
            go to 90
         endif
      endif
      ierr=0
      call refsys(nx,a,n1,n2,n3,v1,v2,v3,vin,d1,d2,ierr)
      if(ierr.ge.1) then
         if(fxd) go to 90
         return
      endif
      mix=.false.
      if(.not.(fxd.or.alt.or.ord.or.nmz.eq.0)) then
         itest=0
         jtest=0
         do 20 i=max(4,nf),nl
            if(ifx(i).and.imz(i).eq.1) itest=itest+1
            if(ifx(i).and.imz(i).eq.-1) jtest=jtest+1
 20      continue
         if(itest*jtest.ne.0) mix=.true.
         if(itest.eq.0.and.jtest.ne.0.and.(.not.(ord.or.alt)))mix=.true.
         if(nf.eq.nl) mix=.false.
      endif
      if(nf.eq.nl.and.ifx(nf).and.imz(nf).eq.-1.and.(.not.ord))
     +     alt=.true.
      if(ifz.and.(.not.mix)) then
         if(.not.alt) write(6,2)
         if(alt) write(6,3)
      elseif(.not.mix) then
         if(.not.alt) write(6,12)
         if(alt) write(6,13)
      else
         write(6,5)
      endif
 2    format(' Symbol  Z      n,   i,   r(n,i)   ,',
     +'  j, tht(n,i,j),  k, phi(n,i,j,k)')
 3    format(' Symbol  Z      n,   i,   r(n,i)   ,',
     +'  j, tht(n,i,j),  k, tht(n,i,k),     sign')
 5    format(' Symbol  Z      n,   i,   r(n,i)   ,',
     +'  j, tht(n,i,j),  k, ang(n,i,(j,) k) sign')
 12   format('                n,   i,   r(n,i)   ,',
     +'  j, tht(n,i,j),  k, phi(n,i,j,k)')
 13   format('                n,   i,   r(n,i)   ,',
     +'  j, tht(n,i,j),  k, tht(n,i,k),     sign')
* Major write-out do-loop starts......................................
      ncount=0
      ncou=1
      n1st=n1
      n2st=n2
      n3st=n3
      do 4 in=nf,nl
         if(.not.ifx(in)) go to 4
         ncount=ncount+1
         if(fxd.or.in.eq.n1st.or.in.eq.n2st) go to 10
         if(imz(in).eq.0.or.cal) then
            call refnuc(nx,ifx,nq,a,n3st,in,n1,n2,n3)
         else
            n1=mz(in,1)
            n2=mz(in,2)
            if((.not.ifx(n1)).or.(.not.ifx(n2))) go to 93
            if(in.ge.4) then
               n3=mz(in,3)
               if(.not.ifx(n3)) go to 93
            endif
         endif
 10      bl=gprm(nx,a,in,n1,nil,nil)
         tht=gprm(nx,a,in,n1,n2,nil)
         phi=gprm(nx,a,in,n1,n2,n3)
         sgn='  '
         if(alt.or.(imz(in).eq.-1.and.(.not.ord))) then
            jt=-idint(dsign(1.D0,phi))
            sgn=' 1'
            if(jt.eq.-1) sgn='-1'
            phi=gprm(nx,a,in,n1,n3,nil)
         endif
         fii=phi
*
         if(ifz) then
            isym=atmnme(nq(in))
            if(in.eq.n1) then
               write(6,80) symb(isym,nnq(in)),nq(in),in
            elseif(in.eq.n2) then
               write(6,81) symb(isym,nnq(in)),nq(in),in,n1,bl
            elseif(in.eq.n3) then
               write(6,82) symb(isym,nnq(in)),nq(in),in,n1,bl,n2,tht
            else
               write(6,83) 
     +              symb(isym,nnq(in)),nq(in),in,n1,bl,n2,tht,n3,fii,sgn
            endif
         else
            if(in.eq.n1) then
               write(6,60) in
            elseif(in.eq.n2) then
               write(6,61) in,n1,bl
            elseif(in.eq.n3) then
               write(6,62) in,n1,bl,n2,tht
            else
               write(6,63) in,n1,bl,n2,tht,n3,fii,sgn
            endif
         endif
         ncou=ncou+1
         if(in.eq.mxn) return
         if(ncou.eq.lnmx.and.nl.eq.mxn) then
            write(6,*)' Press return for more or Type:q to quit'
            read(5,2000,iostat=io) str
 2000       format(a)
            if(index(str,'q ').ne.0) return
            if(ifz) then
               if(.not.alt) write(6,2)
               if(alt) write(6,3)
            else
               if(.not.alt) write(6,12)
               if(alt) write(6,13)
            endif
            ncou=1
         endif
* 
 4    continue
* Write do-loop ends  here .............................
      return
 91   write(6,*) ' Input to Z-matrix display wrong.'
      go to 90
 93   write(6,*) ' Stored z-matrix information is wrong.'
      return
*----------------------------------------Below only format statements-
 80   format(tr2,a5,tr1,i2,4x,i3)
 81   format(tr2,a5,tr1,i2,4x,i3,i5,f12.7)
 82   format(tr2,a5,tr1,i2,4x,i3,i5,f12.7,2x,i3,f11.5)
 83   format(
     +tr2,a5,tr1,i2,4x,i3,i5,f12.7,2x,i3,f11.5,2x,i3,f12.5,tr7,a2)
 60   format(tr14,i3)
 61   format(tr14,i3,i5,f12.7)
 62   format(tr14,i3,i5,f12.7,2x,i3,f11.5)
 63   format(tr14,i3,i5,f12.7,2x,i3,f11.5,2x,i3,f12.5,tr7,a2)
 1000 format('           ON THE Z-MATRIX '/
     +'   A z-matrix is a set of 3N-6 geometric parameters defining'/
     +' the shape (geometry) of an N-atom molecule in a coordinate-'/
     +' independent way. The positions of the atoms (in range nf, nl)'/
     +' are specified in terms of the three geometric parameters:'/
     +'       distance to atom i,                         r(n,i)'/
     +'       angle between the bond to i and bond i-j,   tht(n,i,j)'/
     +'       dihedral angle along i-j to plane i-j-k,    phi(n,i,j,k)'/
     +' and a typical line is:'/
     +' (Z-info)   n   i   r(n,i)   j   tht(n,i,j)   k   phi(n,i,j,k)'/
     +'   No input is needed, just press return and all entered atoms'/
     +' are listed and, unless there is stored z-matrix info, i,j,k'/
     +' are choosen automatically according to a nearest neighbour'/
     +' principle among atoms with numbers less than  n.')
 1001 format(
     +'   As further input one can give: nf,nl (to specify the range)'/
     +' and (if this has been done) also the desired reference atoms:'/
     +' i,j,k (in case the automatic choice is unsuitable).'/
     +'   By giving an: a (for alternative) in the input line the di-'/
     +' hedral angle is replaced by the angle: tht(n,i,k). A sign con-'/
     +' vention is then needed to make the position unique and it is'/
     +' flagged by a 1 or a -1 (rightmost in the line). For sign con-'/
     +' ventions see  ige  and then: hs .           (hh=more help)')
 1002 format(
     +'   If there is stored z-matrix information (reference atoms),'/
     +' then these stored reference atoms will be used in the dis-'/
     +' play. Such information is stored when the program reads a z-'/
     +' matrix file (command  izm). Apart from the reference atoms'/
     +' the type of the specification: ordinary (o) or alternative (a)'/
     +' (i.e. use second bond angle instead of dihedral angle) is also'/
     +' stored. This means that if you read a z-matrix with  izm  and'/
     +' then output it with  ozm,  or display it, it will be exactly'/
     +' as read.  By putting an  a  or  an  o  in the input line you'/
     +' make the display to be of the indicated type independently of'/
     +' stored information. By putting the character  #  in the input'/
     +' line you will get the calculated (automatic) reference atoms.'/
     +' (Thus input of  o #   or   a #  results in a display ignoring'/
     +' stored information.)'/
     +'   Note that you can  input, display  and remove  stored z-'/
     +' matrix reference atoms (and type); give the command  ezm.')
      end
*====================================================================
      subroutine dgecom(str,ifc,nx,ifx,nq,nnq,a,mxn,ifz,lnmx)
      implicit double precision(a-h,o-z)
      logical ifx(nx),ifz,fls,ifh,ifix,ifc
      double precision gprm,a(nx,3),fnum(10),v1(3)
      integer nq(nx),nnq(nx),ivec(10)
      character str*(*),sstr*23,prmp*79
      nil=0
      prmp=' Enter:(-i) (m (n))   (< max-dist.)   (%'//
     +' nf nl)     (h=help, q=skip)'
      prmp(75:79)='*dge*'
 90   write(6,*) prmp
      str='dge>'
      call readin(ifc,str,last,ifh)
      if(ifc.or.str(last:last).eq.'q') return
      if(ifh) then
         write(6,1000) 
         write(6,1001) 
         go to 90
      endif
      nf=1
      nl=mxn
      dm=1.D20
      nint=0
      if(last.eq.0) go to 1
      idpos=index(str(1:last),'<')
      irpos=index(str(1:last),'%')
      if(irpos.ne.0) then
         if(irpos.lt.idpos) go to 91
         if(irpos.eq.last) go to 91
         call intstr(str(irpos+1:last),nint,ivec)
         if(nint.lt.1.or.nint.gt.2) go to 91
         nf=ivec(1)
         if(nint.eq.2) nl=ivec(2)
         if(fls(nx,nf).or.fls(nx,nl).or.nf.gt.nl) go to 91
         last=irpos
      endif
      if(idpos.ne.0) then
         if(idpos.eq.last) go to 91
         call chtflp(str(idpos+1:last),nun,fnum)
         if(nun.eq.0.or.nun.gt.1) go to 91
         dm=fnum(1)
         if(dm.le.0.) go to 91
         last=idpos
      endif
      call intstr(str(1:last),nint,ivec)
 1    m=0
      n=0
      nfi=nf
      nli=nl
      if(nint.ge.1) m=ivec(1)
      if(nint.ge.2) n=ivec(2)
      ifix=.false.
      if(m.lt.0) then
         ifix=.true.
         nint=nint-1
         nfi=iabs(m)
         if(fls(nx,nfi)) go to 91
         if(.not.ifx(nfi)) go to 91
         m=0
         n=0
         if(nint.ge.1) m=ivec(2)
         if(nint.ge.2) n=ivec(3) 
         nli=nfi+1
      endif
      if(nint.gt.2) go to 91
      it=0
      if(m.eq.0.and.n.eq.0) it=1
      if(m.gt.0.and.n.eq.0) it=2
      if(m.gt.0.and.n.gt.0) it=3
      if(it.eq.0.or.m.gt.nx.or.n.gt.nx) go to 91
      if(it.gt.1.and..not.ifx(m)) go to 91
      if(it.eq.3.and..not.ifx(n)) go to 91
      if(it.gt.1.and.m.eq.n) go to 91
      ierr=0
      if(it.eq.3) call bij(nx,a,m,n,v1,dd,ierr)
      if(ierr.gt.0) go to 90
      ncou=1
      if(it.eq.1) then
         if(ifz) write(6,100)
         if(.not.ifz) write(6,101)
      elseif(it.eq.2) then
         if(ifz) write(6,200)
         if(.not.ifz) write(6,201)
      elseif(it.eq.3) then
         if(ifz) write(6,300)
         if(.not.ifz) write(6,301)
      endif
 100  format(
     +'    i    j          r(i,j)         charges   symbols')
 200  format(
     +'    i    m    j        tht(i,m,j)      charges    symbols')
 300  format(
     +'    i    m    n    j    phi(i,m,n,j)     charges        symbols')  
 101  format('    i    j          r(i,j)   ')
 201  format('    i    m    j        tht(i,m,j)  ')
 301  format('    i    m    n    j    phi(i,m,n,j)  ')
*-----------------------------------------------------------------
      do 3 i=nfi,nli-1
         if(.not.ifx(i)) go to 3
         njf=i+1
         if(ifix) njf=nf
         do 2 j=njf,nl
            if(ifix.and.i.eq.j) go to 2
            if(.not.ifx(j)) go to 2
            if(it.gt.1) go to 12
            d1=gprm(nx,a,i,j,nil,nil)
            if(d1.gt.dm) go to 2
            go to 50
 12         if(i.eq.m.or.j.eq.m) go to 2
            d1=gprm(nx,a,i,m,nil,nil)
            d2=gprm(nx,a,j,m,nil,nil)
            if(it.eq.3) go to 13
            if(d1.gt.dm.or.d2.gt.dm) go to 2
            tht=gprm(nx,a,i,m,j,nil)
            go to 50
 13         if(i.eq.n.or.j.eq.n) go to 2
            d3=gprm(nx,a,i,n,nil,nil)
            d4=gprm(nx,a,j,n,nil,nil)
            if(d1.gt.dm.and.d3.gt.dm) go to 2
            if(d2.gt.dm.and.d4.gt.dm) go to 2
            phi=gprm(nx,a,i,m,n,j)
 50         continue
            if(ncou.ge.lnmx) then
               write(6,*)' Press return for more or Type:q to quit'
               read(5,2000,iostat=io) str
 2000          format(a30)
               if(index(str,'q ').ne.0) return
               if(it.eq.1) write(6,100)
               if(it.eq.2) write(6,200)
               if(it.eq.3) write(6,300)
               ncou=1
            endif
            if(ifz) then
               if(it.eq.1) then
                  call symstr(nx,nq,nnq,i,j,nil,nil,sstr)
                  write(6,10) i,j,d1,nq(i),nq(j),sstr
               elseif(it.eq.2) then
                  call symstr(nx,nq,nnq,i,m,j,nil,sstr)
                  write(6,20) i,m,j,tht,nq(i),nq(m),nq(j),sstr
               elseif(it.eq.3) then
                  call symstr(nx,nq,nnq,i,m,n,j,sstr)
                  write(6,30) i,m,n,j,phi,nq(i),nq(m),nq(n),nq(j),sstr
               endif
            else
               if(it.eq.1) write(6,11) i,j,d1
               if(it.eq.2) write(6,21) i,m,j,tht
               if(it.eq.3) write(6,31) i,m,n,j,phi
            endif
            ncou=ncou+1
 2       continue
 3    continue
      return
 10   format(tr2,i3,tr2,i3,tr4,f17.9,tr5,i2,'-',i2,tr4,a23)
 11   format(tr2,i3,tr2,i3,tr4,f17.9)
 20   format(tr2,i3,tr2,i3,tr2,i3,tr4,f14.6,tr5,2(i2,'-'),i2,tr4,a23)
 21   format(tr2,i3,tr2,i3,tr2,i3,tr4,f14.6)
 30   format(4(tr2,i3),tr3,f12.6,tr4,3(i2,'-'),i2,tr4,a23)
 31   format(4(tr2,i3),tr3,f12.6)
 91   write(6,*) ' Input to geometric parameter display wrong.'
      go to 90
 1000 format(
     +'   List of GEOMETRIC PARAMETERS i.e. distances, (unsigned) ang-'/
     +' les between lines and (signed) angles between planes, running'/
     +' over (all) pairs of (entered) atoms i,j is displayed.'//
     +'   For DISTANCES r(i,j) no input is needed (just press return).'/
     +'   For ANGLES tht(i,m,j) at atom m between lines to i and j,',
     +' just enter: m'/
     +'   For DIHEDRAL (or torsion) ANGLES phi(i,m,n,j) between planes',
     +' intersecting at'/' the line from m to n, containing i and j re',
     +'spectively, enter: m, n '/)
 1001 format(
     +'   To get a list over j only, put -i first. E.g. if you  Enter:'/
     +' -5 only you will get a list over distances to atom 5.'/
     +'   To limit output it is optional to specify a RANGE to which'/
     +' i and j are limited by placing: % nf nl (e.g.: % 1 5 for the '/
     +' first five atoms) last in the input line.  '/
     +'   Output can also be limited by giving a MAX-DISTANCE. This is'/
     +' is done by putting: < max-dist. (e.g.: < 1.34698) in the input'/
     +' line, after m and n if these are needed. In the case of angles'/
     +' both distances r(i,m) and r(m,j) (or r(n,j) in the dihedral'/
     +' case) must be less than max-dist for the angle to be listed.'/
     +'   Undefined angles (the) are given as zero while undefined di-'/
     +' hedral angles (phi) are put to -360.')
      end
*=====================================================================
      subroutine ozmcom(str,ifc,nx,ifx,nq,nnq,a,mxn,
     +     imz,mz,symbs,nxz,cnzs)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifxx,alt,ifh,acu,fex,ifc,arenot,ifsymb
      character str*(*),atmnme*2,jct*2,isym*2,prmp*79,filnm*80
      dimension inum(10),v1(3),v2(3),v3(3),nq(nx),a(nx,3),cnzs(nxz)
      character*20 symbs(nxz)
      character blsy*10,thtsy*15,fiisy*20,symstz*20
      integer imz(nx),mz(nx,3),nnq(nx)
      save mem,filnm,acu,ifsymb
      data mem/0/,nil/0/
*     
 69   continue
      if(mem.eq.0) then
         prmp=
     +        ' Type: file-name (*=symb. +=accu.) for z-matrix output'//
     +        '     (h = help)'
         prmp(74:79)='f*ozm*'
 70      write(6,*) prmp
         str='file'
         call readin(ifc,str,last,ifh)
         if(last.eq.0) return
         if(ifh) then
            if(last.eq.2.and.str(1:last).eq.'hh') then 
               write(6,1000)
               write(6,1010)
            else
               write(6,1001)
            endif
            go to 70
         endif
         acu=.false.
         ifsymb=.false.
         if(str(last:last).eq.'+') acu=.true.
         if(str(last:last).eq.'*') ifsymb=.true.
         if(acu.or.ifsymb) last=last-1
         filnm=str(1:last)
         inquire(file=filnm,exist=fex)
         if(fex) then
            write(6,*) ' File already exists, choose another name.'
            go to 70
         endif
         open(4,file=filnm,err=93,status='new')
         mem=mem+1
      else
         prmp=' Output will be added to file:'//filnm
         write(6,*) prmp
      endif
*     
 90   prmp=
     +     ' Enter: ( nf ( nl ( i j k)))  (a)   (/=c'//
     +     'lose, *, +)  (h=help, q=skip)     *ozm*'
      write(6,*) prmp
      str='ozm>'
      call readin(ifc,str,last,ifh)
      if(ifc.or.str(last:last).eq.'q') return
      if(last.eq.0) go to 92
      if(last.eq.1.and.str(1:1).eq.'/') then
         close (4)
         mem=0
         write(6,*)' File closed. You can now give a new name.'
         go to 69
      endif
      if(index(str,'*').ne.0) ifsymb=.true.
      if(index(str,'+').ne.0) then
         acu=.true.
         if(ifsymb) acu=.false.
      endif
      if(ifh) then
         if(last.eq.2.and.str(1:last).eq.'hh') then
            write(6,1001)
         else
            write(6,1000)
            write(6,1010)
         endif
         go to 90
      endif
 92   call trefst(nx,ifx,n1,n2,n3)
      if(n3.eq.0) then
         write(6,*) ' Z-matrix requires three reference atoms.'
         return
      endif
      call intstr(str(1:last),nun,inum)
      if(nun.eq.3.or.nun.eq.4.or.nun.gt.5) go to 91
      nf=1
      if(nun.ge.1) nf=inum(1)
      nl=mxn
      if(nun.ge.2) nl=inum(2)
      if(fls(nx,nf).or.fls(nx,nl).or.nl.lt.nf) then
         write(6,*) ' Range input (nf,nl) wrong.'
         go to 90
      endif
      irt=0
      alt=.false.
      if(index(str,'a').ne.0) alt=.true.
      if(nun.ge.5) then
         n1=inum(3)
         n2=inum(4)
         n3=inum(5)
         irt=1
         if(arenot(nx,ifx,n1,n2,n3)) then
            write(6,*) ' Entered reference atom numbers not accepted.'
            go to 90
         endif
      endif
      ierr=0
      call refsys(nx,a,n1,n2,n3,v1,v2,v3,vin,d1,d2,ierr)
      if(ierr.ge.1) then
         if(irt.eq.1) go to 90
         return
      endif
*     Major write-out do-loop starts......................................
      n3st=n3
      ncount=0
      nvix=1
      do 4 in=nf,nl
         if(.not.ifx(in)) then
            ifxx=.false.
            if(in.lt.nl) ifxx=ifx(in+1)
            if(ifxx) write(6,400) in
 400        format(' Warning, output z-matrix ruined by gap at:',i2)
            go to 4
         endif
         ncount=ncount+1
         if(irt.ne.0.or.ncount.le.2) go to 10
         if(imz(in).eq.0) then
            call refnuc(nx,ifx,nq,a,n3st,in,n1,n2,n3)
         else
            n1=mz(in,1)
            n2=mz(in,2)
            if((.not.ifx(n1)).or.(.not.ifx(n2))) go to 94
            if(in.ge.4) then
               n3=mz(in,3)
               if(.not.ifx(n3)) go to 94
            endif
         endif
 10      bl=gprm(nx,a,in,n1,nil,nil)
         tht=gprm(nx,a,in,n1,n2,nil)
         phi=gprm(nx,a,in,n1,n2,n3)
         if(ifsymb.and.in.gt.n1) then
            if(nvix.gt.nxz) go to 95
            cnzs(nvix)=bl
            blsy=symstz(nx,nq,nnq,in,n1,nil,nil)
            symbs(nvix)=blsy
            nvix=nvix+1
            if(nvix.gt.nxz) go to 95
            if(in.gt.n2) then
               cnzs(nvix)=tht
               thtsy=symstz(nx,nq,nnq,in,n1,n2,nil)
               symbs(nvix)=thtsy
               nvix=nvix+1
            endif
            if(in.gt.n3) then
               if(nvix.gt.nxz) go to 95
               symbs(nvix)=symstz(nx,nq,nnq,in,n1,n2,n3)
            endif
         endif
         if((dabs(tht).lt.1.D-3.or.dabs(tht-180.D0).lt.1.D-3).and.
     +        in.ge.3) 
     +        write(6,*)' Warning,angle 0 or 180 degrees.'
         jt=0
         if(alt.or.imz(in).eq.-1) then
            jt=-idint(dsign(1.D0,phi))
            phi=gprm(nx,a,in,n1,n3,nil)
            if(ifsymb.and.in.gt.n3)
     +           symbs(nvix)=symstz(nx,nq,nnq,in,n1,n3,nil)
         endif
         fii=phi
         if(ifsymb.and.in.gt.n3) then
            cnzs(nvix)=phi
            fiisy=symbs(nvix)
            nvix=nvix+1
         endif
         isym=atmnme(nq(in))
*     
*     Output to logical unit 4 written here:
*     
         if(acu) then
            if(in.eq.n1) write(4,100) isym
            if(in.eq.n1) write(6,100) isym
            if(in.eq.n2) write(4,301) isym,n1,bl
            if(in.eq.n2) write(6,301) isym,n1,bl
            if(in.eq.n3) write(4,302) isym,n1,bl,n2,tht
            if(in.eq.n3) write(6,302) isym,n1,bl,n2,tht
            if(in.ne.n1.and.in.ne.n2.and.in.ne.n3) then
               if(jt.eq.0) jct='  '
               if(jt.eq.1) jct=' 1'
               if(jt.eq.-1) jct='-1'
               write(4,303) isym,n1,bl,n2,tht,n3,fii,jct
               write(6,303) isym,n1,bl,n2,tht,n3,fii,jct
               if(in.lt.n1.or.in.lt.n2.or.in.lt.n3) write(6,*)
     +' Warning:reference nucleus number greater than line number!'
            endif
         elseif(ifsymb) then
            if(in.eq.n1) write(4,100) isym
            if(in.eq.n1) write(6,100) isym
            if(in.eq.n2) write(4,201) isym,n1,blsy
            if(in.eq.n2) write(6,201) isym,n1,blsy
            if(in.eq.n3) write(4,202) isym,n1,blsy,n2,thtsy
            if(in.eq.n3) write(6,202) isym,n1,blsy,n2,thtsy
            if(in.ne.n1.and.in.ne.n2.and.in.ne.n3) then
               if(jt.eq.0) jct='  '
               if(jt.eq.1) jct=' 1'
               if(jt.eq.-1) jct='-1'
               write(4,203) isym,n1,blsy,n2,thtsy,n3,fiisy,jct
               write(6,203) isym,n1,blsy,n2,thtsy,n3,fiisy,jct
               if(in.lt.n1.or.in.lt.n2.or.in.lt.n3) write(6,*)
     +' Warning:reference nucleus number greater than line number!'
            endif            
         else
            if(in.eq.n1) write(4,100) isym
            if(in.eq.n1) write(6,100) isym
            if(in.eq.n2) write(4,101) isym,n1,bl
            if(in.eq.n2) write(6,101) isym,n1,bl
            if(in.eq.n3) write(4,102) isym,n1,bl,n2,tht
            if(in.eq.n3) write(6,102) isym,n1,bl,n2,tht
            if(in.ne.n1.and.in.ne.n2.and.in.ne.n3) then
               if(jt.eq.0) jct='  '
               if(jt.eq.1) jct=' 1'
               if(jt.eq.-1) jct='-1'
               write(4,103) isym,n1,bl,n2,tht,n3,fii,jct
               write(6,103) isym,n1,bl,n2,tht,n3,fii,jct
               if(in.lt.n1.or.in.lt.n2.or.in.lt.n3) write(6,*)
     +' Warning:reference nucleus number greater than line number!'
            endif
         endif
*     
 4    continue
*     Write do-loop ends  here .............................
      if(ifsymb) then
         write(6,*) '      '
         write(4,*) '      '
         do 5 kk=1,nvix-1
            write(6,500) symbs(kk), cnzs(kk)
            write(4,500) symbs(kk), cnzs(kk)
 5       continue
         close (4)
         mem=0
      endif
      return
 500  format(tr2,a20,tr2,f20.13)
 91   write(6,*) ' Input to Z-matrix output wrong.'
      go to 90
 93   write(6,*) ' File could not be opened on unit 4.'
      return
 94   write(6,*) ' Stored z-matrix information wrong.'
      return
 95   write(6,*) ' More than', nxz, ' symbols needed.'
      return
 96   write(6,*) ' When symbolic z-matrix nf and n1 must be the same.'
      go to 90
*----------------------------------------Below only format statements-
*     Note:output z-matrix not compatible with Gaussian76. Chemical
*     symbol must be changed to atomic number and the format must
*     also be changed. This output for Gaussian80,90  and EVCLID.
 100  format(tr1,a2)
 101  format(tr1,a2,tr1,i3,tr1,f10.7)
 102  format(tr1,a2,tr1,i3,tr1,f10.7,tr1,i3,tr1,f11.5)
 103  format(tr1,a2,tr1,i3,tr1,f10.7,2(tr1,i3,tr1,f11.5),3x,a2)
 201  format(tr1,a2,tr1,i3,tr1,a10)
 202  format(tr1,a2,tr1,i3,tr1,a10,tr1,i3,tr1,a15)
 203  format(tr1,a2,tr1,
     +i3,tr1,a10,tr1,i3,tr1,a15,tr1,i3,tr1,a20,tr2,a2)
 301  format(tr1,a2,tr1,i3,tr1,f17.11)
 302  format(tr1,a2,tr1,i3,tr1,f17.11,tr1,i3,tr1,f14.8)
 303  format(tr1,a2,tr1,i3,tr1,f17.11,2(tr1,i3,tr1,f14.8),3x,a2)
 1001 format('    ON THE  OUTPUT FILE FORMAT '/
     +' For general information about the z-matrix type: ;dzm;h .'/
     +'   The file-name for Z-MATRIX OUTPUT must be a new name.'/
     +' A star:* after the file-name results in a symbolic z-matrix.'/
     +' A plus:+ after the file name gives more digits in the nume-'/
     +' cal z-matrix, which is default. When the output is numerical'/
     +' the file, once opened, stays open until you close it. If you'/
     +' give the command:ozm again during the run you will only get a'/
     +' message saying which file is open and then the output prompt.'/
     +' What is then displayed is added to the end of the file.'/
     +'   The output prompt gives you the alternative to close the'/
     +' currently open output file so that a new one can be opened.'/
     +'   The default output FORMAT for a typical line is:'/
     +'     format(tr1,a2,tr1,i3,tr1,f10.7,2(tr1,i3,tr1,f11.5),3x,a2)'/
     +' Here the first a2-field is for the chemical symbol Xx, the'/
     +' last a2-field contains the sign for the second angle when the'/
     +' alternative z-matrix is used, otherwise it is blank.'/
     +' If you put  format = +  you get'/
     +'     format(tr1,a2,tr1,i3,tr1,f17.11,2(tr1,i3,tr1,f14.8),3x,a2)'/
     +'                                                 (hh=more help)')
 1000 format('     ON  THE  OUTPUT  Z-MATRIX '/
     +'   The positions of the atoms (in range nf,nl) are given'/
     +' in terms of the three geometric parameters as follows:'/
     +' (Z-info)   n   i   r(n,i)   j   tht(n,i,j)   k   phi(n,i,j,k)'/
     +'   No input is needed, just press return and all entered atoms'/
     +' are listed and, for each atom n,  i,j,k are either taken from'/
     +' stored information, or are choosen automatically by a nearest'/
     +' neighbour algorithm, if no such information is stored.'/
     +'   The reference atoms for atom n are choosen among atoms with'/
     +'  numbers less than n (an output z-matrix can be read with: izm'/
     +' only when this is the case). No gaps are allowed, all atoms'/
     +' from 1 and up to the last atom are needed for consistency.'/
     +'   As optional input one can give: nf,nl (to specify the range)'/
     +' and (if this has been done) also the desired reference atoms:'/
     +' i,j,k (in case the automatic or stored choice is unsuitable).'/
     +'   By giving an: a (for alternative) in the input line the di-'/
     +' hedral angle is replaced by the angle: tht(n,i,k).')
 1010 format(
     +'   An output file created by (repeatedly) using ozm is accept-'/
     +' able geometry input for the  Gaussian 80 - 90  programs.'/
     +' Use dzm and ezm first to make the desired z-matrix, then ozm.'/
     +'   To close the currently open output file type:/   You will'/
     +' then be prompted for a new output filename.    (hh=more help)')
      end
*========================================================================
      subroutine ocfcom(str,ifc,nx,ifx,nq,nnq,a,mxn,ifz)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifh,acu,aft,fex,ifc,ifz
      character str*(*),atmnme*2,iq*2,prmp*79,filnm*80,symb*5,smb*5
      integer nq(nx),nnq(nx),inum(10)
      dimension a(nx,3)
      save mem,filnm,acu,aft
      data mem/0/
*
 70   continue
      if(mem.eq.0) then
         prmp=
     +' Type: file-name  (format)         '//
     +' (h = help)'
         prmp(74:79)='f*ocf*'
 80      write(6,*) prmp
         str='file'
         call readin(ifc,str,last,ifh)
         if(last.eq.0) return
         if(ifh) then
            write(6,1001)
            go to 80
         endif
         acu=.false.
         aft=.false.
         if(str(last:last).eq.'+') acu=.true.
         if(str(last:last).eq.'!') aft=.true.
         if(acu.and.aft) then
            write(6,*) ' Input wrong. Choose one format only.'
            go to 80
         endif
         if(acu.or.aft) last=last-1
         filnm=str(1:last)
         inquire(file=filnm,exist=fex)
         if(fex) then
            write(6,*) ' File already exists, choose another name.'
            go to 80
         endif
         open(2,file=filnm,err=92,status='new')
         mem=mem+1
      else
         prmp=' Output will be added to file:'//filnm
         write(6,*) prmp
      endif
      nf=1
      nl=mxn
      nuno=0
      prmp=
     +' Enter:(nf (nl))    (/=close)           '//
     +'        (h=help, q=skip)          *ocf*'
 90   write(6,*) prmp
      str='ocf>'
      call readin(ifc,str,last,ifh)
      if(ifc) return
      if(last.eq.0) go to 60
      if(ifh) then
         write(6,1000) 
 1000    format('        COORDINATE FILE OUTPUT'//
     +'   To output all atoms, just push return. '/
     +'   If you wish to output a sub-range nf, nl  then enter it.'/
     +' What appears on the screen is added at the end of the'/
     +' output file with the format displayed.'//
     +'   To cancel (skip) output type:q.'/
     +'   To close the currently open output file, type:/   You will'/
     +' then be prompted for a new filename.'/)
         go to 90
      endif
      if(last.eq.1.and.str(1:1).eq.'q') return
      if(last.eq.1.and.str(1:1).eq.'/') then
         close (2)
         mem=0
         write(6,*)' File closed. You can now give a new name.'
         go to 70
      endif
      call intstr(str(1:last),nuno,inum)
      if(nuno.gt.2) go to 91
      if(nuno.eq.1) then
         nf=inum(1)
      elseif(nuno.eq.2) then
         nf=inum(1)
         nl=inum(2)
      endif
      if(nf.eq.0) return
      if(fls(nx,nf).or.fls(nx,nl)) go to 91
      if(nf.gt.nl) go to 91
 60   continue
      do 6 i=nf,nl
         if(.not.ifx(i)) go to 6
         x=a(i,1)
         y=a(i,2)
         z=a(i,3)
         iq=('  ')
         if(ifz) iq=atmnme(nq(i))
         smb=('     ')
         if(ifz) smb= symb(iq,nnq(i))
*
*  output to logical unit 2 written here
*
         if(aft) then
            write(2,2) smb,x,y,z
            write(6,2) smb,x,y,z
         elseif(acu) then
            write(2,22) smb,x,y,z
            write(6,22) smb,x,y,z
         else
            write(2,222) iq,x,y,z
            write(6,222) iq,x,y,z
         endif
 2       format(tr1,a5,3f15.9)
 22      format(tr1,a5,3f24.12)
 222     format(tr1,a2,tr1,3f10.6)
*
 6    continue
      return
 91   write(6,*) ' Input wrong to coordinate output type:q to cancel.'
      go to 90
 92   write(6,*) ' File could not be opened on unit 2.'
      return
 1001 format(' The file-name for COORDINATE OUTPUT must be a new name.'/
     +' Once opened the file stays open until you close it. If you'/
     +' give the command:ocf again during the run you will only get a'/
     +' message saying which file is open and then the output prompt.'/
     +' What is then displayed is added to the end of the file.'/
     +'   The output prompt gives you the alternative to close the'/
     +' currently open output file. You will then be prompted for a'/
     +' new name.'/
     +'   The default FORMAT (needs no input) is as follows:'/
     +'          format(tr1,a2,tr1,3f10.6)'/
     +' Here a2 is a field for the chemical symbol Xx of the atom. If'/
     +' you put  format = !  you get'/
     +'          format(tr1,a5,3f15.9)'/
     +' Here a5 is a field containing chemical symbol Xx followed by'/
     +' the number No(Z) of this atom counting atoms with same Z.'/
     +' If you put  format = +  you get'/
     +'          format(tr1,a5,3f24.12)'/
     +' This high accuracy output is recommended when a high symmetry'/
     +' output is intended to be read in again (command:icf).'/
     +'   If there is no atomic number input the a-fields are blank.')
      end
*====================================================================
      subroutine osacom(nx,ifx,nq,a,noa,mxn,ifz,nnh,nmz,imz,mz)
      double precision a(nx,3),x,y,z
      logical ifx(nx),ifz,fex, open, empty, ifzr
      integer nx,nq(nx),nnh, noa,mxn,  nmz,mz(nx,3),imz(nx), nn, jmx
      integer  i,j, io, ncount, n(5), mcount, k(15), kcount, nsaves
      integer noar, nnhr, nmzr, izn(15), in(5),jn(5),kn(5)
      character line*10
*
      empty=(noa.eq.0).and.(.not.ifz).and.(nmz.eq.0)
      if(empty) then
         write(6,*) ' There is nothing to save.'
         return
      endif
      inquire(unit=3,opened=open)
      if(open) go to 80
      inquire(file='ecldsv.dat',exist=fex)
      if(fex) open (3,file='ecldsv.dat', err=91,status='old')
      if(.not.fex) then
         open (3,file='ecldsv.dat', err=91,status='new')
         nsaves=1
         go to 100
      endif
 80   continue
      rewind 3
*
*  Read file till end of file and check it:
      do 6 nsaves=0,900
*
         read(3,30,iostat=io) noar,ifzr,nnhr, nmzr
         if(io.ne.0) go to 90
*
         if(noar.gt.0) then
            do 2 i=1,noar
               read(3,40,iostat=io) nn, x, y, z
               if(io.ne.0) go to 92
 2          continue
         endif
         if(nnhr.gt.0) then
            do 3 i=1,nnhr/15
               read(3,50,iostat=io) (k(j),izn(j), j=1,15)
               if(io.ne.0) go to 92
 3          continue
            jmx=nnhr-15*(nnhr/15)
            read(3,50,iostat=io) (k(j),izn(j), j=1,jmx)
            if(io.ne.0) go to 92
         endif
         if(nmzr.gt.0) then
            do 4 i=1,nmzr/5
               read(3,60,iostat=io)(k(j),in(j),jn(j),kn(j),izn(j),j=1,5)
               if(io.ne.0) go to 92
 4          continue
            jmx=nmzr-5*(nmzr/5)
            read(3,60,iostat=io)(k(j),in(j),jn(j),kn(j),izn(j),j=1,jmx)
            if(io.ne.0) go to 92
         endif
         read(3,70,iostat=io) line
         if(io.ne.0) go to 90
*     
 6    continue
      write(6,*) ' Too many saves on file (more than 900)!'
      return
*
*               Start writing save:
*
 90   continue
      if(io.gt.0) then
         write(6,*) 
     +        ' Check of save file ended with IO status',io,' after',
     +        nsaves,' saves (or one less).'
         return
      endif    
      nsaves=nsaves+1
 100  ncount=0
      write(3,30,iostat=io) noa,ifz,nnh, nmz
      if(io.ne.0) go to 192
      do 101 i=1,mxn
         if(.not.ifx(i)) go to 101
         write(3,40,iostat=io) i,a(i,1),a(i,2),a(i,3)
         if(io.ne.0) go to 192
         ncount=ncount+1
 101  continue
      if(nnh.gt.0) then
         kcount=0
         do 102 i=1,nx
            if(nq(i).ne.1) then
               kcount=kcount+1
               k(kcount)=i
            endif
            if(kcount.eq.15) then
               kcount=0
               write(3,50,iostat=io) (k(j),nq(k(j)), j=1,15)
               if(io.ne.0) go to 192
            endif
 102     continue
         write(3,50,iostat=io) (k(j),nq(k(j)), j=1,kcount)
         if(io.ne.0) go to 192
      endif
      if(nmz.eq.0) go to 104
      mcount=0
      do 103 i=1,nx
         if(imz(i).ne.0) then
            mcount=mcount+1
            n(mcount)=i
         endif
         if(mcount.eq.5) then
            mcount=0
            write(3,60,iostat=io) 
     +      (n(j),mz(n(j),1),mz(n(j),2),mz(n(j),3),imz(n(j)),j=1,5)
            if(io.ne.0) go to 192
         endif
 103  continue
      write(3,60,iostat=io) 
     +    (n(j),mz(n(j),1),mz(n(j),2),mz(n(j),3),imz(n(j)),j=1,mcount)
      if(io.ne.0) go to 192
 104  write(3,70,iostat=io) '----------'
      if(io.ne.0) go to 192
*
      write(6,10) ncount, nsaves
 10   format(
     +tr1,i3,' atoms written to file:  ecldsv.dat  as save number',i4)
      return
 30   format(i4,l2,i4,i4)
 40   format(tr1,i3,3f25.18)
 50   format(tr3,15(i3,i2))
 60   format(5(tr1,4i3,i2))
 70   format(a10)
 91   write(6,*) ' File: ecldsv.dat  could not be opened on unit 3.'
      return
 92   continue
      if(io.eq.-1) then
         write(6,*) ' Unexpected end-of-file on file:  ecldsv.dat .'
      else
         write(6,*) ' Error in access when reading:  ecldsv.dat .'
         write(6,*) ' IO status is ', io,' nsaves',  nsaves
      endif
      return
 192  write(6,*)' Error when writing to file, IO status is ', io,
     +' number of atoms written are:', ncount
      return
      end














