*=====================================================================
* FIRST GEOMETRIC ROUTINES: GPRCOM GTRCOM GROCOM GSCCOM GRECOM GINCOM
* THEN EDIT ROUTINES: EMOCOM EORCOM EPECOM ERRCOM ERMCOM ECPCOM EZMCOM
*=====================================================================
      subroutine gprcom(str,ifc,nx,ifx,nq,a,mxn)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifh,ifc
      character str*(*)
      dimension nq(nx),a(nx,3),ivec(10)
      dimension r(3),x(3),y(3),z(3),w(3),xm(3),ym(3)
 1    write(6,*)
     +' Enter: nf ( nl ( ndf ( ndl)))   to tran',
     +'sform        (h = help)           *gpr*'
      str='gpr>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         if(last.eq.1) then
            write(6,1000)
            write(6,1100)
            go to 1
         else
            write(6,1200)
            go to 1
         endif
      endif
      call intstr(str(1:last),nun,ivec)
      if(nun.eq.0) return
      nf=ivec(1)
      nl=mxn
      if(nun.ge.2) nl=ivec(2)
      ndf=nf
      ndl=nl
      if(nun.ge.3) then
         ndf=ivec(3)
         ndl=mxn
         if(nun.eq.4) ndl=ivec(4)
      endif
      if(fls(nx,nf).or.fls(nx,nl).or.fls(nx,ndf).or.fls(nx,ndl))goto 200
      if(nun.gt.4.or.nf.gt.nl.or.ndf.gt.ndl) go to 200
      if(noar(nx,ifx,ndf,ndl).eq.0) go to 200
      qnq=0.
      do 10 i=1,3
         r(i)=0.D0
         xm(i)=0.D0
 10      ym(i)=0.D0
      do 21 i=ndf,ndl
         if(.not.ifx(i).or.nq(i).le.0) go to 21
         qnq=qnq+dble(nq(i))
         do 20 j=1,3
 20         r(j)=r(j)+dble(nq(i))*a(i,j)
 21   continue
      if(qnq.eq.0.) go to 5
      do 24 j=1,3
 24      r(j)=r(j)/qnq
      do 25 i=nf,nl
         if(.not.ifx(i)) go to 25
         do 26 j=1,3
 26         a(i,j)=a(i,j)-r(j)
 25   continue
      dist=abv(r)
      if(dist.gt.1.d-6) write(6,1001) dist
 1001 format('  Translation done;  distance:',f10.6)
      if(ndf.eq.ndl) return
      do 30 i=ndf,ndl
         if(.not.ifx(i).or.nq(i).le.0) go to 30
         do 31 j=1,3
            if(dabs(a(i,j)).lt.1.D-7) a(i,j)=0.D0
            xm(j)=xm(j)+dble(nq(i))*a(i,j)*a(i,j)
 31         ym(j)=ym(j)+dble(nq(i))*a(i,j)*a(i,mod(j,3)+1)
 30      continue
      tenor=(xm(1)+xm(2)+xm(3))/3.D0
      do 32 j=1,3
         xm(j)=xm(j)/tenor
 32      ym(j)=ym(j)/tenor
      call eigval(xm,ym,w,ndeg)
      write(6,101) ndeg
 101  format('  Degeneracy:',i3)
      if(ndeg.eq.3) go to 90
      call eigvec(w(3),xm,ym,z,ier)
      if(ier.ne.0) go to 5
      if(ndeg.eq.2) go to 80
      call eigvec(w(2),xm,ym,y,ier)
      if(ier.ne.0) go to 5
      call vpr(0,y,z,x)
      if(scp(x,x).lt.0.999999d0) go to 5
      go to 81
 90   continue
      do 65 j=ndf,ndl
         if(.not.ifx(j)) go to 65
         do 60 i=1,3
 60         z(i)=a(j,i)
         if(scp(z,z).gt.1.D-2) go to 80
 65   continue
      go to 5
 80   continue
      do 45 j=ndf,ndl
         if(.not.ifx(j)) go to 45
         do 40 i=1,3
 40         x(i)=a(j,i)
         call vpr(1,z,x,y)
         if(scp(y,y).gt.0.9999999d0) go to 41
 45   continue
      x(1)=z(1)+1.D0
      x(2)=z(2)
      x(3)=z(3)
      if(z(2)**2+z(3)**2.lt.1.D-6) x(2)=x(2)+1.D0
      call vpr(1,z,x,y)
 41   call vpr(1,y,z,x)
 81   call vpr(0,x,y,z)
      do 50 i=nf,nl
         if(.not.ifx(i)) go to 50
         w(1)=x(1)*a(i,1)+x(2)*a(i,2)+x(3)*a(i,3)
         w(2)=y(1)*a(i,1)+y(2)*a(i,2)+y(3)*a(i,3)
         w(3)=z(1)*a(i,1)+z(2)*a(i,2)+z(3)*a(i,3)
         do 53 j=1,3
 53         a(i,j)=w(j)
 50   continue
      write(6,*) ' Transformation done.'
      return
 5    write(6,*) ' Something wrong, transformation incomplete or none.'
      return
 200  write(6,*) ' Input to principal axes transformation wrong.'
      go to 1
 1000 format(
     +'        PRINCIPAL AXES TRANSFORMATION'/
     +'   This option transforms the Cartesian coordinates to a new'/
     +' coordinate system through a combined translation and rotation.'/
     +' The origin is translated to the center of charge and the axes'/
     +' are rotated to eigen-vectors of the symmetric tensor formed by'/
     +' summing: Z*(xx,yy,zz,xy,xz,yz)  Here Z are the atomic numbers'/
     +' and unless you have changed them they are =1. A molecule with'/
     +' a threefold or higher symmetry axis has a tensor with two'/
     +' equal (degenerate) eigenvalues. The eigenvectors corresponding'/
     +' to these will define the new xy-plane so the symmetry axis be-'/
     +' comes z-axis. A planar object will be rotated to the xy-plane'/
     +' since the zero eigen-value vector is choosen as z-axis.')
 1100 format(
     +'        INPUT   '/
     +'   Minimum input is: nf. The transformation then affects and is'/
     +' defined by all entered atoms from nf and up (nf,nf+1,..).'/
     +'   If you give: nf, nl  it affects and is defined by the atoms'/
     +' that have been entered in the range with first atom at nf and'/
     +' last atom at nl  (nf,nf+1,...,nl-1,nl).'/
     +'   A different DEFINING range: ndf, ndl can be given after the'/
     +' nf, nl  input (both must then be given). The first range (nf,'/
     +' nl) then represents the atoms to be affected by the transfor-'/
     +' mation and the second those that define it.        (hh = more',
     +' help)')
 1200 format('     HOW THE NEW AXES ARE DETERMINED'/
     +'   When there is no degeneracy (i.e. all eigenvalues are diffe-'/
     +' rent) the x-axis is taken to be the eigenvector corresponding'/
     +' to the largest eigenvalue, the y-axis to the middle one and'/
     +' the z-axis to the smallest eigenvalue.'/
     +'   In case of two equal eigenvalues (degeneracy 2) the eigen-'/
     +' vector corresponding to the different eigenvalue is taken as'/
     +' z-axis. The x-axis is then determined so that the lowest num-'/
     +' bered atom (in the defining range) not on the z-axis is in the'/
     +' xz-plane with positive x-coordinate.'/
     +'   Finally in case of triple degeneracy (all eigenvalues equal)'/
     +' the lowest numbered atom (in defining range) not at the origin'/
     +' is taken to define the z-axis. The x-axis is then determined'/
     +' in the same way as for degeneracy 2.'/
     +'   The coordinate system is kept right-handed. Any atoms with'/
     +' zero or negative charge are ignored in the definition of the'/
     +' transformation (except to resolve degeneracy).')
      end
*====================================================================
      subroutine gtrcom(str,ifc,nx,ifx,nq,a,mxn)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,cpy,ifh,ifc,dst
      character str*(*)
      dimension tv(3),fnum(10),a(nx,3),nq(nx),ivec(10)
 99   write(6,*)
     +' Enter: vector (# dist)     (% nf nl)   ', 
     +' (+ = copy)        (h = help)     *gtr*'
      str='gtr>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1000)
         go to 99
      endif
      cpy=.false.
      dst=.false.
      if(index(str,'+').ne.0) cpy=.true.
      idist=index(str,'#')
      nf=1
      nl=mxn
      irng=index(str,'%')
      if(irng.gt.0) then
         if(irng.lt.idist) go to 91
         call intstr(str(irng:last),nint,ivec)
         if(nint.ge.1) nf=ivec(1)
         if(nint.eq.2) nl=ivec(2)
         if(nint.gt.2.or.fls(nx,nf).or.fls(nx,nl)) go to 91
      endif
      nl=min(nl,mxn)
      if(cpy.and.mxn+nl-nf+1.gt.nx) then
         write(6,*) ' Not room enough for copy.'
         return
      endif
      if(irng.eq.0) irng=last
      if(idist.gt.0) then
         call chtflp(str(idist:irng),nun,fnum)
         if(nun.lt.1.or.nun.gt.1) go to 91
         dst=.true.
         dist=fnum(1)
      endif
      if(idist.gt.0) irng=idist
      call chtflp(str(1:irng),nun,fnum)
      if(nun.eq.1) then
         n1=int(fnum(1))
         if(fnum(1)-float(n1).ne.0.) go to 91
         if(n1.eq.0) then
            do 14 i=1,3
 14            tv(i)=0.D0
            qnq=0.
            do 11 i=1,mxn
               if(.not.ifx(i).or.nq(i).le.0) go to 11
               do 12 j=1,3
 12               tv(j)=tv(j)+a(i,j)*dble(nq(i))
               qnq=qnq+dble(nq(i))
 11            continue
            if(qnq.eq.0.) then
               write(6,*)
     +              ' Center of charge undefined. No translation done.'
               return
            endif
            do 13 i=1,3
               tv(i)=-tv(i)/qnq
 13            if(dabs(tv(i)).lt.1.D-8) tv(i)=0.D0
         elseif(n1.ne.0) then
            if(fls(nx,n1)) go to 91
            if(.not.ifx(n1)) go to 91
            do 3 i=1,3
 3             tv(i)=-a(n1,i)
         endif
      elseif(nun.eq.2) then
         n1=int(fnum(1))
         if(fnum(1)-float(n1).ne.0.) go to 91
         n2=int(fnum(2))
         if(fnum(2)-float(n2).ne.0.) go to 91
         if(fls(nx,n1).or.fls(nx,n2)) go to 91
         if(.not.ifx(n1).or..not.ifx(n2)) go to 91
         do 4 i=1,3
 4          tv(i)=a(n2,i)-a(n1,i)
      elseif(nun.eq.3) then
         tv(1)=fnum(1)
         tv(2)=fnum(2)
         tv(3)=fnum(3)
      else
         go to 91
      endif
      vdist=abv(tv)
      if(vdist.lt.1.d-6) then
         write(6,*) ' Input specifies null translation.'
         return
      endif
      if(dst) then
         do 7 i=1,3
 7          tv(i)=tv(i)*dist/vdist
      endif
      do 5 i=nf,nl
         if(.not.ifx(i)) go to 5
         inew=i
         if(cpy) inew=mxn+i-nf+1
         do 6 j=1,3
 6          a(inew,j)=a(i,j)+tv(j)
         ifx(inew)=ifx(i)
         nq(inew)=nq(i)
 5       continue
      if(nl.ge.nf) then
         write(6,100) abv(tv),tv(1),tv(2),tv(3)
 100     format(
     +' Translation done; distance:',f9.4,'   vector=(',3f10.4,')')
         if(cpy) write(6,200) mxn+1,inew
 200     format(' Translated atoms in range:',i3,' -',i3)
      endif
      return
 91   write(6,*) ' Input to translate wrong.'
      go to 99
 1000 format(
     +'  A TRANSLATION vector can be specified in one of these ways:'//
     +' 1. vector= x  y  z  the Cartesian components (to be added to'/
     +'                     the position vectors of the atoms).'/
     +' 2. vector= m        the vector from atom m to the origin (i.e.'/
     +'                     move m to origin).'/
     +' 3. vector= 0        translate center of charge (as calculated'/
     +'                     from all entered atoms) to origin. (Note:'/
     +'                     charge is 1 unless changed).'/
     +' 4. vector= m  n     vector from atom m to atom n.'//
     +'   The translation DISTANCE along the vector can be changed by'/
     +' giving: # dist    A dist < 0 gives the opposite direction.'//
     +'   The Translation affects All entered atoms unless you spe-'/
     +' cify a sub-RANGE by placing: % nf nl  in the input line'/
     +' after the other input.'//
     +'   You can make the translation act on a COPY (of the specified'/
     +' range if given, else of all atoms) by putting a  + sign in the'/
     +' input line (the position does not matter).'/)
      end
*===================================================================
      subroutine grocom(str,ifc,nx,ifx,nq,a,mxn)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,cpy,ifh,iflg,ifc
      character str*(*)
      dimension qr(4),qv(4),qw(4),v1(3),fnum(10)
      dimension nq(nx),a(nx,3),ivec(10)
      common/pp/conv
 99   write(6,*)
     +' Enter: axis  angle   or: (-)m (n)    (%',
     +' nf nl) (+ = copy)   (h = help)   *gro*'
      str='gro>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1000)
         write(6,1001)
         go to 99
      endif 
      cpy=.false.
      if(index(str,'+').ne.0) cpy=.true.
      nf=1
      nl=mxn
      irng=index(str,'%')
      if(irng.ne.0) then
         if(irng.eq.last) go to 91
         call intstr(str(irng+1:last),nint,ivec)
         if(nint.eq.0) go to 91
         if(nint.ge.1) nf=ivec(1)
         if(nint.eq.2) nl=ivec(2)
         if(nint.gt.2.or.fls(nx,nf).or.fls(nx,nl)) go to 91
         if(nf.gt.nl) go to 91
         last=irng
      endif
      if(nf.gt.nl) go to 91
      nl=min(nl,mxn)
      if(cpy.and.mxn+nl-nf+1.gt.nx) then
         write(6,*) ' Not room enough for copy.'
         return
      endif
      call chtflp(str(1:last),nun,fnum)
      if(nun.eq.0.or.nun.gt.4) go to 91
      ivpos=index(str(1:last),'!')
* 
      if(ivpos.eq.0) then
         if(nun.le.3) then
            if(nun.eq.3) then
               th=fnum(1)
               ph=fnum(2)
               alf=fnum(3)
            elseif(nun.lt.3) then
               m=int(fnum(1))
               n=0
               if(nun.eq.2) n=int(fnum(2))
               iflg=.true.
               if(m.lt.0) iflg=.false.
               m=iabs(m)
               if(fls(nx,m).or.n.gt.nx.or.n.lt.0) go to 91
               if(.not.ifx(m)) go to 91
               if(n.eq.0) then
                  do 1 i=1,3
 1                   v1(i)=a(m,i)
                  d1=abv(v1)
                  if(d1.lt.1.D-4) go to 91
                  do 2 i=1,3
 2                   v1(i)=v1(i)/d1
               else
                  if(.not.ifx(n)) go to 91
                  ierr=0
                  call bij(nx,a,m,n,v1,d1,ierr)
                  if(ierr.eq.1) go to 91
               endif
               call spheri(v1(1),v1(2),v1(3),r,ph,th,-1)
               if(iflg) then
                  alf=th
                  ph=ph-90.D0
                  th=90.D0
               else
                  th=0.D0
                  alf=-ph
                  ph=0.D0
               endif
            endif
            call spheri(x,y,z,1.D0,ph,th,1) 
         elseif(nun.eq.4) then
            v1(1)=fnum(1)
            v1(2)=fnum(2)
            v1(3)=fnum(3)
            d1=abv(v1)
            if(d1.lt.1.D-4) go to 91
            x=v1(1)/d1
            y=v1(2)/d1
            z=v1(3)/d1
            alf=fnum(4)
         endif
      elseif(ivpos.ne.0) then
         if(nun.eq.3) then
            m=int(fnum(1))
            n=int(fnum(2))
            if(fls(nx,m).or.fls(nx,n)) go to 91
            if(.not.ifx(m).or..not.ifx(n)) go to 91
            ierr=0
            call bij(nx,a,m,n,v1,d1,ierr)
            if(ierr.eq.1) go to 91
            x=v1(1)
            y=v1(2)
            z=v1(3)
            alf=fnum(3)
         elseif(nun.eq.2) then
            m=int(fnum(1))
            if(fls(nx,m)) go to 91
            if(.not.ifx(m)) go to 91
            do 3 i=1,3
 3             v1(i)=a(m,i)
            d1=abv(v1)
            if(d1.lt.1.D-4) go to 91
            x=v1(1)/d1
            y=v1(2)/d1
            z=v1(3)/d1
            alf=fnum(2)
         else
            go to 91
         endif
      endif
*
      qr(1)=dcos(alf*conv/2.D0)
      s=dsin(alf*conv/2.D0)
      qr(2)=s*x
      qr(3)=s*y
      qr(4)=s*z
      qv(1)=0.D0
      do 6 i=nf,nl
         if(.not.ifx(i)) go to 6
         do 4 j=1,3
 4          qv(j+1)=a(i,j)
         call rot(qr,qv,qw)
         inew=i
         if(cpy) inew=mxn+i-nf+1
         do 5 j=1,3
 5         a(inew,j)=qw(j+1)
         ifx(inew)=ifx(i)
         nq(inew)=nq(i)
 6    continue
      if(nl.ge.nf) then
         write(6,100) alf,x,y,z
 100     format(
     +' Rotation done; angle:',f9.4,'   axis=(x,y,z)=(',3f9.5,' )')
         if(cpy) write(6,200) mxn+1,inew
 200     format(' Rotated atoms in range:',i3,' -',i3)
         endif
      return
*
 91   write(6,*) ' Input to rotate wrong.'
      go to 99
 1000 format(
     +'  A ROTATION around an axis through the origin can be specified'
     +//' Either (I) in terms of AXIS and ANGLE',
     +' (of rotation) as follows:'/
     +' 1.  axis = tht  phi    spherical coordinates of axis-direction'/
     +' 2.  axis = x  y  z     cartesian coordinates of axis'/
     +' 3.  axis = ! m         axis to be position vector of atom m'/
     +' 4.  axis = ! m  n      axis parallel to vector from m to n'/
     +' and  angle = angle of rotation around the axis.'//
     +'  Or (II) in terms of desired RESULT of rotation as follows:'/
     +' A. Enter:m       rotate position vector of atom m to z-axis'/
     +' B. Enter:m n     make vector from m to n parallel to z-axis'/
     +' C. Enter:-m      rotate position vector to xz-plane (y=0)'/
     +' D. Enter:-m n    make vector from m to n parallel to xz-plane')
 1001 format(
     +' (Rotations in case A and B will be around an axis in the xy-'/
     +' plane. In case C and D the axis is the z-axis.)'//
     +'   To specify the RANGE, place the symbol % followed by nf,nl'/
     +' last in the input line, otherwise all atoms will be rotated.'/
     +'   If you place a + sign (anywhere) in the input line, a COPY'/
     +' of the range will be made and the rotation will act on it.'/)
      end
*====================================================================
      subroutine gsccom(str,ifc,nx,ifx,nq,a,mxn)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifh,ncp,ifc,cnv
      dimension nq(nx),a(nx,3),fnum(10),ivec(10),v1(3)
      character str*(*),prmp*79
      parameter ( bohr = 0.529177249d0 )
      prmp=' Enter:(i) factor (% nf ( nl))     (+ = copy)  (h = help)'
      prmp(75:79)='*gsc*'
 99   write(6,*) prmp
      str='gsc>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1000)
         go to 99
      endif
      cnv=.false.
      ncp=.true.
      if(index(str,'+').ne.0) ncp=.false.
      if(index(str,'#').ne.0) cnv=.true.
      nf=1
      nl=mxn
      irpos=index(str,'%')
      if(irpos.ne.0) then
         call intstr(str(irpos:last),nint,ivec)
         if(nint.gt.2) go to 90
         if(nint.ge.1) nf=ivec(1)
         if(nint.ge.2) nl=ivec(2)
         if(fls(nx,nf).or.fls(nx,nl)) go to 90
         if(nl.eq.0) return
         if(nf.gt.nl) go to 90
         nl=min(mxn,nl)
      else
         irpos=last
      endif
      call chtflp(str(1:irpos),nun,fnum)
      call intstr(str(1:irpos),nint,ivec)
      n1=0
      n2=0
      fact=0.
      if(nun.ne.nint) then
         if(nun.gt.1) go to 90
         if(cnv) go to 90
         fact=fnum(1)
      else
         if(nint.eq.0.and.(.not.cnv)) go to 90
         if(nint.eq.0.and.cnv) fact=1.d0/bohr
         if(nint.eq.1.and.ivec(1).ge.0) fact=ivec(1)
         if(nint.ge.1.and.ivec(1).lt.0) n1=-ivec(1)
         if(nint.eq.2) n2=ivec(2)
         if(nint.gt.3) go to 90
      endif
      if(n1.ne.0.or.n2.ne.0) then
         if(fls(nx,n1).or.(.not.ifx(n1))) go to 90
         if(n2.eq.0) fact=1.d0/dsqrt(a(n1,1)**2+a(n1,2)**2+a(n1,3)**2)
         if(n2.ne.0) then
            if(fls(nx,n2)) go to 90
            if(.not.ifx(n2)) go to 90
            ierr=0
            call bij(nx,a,n1,n2,v1,fact,ierr)
            if(ierr.gt.0) fact=0.
            fact=1.d0/fact
         endif
      endif
      if(fact.le.0.) go to 90
      if(str(1:1).eq.'i') fact=1./fact
      if(.not.ncp.and.mxn+nl-nf+1.gt.nx) then
         write(6,*) ' Not room enough for copy.'
         return
      endif
      nc=0
      ncount=0
      do 4 nn=nf,nl
         nc=nc+1
         if(.not.ifx(nn)) go to 4
         in=mxn+nc
         if(ncp) in=nn
         nq(in)=nq(nn)
         do 3 j=1,3
 3          a(in,j)=fact*a(nn,j)
         ifx(in)=.true.
         ncount=ncount+1
 4    continue
      if(ncp) mxn=nf-1
      if(ncount.gt.0) write(6,5) ncount,mxn+1,in,fact
 5    format(tr1,i3,tr1,'scaled atoms in range:',i3,' - ',i3,
     +'  scale factor:',f15.9)
      return
 90   write(6,*) ' Input wrong to scale.'
      go to 99
 1000 format('                 ON SCALING'/
     +'   To scale, i.e. multiply all lengths (cartesian coordinates)'/
     +' with a given numerical factor, just enter this factor. This'/
     +' can be done in one of the following ways:'//
     +' 1. factor= real number > 0.'/
     +' 2. factor= -m    distance from origin to atom m scaled to 1.'/
     +' 3. factor= -m n  distance from atom m to atom n scaled to 1.'/
     +' 4. factor= #     converts from angstroms to atomic units.'//
     +' The scaling will be performed on all entered atoms unless a'/
     +' sub-range is given by placing % nf (nl) at the end of the in-'/
     +' put line. If nl is not given all atoms with numbers nf and'/
     +' higher will be scaled.'/
     +'   A  + sign in the input line will make the scaling act on a'/
     +' COPY (of the sub-range, if specified). The copy is placed '/
     +' immediately after the entered atoms.'/
     +'   If the letter: i  is placed first in the input line the '/
     +' multiplying factor will be the INVERSE of the given number.'/)
      end
*======================================================================
      subroutine grecom(str,ifc,nx,ifx,nq,a,mxn)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifh,ncp,ifc
      integer nq(nx),ivec(10)
      double precision scp,v(3),a(nx,3),fnum(10)
      character str*(*)
 90   write(6,*) 
     +' Enter: plane   or: x  y  z (% nf (nl)) ',
     +'(/ = no copy)  (h = help)         *gre*'
      str='gre>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1000)
         go to 90
      endif
      ncp=.false.
      if(index(str,'/').ne.0) ncp=.true.
      nf=1
      nl=mxn
      irpos=index(str,'%')
      if(irpos.ne.0) then
         call intstr(str(irpos:last),nint,ivec)
         if(nint.gt.2) go to 91
         if(nint.ge.1) nf=ivec(1)
         if(nint.ge.2) nl=ivec(2)
         if(fls(nx,nf).or.fls(nx,nl)) go to 91
         if(nl.eq.0) return
         if(nf.gt.nl) go to 91
         nl=min(mxn,nl)
      endif
      if(.not.ncp.and.mxn+nl-nf+1.gt.nx) then
         write(6,*) ' Not room enough for copy.'
         return
      endif
      if(irpos.eq.0) irpos=last
      call chtflp(str(1:irpos),nun,fnum)
      if(nun.eq.0) then
         v(1)=0.
         v(2)=0.
         v(3)=0.
         if(index(str(1:irpos),'x').ne.0) then
            v(1)=1.
         elseif(index(str(1:irpos),'y').ne.0) then
            v(2)=1.
         elseif(index(str(1:irpos),'z').ne.0) then
            v(3)=1.
         else
            go to 91
         endif
      elseif(nun.eq.3) then
         do 1 i=1,3
 1          v(i)=fnum(i)
      else
         go to 91
      endif
      dd=scp(v,v)
      if(dd.lt.1.D-10) go to 91
      ncount=0
      do 2 i=nf,nl
         ncount=ncount+1
         if(.not.ifx(i)) go to 2
         sc=0.D0
         do 3 j=1,3
 3          sc=sc+a(i,j)*v(j)
         ii=mxn+ncount
         if(ncp) ii=i
         do 4 j=1,3
 4          a(ii,j)=a(i,j)-2.D0*v(j)*sc/dd
         ifx(ii)=.true.
         nq(ii)=nq(i)
 2    continue
      if(ncount.eq.0) return
      if(ncp) then
         write(6,*) ' Reflection done.'
      else
         write(6,100) mxn+1,ii
      endif
 100  format(' Reflected atoms in range:',i3,' -',i3)
      return
 91   write(6,*) ' Input to reflect wrong.'
      go to 90
 1000    format(
     +'   To copy and REFLECT all atoms or a subrange nf,nl in a plane'/
     +' through the origin you can use one of the following options:'/
     +'   (1) To get a reflection in one of the coordinate-planes,'/
     +' (x=0, y=0, or z=0), just type:x  or y  or  z    for the corre-'/
     +' sponding plane.'/
     +'   (2) To get a reflection in the plane through the origin per-'/
     +' pendicular to the vector with components x,y,z, just enter'/
     +' these components. Example: input of the line:  1.0  1.0  0.0'/
     +' results in a reflection in a plane through the z-axis making'/
     +' the angle 45 degrees with both the planes x=0 (the yz-plane)'/
     +' and the plane y=0 (the xz-plane).'/
     +'   All entered atoms will automatically be copied and reflected'/
     +' unless you specify a sub-RANGE by putting % nf nl  last in the'/
     +' input line. (Example: % 3 8).'/
     +'   If you put a: /  in the input line you avoid the copying.')
      end
*===================================================================
      subroutine gincom(str,ifc,nx,ifx,nq,a,mxn)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifh,ncp,ifc
      dimension nq(nx),a(nx,3),ivec(10)
      character str*(*),prmp*79
      prmp=' Enter: nf  (nl)    (/ = no copy)  to invert    (h = help)'
      prmp(75:79)='*gin*'
 99   write(6,*) prmp
      str='gin>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1000)
 1000    format(
     +'     ON  INVERSION'//
     +'   This option allows you to invert the atoms in the specified'/
     +' range nf, nl  through the origin i.e. the position vectors of'/
     +' the atoms change sign.'/
     +'   Minimum input is  nf .  All entered atoms with numbers nf'/
     +' and higher are then copied and inverted. The copy is placed'/
     +' after the highest numbered entered atom.'/
     +'   If you do not want a copy, place a  /  in the input line.'/)
         go to 99
      endif
      ncp=.false.
      if(index(str,'/').ne.0) ncp=.true.
      nf=1
      nl=mxn
      call intstr(str(1:last),nun,ivec)
      if(nun.lt.1.or.nun.gt.2) go to 90
      nf=ivec(1)
      if(nun.eq.2) nl=ivec(2)
      if(fls(nx,nf).or.fls(nx,nl)) go to 90
      if(nf.gt.nl) go to 90
      nl=min(nl,mxn)
      if(.not.ncp.and.mxn+nl-nf+1.gt.nx) then
         write(6,*) ' Not room enough for copy.'
         return
      endif
      nc=0
      ncount=0
      do 4 nn=nf,nl
         nc=nc+1
         if(.not.ifx(nn)) go to 4
         in=mxn+nc
         if(ncp) in=nn
         nq(in)=nq(nn)
         do 3 j=1,3
 3          a(in,j)=-1.d0*a(nn,j)
         ifx(in)=.true.
         ncount=ncount+1
 4    continue
      if(ncp) mxn=nf-1
      if(ncount.gt.0) write(6,5) ncount,mxn+1,in
 5    format(tr1,i3,tr1,'inverted atoms in range:',i3,' - ',i3)
      return
 90   write(6,*) ' Input wrong to invert.'
      go to 99
      end
*===================================================================
* BELOW EDIT ROUTINES 
*===================================================================
      subroutine emocom(str,ifc,nx,ifx,nq,a,mxn)
      logical ifx(nx),ifxst,fls,ifh,ifc
      character str*(*)
      integer nq(nx),nqst,ivec(10)
      double precision a(nx,3),ast(3)
* This routine moves the range nf,nl to the position after nfnew directly
* if possible, otherwise by one step of store, shift, retrieve at the time.
 90   write(6,*) 
     +' Enter: nnew  n   or: nfnew  nf  nl    t',
     +'o move atom or range  (h = help)  *emo*'
      str='emo>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      call intstr(str(1:last),nint,ivec)
      if(ifh) then
         write(6,1010)
         write(6,1011)
         go to 90
      endif
      if(nint.lt.2.or.nint.gt.3) go to 91
      n=ivec(1)-1
      nf=ivec(2)
      nl=nf
      if(nint.eq.3) nl=ivec(3)
      if(fls(nx,ivec(1)).or.fls(nx,nf).or.fls(nx,nl).or.nl.lt.nf)
     +go to 91
      nu=nl-nf+1
      if(n+nu.gt.nx) then
         write(6,*) ' Not room enough for requested move!'
         go to 90
      endif
      nstep=n+1-nf
      if(nstep.eq.0) then
         write(6,*) ' Input implies no move!'
         go to 90
      endif
* Find out if new range is empty
      nfn=n+1
      nln=n+nu
      noan=noar(nx,ifx,nfn,nln)
      if(noan.eq.0.and.(nfn.gt.nl.or.nln.lt.nf)) then
         noam=noar(nx,ifx,nf,nl)
         if(noam.eq.0) then
            write(6,*) ' Input has no effect.'
            return
         endif
         do 11 i=1,nu
            ifx(n+i)=ifx(nf-1+i)
            ifx(nf-1+i)=.false.
            nq(n+i)=nq(nf-1+i)
            nq(nf-1+i)=1
            a(n+i,1)=a(nf-1+i,1)
            a(n+i,2)=a(nf-1+i,2)
            a(n+i,3)=a(nf-1+i,3)
 11      continue
         if(nf.lt.nl) write(6,1000) nf,nl,noam,nfn,nln
         if(nf.eq.nl) write(6,1002) nf,noam,nfn
         return
      endif
      nfn=nf
      nln=nl
      do 1 i=1,iabs(nstep)
         if(nstep.lt.0) then
            ifxst=ifx(nfn-1)
            if(ifxst) then
               nqst=nq(nfn-1)
               ast(1)=a(nfn-1,1)
               ast(2)=a(nfn-1,2)
               ast(3)=a(nfn-1,3)
            endif
            do 2 j=1,nu
               nnext=nfn-1+(j-1)
               ifx(nnext)=ifx(nnext+1)
               if(ifx(nnext)) then
                  nq(nnext)=nq(nnext+1)
                  a(nnext,1)=a(nnext+1,1)
                  a(nnext,2)=a(nnext+1,2)
                  a(nnext,3)=a(nnext+1,3)
               endif
 2             continue
            ifx(nln)=ifxst
            if(ifx(nln)) then
               nq(nln)=nqst
               a(nln,1)=ast(1)
               a(nln,2)=ast(2)
               a(nln,3)=ast(3)
            endif
            nfn=nfn-1
            nln=nln-1
         elseif(nstep.gt.0) then
            ifxst=ifx(nln+1)
            if(ifxst) then
               nqst=nq(nln+1)
               ast(1)=a(nln+1,1)
               ast(2)=a(nln+1,2)
               ast(3)=a(nln+1,3)
            endif
            do 3 j=1,nu
               nnext=nln+1-(j-1)
               ifx(nnext)=ifx(nnext-1)
               if(ifx(nnext)) then
                  nq(nnext)=nq(nnext-1)
                  a(nnext,1)=a(nnext-1,1)
                  a(nnext,2)=a(nnext-1,2)
                  a(nnext,3)=a(nnext-1,3)
               endif
 3             continue
            ifx(nfn)=ifxst
            if(ifx(nfn)) then
               nq(nfn)=nqst
               a(nfn,1)=ast(1)
               a(nfn,2)=ast(2)
               a(nfn,3)=ast(3)
            endif
            nfn=nfn+1
            nln=nln+1
         endif
 1       continue
      if(nstep.lt.0) then
         nfold=n+1
         nlold=nf-1
         nfnew=nfold+nu
         nlnew=nl
      elseif(nstep.gt.0) then
         nfold=nl+1
         nlold=n+nu
         nfnew=nf
         nlnew=n
      endif
      noam=noar(nx,ifx,nfn,nln)
      noan=noar(nx,ifx,nfnew,nlnew)
      if(nf.lt.nl) write(6,1000) nf,nl,noam,nfn,nln
      if(nf.eq.nl) write(6,1002) nf,noam,nfn
      if(noan.gt.0) write(6,1001) nfold,nlold,nfnew,nlnew,noan
 1000 format(
     +'       The input range:',i3,' -',i3,'   containing:',i3,' atoms'/
     +' is now moved to range:',i3,' -',i3,'.')
 1002 format(
     +'     Input atom number:',i3,'        with contents:',i3,' atom'/
     +'  is now renumbered to:',i3)
 1001 format(
     +' As a result the range:',i3,' -',i3,'  has become',/
     +' renumbered and is now:',i3,' -',i3,'  It contains:',i3,' atoms')
      return
 91   write(6,*) ' Input to move wrong or inconsistent.'
      go to 90
 1010 format(
     +'   Use this option to MOVE a range:  nf - nl (a single atom n'/
     +' is considered as a range of length one with n=nf=nl) to a new'/
     +' starting position: nfnew  (or new position: nnew  for atom).'/
     +'   If the move does not affect other atoms it is done in the'/
     +' obvious way. Otherwise the move is done as follows:'/
     +'  (1) The direction of motion is determined; if nfnew is less'/
     +' than nf it is backwards, if greater forwards.'/
     +'  (2) The contents of the first atom outside the range in the'/
     +' direction of motion is stored.')
 1011 format(
     +'  (3) All atoms in the range are shifted one number in the di-'/
     +' rection of motion, starting with the one nearest the space'/
     +' freed by the storing.'/
     +'  (4) When this has been done there is a new empty space left'/
     +' behind by the shifted range. The stored atom is retrieved to'/
     +' this space.'/
     +'   This sequence of store - shift - retrieve steps is repeated'/
     +' until the desired position for the range is reached.'/
     +'   Note that atoms between the two range positions become re-'/
     +' numbered but no atoms are lost. There is a message if atoms'/
     +' outside the input range are renumbered.')
      end
*======================================================================
      subroutine eorcom(str,ifc,nx,ifx,nq,a,mxn,ifz)
      logical ifx(nx),ifxst,fls,ifh,ifz,ifc
      double precision a(nx,3),ast(3)
      integer nq(nx),nqst,ivec(10)
      character str*(*),sst*10
 90   write(6,*) 
     +' Enter: ct   or: an   or: rv   (% nf nl)',
     +'  to order   (h = help)           *eor*'
      str='eor>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1000)
         go to 90
      endif
      nf=1
      nl=mxn
      irng=index(str,'%')
      if(irng.ne.0) then
         call intstr(str(irng:last),nint,ivec)
         if(nint.ge.1) nf=ivec(1)
         if(nint.ge.2) nl=ivec(2)
         if(nint.gt.2.or.fls(nx,nf).or.fls(nx,nl)) go to 91
         if(nf.ge.nl) go to 91
      endif
      nl=min(nl,mxn)
      if(irng.ne.0) irng=irng-1
      if(irng.eq.0) irng=last
      sst=str(1:irng)
      if(index(sst,'x').ne.0) then
         call order(nx,ifx,nq,a,nf,nl,1,ifz)
      elseif(index(sst,'c').ne.0) then
         call order(nx,ifx,nq,a,nf,nl,2,ifz)
      elseif(index(sst,'s').ne.0) then
         call order(nx,ifx,nq,a,nf,nl,3,ifz)
      elseif(index(sst,'an').ne.0) then
         if(.not.ifz) go to 91
         do 2 i=nf,nl-1
            if(.not.ifx(i)) go to 2
            do 1 j=i+1,nl
               if(.not.ifx(j)) go to 1
               if(nq(j).lt.nq(i)) then
                  nqst=nq(i)
                  do 10 k=1,3
 10                  ast(k)=a(i,k)
                  nq(i)=nq(j)
                  do 11 k=1,3
 11                  a(i,k)=a(j,k)
                  nq(j)=nqst
                  do 12 k=1,3
 12                  a(j,k)=ast(k)
               endif
 1          continue
 2       continue
      elseif(index(sst,'rv').ne.0) then
         do 3 i=nf,nf-1+(nl-(nf-1))/2
            j=nl+nf-i
            ifxst=ifx(i)
            nqst=nq(i)
            do 30 k=1,3
 30            ast(k)=a(i,k)
            ifx(i)=ifx(j)
            nq(i)=nq(j)
            do 31 k=1,3
 31            a(i,k)=a(j,k)
            ifx(j)=ifxst
            nq(j)=nqst
            do 32 k=1,3
 32            a(j,k)=ast(k)
 3       continue
      else
         go to 91
      endif
      write(6,1001) nf,nl,sst
 1001 format(' Range:',i3,' -',i3,' ordered after type:',a10)
      return
 91   write(6,*) ' Input to order wrong or meaningless.'
      go to 90
 1000 format(
     +'       ORDER renumbers the atoms in some standard order.'/
     +'  (1) Input of ct (=x, c, or s) gives an order based on the'/
     +' choosen coordinate-type (ct=x: 1st after x, then after y, then'/
     +' z; ct=c: 1st after -z, then after rh, then phi; ct=s: 1st af-'/
     +' ter r, then after theta, then phi).'/
     +'  (2) Input of the letters: an =(atomic number) will order the'/
     +' atoms in the range after increasing Z-value.'/
     +'  (3) If you desire the reverse order type the letters: rv in '/
     +' the input line and the order will be reversed.'/
     +'   All atoms entered will be affected by the ordering unless'/
     +' the RANGE nf,nl is given by putting % nf nl  last in the in-'/
     +' put line.'/
     +'   Empty positions are unaffected by (1) and (2) but the'/
     +' reversing is done with a reflection in the middle of the range'/
     +' so that empty places may get occupied and vice versa.')
      end
*=====================================================================
      subroutine epecom(str,ifc,nx,ifx,nq,a,mxn)
      logical ifx(nx),ifxst,fls,ifh,ifc
      double precision a(nx,3),ast(3)
      integer nx,nq(nx),nqst,nint,ivec(10),nm,nb,n1,nl,last
      character str*(*)
 90   write(6,*) 
     +' Enter: n1  n2 ( n3..( n10)..))    to pe',
     +'rmute        (h = help)           *epe*'
      str='epe>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1001)
         go to 90
      endif
      call intstr(str(1:last),nint,ivec)
      if(nint.lt.2.or.nint.gt.10) go to 91
      do 1 i=1,nint
         if(fls(nx,ivec(i))) go to 91
 1    continue
      nl=ivec(nint)
      ifxst=ifx(nl)
      nqst=nq(nl)
      do 20 j=1,3
 20      ast(j)=a(nl,j)
      do 2 i=nint-1,1,-1
         nm=ivec(i+1)
         nb=ivec(i)
         ifx(nm)=ifx(nb)
         nq(nm)=nq(nb)
         do 21 j=1,3
 21         a(nm,j)=a(nb,j)
         if(ifx(nb).and.nb.ne.nm) write(6,1000) nb,nm
 1000    format(' Atom ',i3,'  renumbered to ',i3)
 2    continue
      n1=ivec(1)
      ifx(n1)=ifxst
      nq(n1)=nqst
      do 22 j=1,3
 22      a(n1,j)=ast(j)
      if(ifxst.and.nl.ne.n1) write(6,1000) nl,n1
      if(nint.eq.2) go to 90
      return
 91   write(6,*) ' Input to permute wrong!'
      go to 90
 1001 format(
     +'   Input to PERMUTE is a string of atom numbers: n1 n2 .. nn '/
     +' The maximum possible number of numbers is 10. '/
     +'   The result of this input is a permutation as follows:'/
     +'          atom number n1 is renumbered to n2'/
     +'          atom number n2 is renumbered to n3'/
     +'          ..................................'/
     +'          atom number nn is renumbered to n1'/
     +' I.e. a cyclic permutation. It is not necessary that the cor-'/
     +' responding atoms have been entered. If any entered atoms'/
     +' have become renumbered as a result of the permutation this is'/
     +' explicitly reported. Note that the permutation of an entered'/
     +' atom to a free number is equivalent to a move that does not'/
     +' affect the numberings of any other atoms.'/
     +'   Charges (atomic numbers), if entered, are automatically re-'/
     +' numbered together with their atoms. Z-matrix information is,'/
     +' however, not affected and will usually become obsolete after'/
     +' renumbering.'/)
      end
*=====================================================================
      subroutine errcom(str,ifc,nx,ifx,nq,a,mxn)
      logical ifx(nx),fls,ifh,ifc
      integer nq(nx),ivec(10)
      double precision a(nx,3)
      character str*(*)
 90   continue
      write(6,*)
     +' Enter: n   or: nf  nl    to remove & re',
     +'number          (h = help)        *err*'
      str='err>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1001)
         go to 90
      endif
      call intstr(str(1:last),nint,ivec)
      if(nint.eq.0.or.nint.gt.2) go to 91
      nf=ivec(1)
      nl=nf
      if(nint.eq.2) nl=ivec(2)
      if(fls(nx,nf).or.fls(nx,nl)) go to 91
      if(nf.gt.nl) go to 91
      natmd=0
      do 1 i=nf,nl
         if(ifx(i)) natmd=natmd+1
         ifx(i)=.false.
         nq(i)=1
 1    continue
      if(natmd.eq.0) write(6,*) '  Zero atoms removed.'
      if(natmd.eq.1) write(6,*) '  One atom removed.'
      if(natmd.gt.1) write(6,1000) natmd
 1000 format(1x,i5,' atoms removed.')
      natmr=0
      ncount=0
      do 2 i=nl+1,mxn
         if(.not.ifx(i)) go to 4
         natmr=natmr+1
         ifx(nf+ncount)=ifx(i)
         ifx(i)=.false.
         nq(nf+ncount)=nq(i)
         nq(i)=1
         do 3 j=1,3
 3          a(nf+ncount,j)=a(i,j)
 4       ncount=ncount+1
 2    continue
      if(natmr.eq.0) write(6,*) '  Zero atoms renumbered.'
      if(natmr.gt.0) write(6,1020) nl+1,mxn,nf,nf+ncount-1,natmr
 1020 format(' Range ',i3,' -',i3,'  has become renumbered and is',/
     +     '   now ',i3,' -',i3,'  It contains ',i3,'  atoms.')
      nln=min(mxn,nl)
      mxn=mxn-max(0,(nln-nf+1))
      if(nint.eq.1) go to 90
      return
 91   write(6,*)' Input to remove&renumber wrong.'
      go to 90
 1001 format(/
     +'   The atom n or the range nf,nl is REMOVED and all entered'/
     +' atoms with numbers larger than n or nl respectively are then'/
     +' RENUMBERED to numbers 1 or nl-nf+1(=length of range) less '/
     +' respectively. This happens irrespectively of whether the re-'/
     +' moved atoms actually existed (were entered) or not.'/
     +'   Charges (atomic numbers), if entered, are automatically re-'/
     +' numbered together with their atoms. Z-matrix information is,'/
     +' however, not affected and will usually become obsolete after'/
     +' renumbering.'/
     +'   If this option is used repeatedly it must be remembered'/
     +' that one should remove from higher numbers downwards since'/
     +' otherwise the renumberings make it hard to knows which atoms'/
     +' that actually are removed.'/)
      end
*========================================================================
      subroutine ermcom(str,ifc,nx,ifx,nq,mxn,ifz,imz)
      logical ifx(nx),fls,ifh,ifc,ifz
      integer nq(nx),imz(nx),ivec(10)
      character str*(*)
 90   continue
      write(6,*)
     +' Enter: n    or: nf nl    or: xy      to',
     +' remove       (h = help)          *erm*'
      str='erm>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1002)
         go to 90
      endif
      call intstr(str(1:last),nint,ivec)
      if(last.eq.2.and.str(1:2).eq.'an') then
         if(ifz) write(6,*) ' Atomic numbers, Z (=charges), removed.'
         ifz=.false.
         do 2 i=1,nx
 2          nq(i)=1
         return
      elseif(last.eq.2.and.str(1:2).eq.'zm') then
         nzm=0
         do 3 i=1,nx
            nzm=nzm+iabs(imz(i))
 3          imz(i)=0
         if(nzm.ne.0) write(6,*) ' Stored z-matrix information removed.'
         return
      elseif(last.eq.2.and.str(1:2).eq.'rm') then
         nf=1
         nl=nx
         go to 10
      elseif(last.eq.2.and.str(1:3).eq.'et') then
         if(ifz) write(6,*) ' Atomic numbers, Z (=charges), removed.'
         ifz=.false.
         do 20 i=1,nx
 20         nq(i)=1
         nzm=0
         do 30 i=1,nx
            nzm=nzm+iabs(imz(i))
 30         imz(i)=0
         if(nzm.ne.0) write(6,*) ' Stored z-matrix information removed.'
         nf=1
         nl=nx
         go to 10
      endif
      if(nint.eq.0.or.nint.gt.2) go to 91
      nf=ivec(1)
      nl=nf
      if(nint.eq.2) nl=ivec(2)
      if(fls(nx,nf).or.fls(nx,nl)) go to 91
      if(nf.gt.nl) go to 91
 10   natmd=0
      do 1 i=nf,nl
         if(ifx(i)) natmd=natmd+1
         ifx(i)=.false.
 1    continue
      if(natmd.eq.0) write(6,*) '  Zero atoms removed.'
      if(natmd.eq.1) write(6,*) '  One atom removed.'
      if(natmd.gt.1) write(6,1000) natmd
 1000 format(1x,i5,' atoms removed.')
      if(nint.eq.1) go to 90
      return
 91   write(6,*)' Input to remove wrong.'
      go to 90
 1002 format('          REMOVE  OPTIONS  '//
     +' ATOMS: n   or   nf  nl'/
     +'   The atom  n or the atoms in the range  nf,nl are removed,'/
     +' i.e. they are flagged as not having been entered. Their'/
     +' atomic numbers Z are unchanged (i.e. Z=1 unless changed by'/
     +' you.) Stored z-matrix information also remains unchanged.'/
     +' ATOMS: xy = rm'/
     +'   All atoms are removed; same as nf=1, nl=max n.'//
     +' ATOMIC NUMBERS (charges): xy = an'/
     +'   Atomic numbers are marked as off and reset to Z = 1.'/
     +' They will then not be displayed by  dco, dzm, or dge.'//
     +' Z-MATRIX INFORMATION: xy = zm'/
     +'   All z-matrix information that has been read from file or'/
     +' entered with  ezm  is removed. Display and output z-matrices'/
     +' will then be calculated using nearest neighbours.'//
     +' EVERYTHING: xy = et'/
     +' Same as rm, an, and, zm  given in succesion. Use to restart.'/)
      end
*========================================================================
      subroutine ecpcom(str,ifc,nx,ifx,nq,a,mxn)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,ifh,ifc,plc
      dimension nq(nx),a(nx,3),ivec(10)
      character str*(*),prmp*79
      prmp=' Enter: n   or: nf nl   (# nfc)   to copy atom or range'//
     +'  (h = help)'
      prmp(75:79)='*ecp*'
 99   write(6,*) prmp
      str='ecp>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1000)
         go to 99
      endif
      inpos=index(str,'#')
      plc=.false.
      if(inpos.ne.0) then
         call intstr(str(inpos:last),nun,ivec)
         if(nun.ne.1) go to 90
         nfc=ivec(1)
         if(fls(nx,nfc)) go to 90
         if(ifx(nfc)) go to 90
         plc=.true.
      endif
      if(inpos.eq.0) inpos=last
      call intstr(str(1:inpos),nun,ivec)
      if(nun.lt.1.or.nun.gt.2) go to 90
      nf=ivec(1)
      nl=nf
      if(nun.eq.2) nl=ivec(2)
      if(fls(nx,nf).or.fls(nx,nl)) go to 90
      if(nun.eq.2) nl=min(nl,mxn)
      if(nf.gt.nl) go to 90
      if(noar(nx,ifx,nf,nl).eq.0) then
         write(6,*) ' You can not copy an empty range.'
         return
      endif
      if(plc.and.noar(nx,ifx,nfc,nfc+nl-nf).ne.0) go to 91
      if(mxn+nl-nf+1.gt.nx) go to 91
      nc=0
      ncount=0
      do 4 nn=nf,nl
         nc=nc+1
         if(.not.ifx(nn)) go to 4
         in=mxn+nc
         if(plc) in=nfc-1+nc
         nq(in)=nq(nn)
         do 3 j=1,3
 3          a(in,j)=a(nn,j)
         ifx(in)=.true.
         ncount=ncount+1
 4    continue
      if(nun.eq.2) then
         if(.not.plc) then
            if(ncount.gt.0) write(6,5) ncount,mxn+1,in
 5          format(tr1,i3,tr1,'copied atoms in range:',i3,' - ',i3)
         else
            if(ncount.gt.0) write(6,5) ncount,nfc,nfc+nl-nf
         endif
      else
         if(.not.plc) then
            if(ncount.gt.0) write(6,6) ncount,mxn+1
 6          format(tr1,i3,tr1,'copied atom at:',i3)
         else
            if(ncount.gt.0) write(6,6) ncount,nfc
         endif
      endif
      return
 90   write(6,*) ' Input wrong to copy.'
      go to 99
 91   write(6,*) ' Not room enough for copy.'
      go to 99
 1000 format(
     +'      ON COPYING'//
     +'   This option allows you to copy an atom, n, or a given range:'/
     +' nf, nl. The copy is placed after the highest numbered en-'/
     +' tered atom. If you wish it somewhere else enter: # nfc  to '/
     +' have the copy placed in range starting with nf=nfc. The place'/
     +' thus specified for the copy must be free.'/
     +'   Atomic numbers (charges), if entered, are also copied to-'/
     +' gether with the atoms (z-matrix information is not affected).'/)
      end
*========================================================================
      subroutine ezmcom(str,ifc,nx,ifx,nq,a,mxn,nmz,imz,mz)
      logical ifx(nx),ifh,fls,ifc,cal
      integer nq(nx),imz(nx),mz(nx,3),ivec(10)
      double precision a(nx,3)
      character str*(*), alt*1
 90   continue
      write(6,*)
     +' Enter: n i j k (a)    or: x (n)(#)   to',
     +' edit z-matrix   (h = help)       *ezm*'
      str='ezm>'
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         write(6,1000)
         go to 90
      endif
      cal=.false.
      if(index(str,'#').ne.0) cal=.true.
      if(last.eq.1.and.str(1:1).eq.'r') then
         if(nmz.gt.0) write(6,*) ' All z-matrix information removed'
         do 2 i=1,nx
            imz(i)=0
 2       continue
         nmz=0
         go to 90
      elseif(last.eq.1.and.str(1:1).eq.'d') then
         if(nmz.eq.0) then
            write(6,*) ' No z-matrix information is stored'
            go to 90
         endif
         write(6,1010)
         do 3 i=3,nx
            if(imz(i).ne.0) then
               alt=' '
               if(imz(i).eq.-1) alt='a'
               if(i.eq.3) write(6,1001) i, mz(i,1),mz(i,2)
               if(i.gt.3) write(6,1002) i, mz(i,1),mz(i,2),mz(i,3), alt
            endif
 3       continue
 1001    format(tr1,i3,tr3,i3,tr2,i3)
 1002    format(tr1,i3,tr3,i3,tr2,i3,tr2,i3,tr3,a1) 
 1010    format( '   n,    i,   j,   k,  alt ')
         go to 90
      endif
      call intstr(str(1:last),nint,ivec)
      if(nint.eq.0.and.(.not.cal)) go to 91
      if(cal) then
         if(str(1:2).ne.'d ') go to 91
         call trefst(nx,ifx,n1,n2,n3)
         if(n3.eq.0) go to 90
      endif
      if(cal.and.nint.eq.0) then
         write(6,1010)
         do 12 n=n3,mxn
            if(ifx(n)) then
               call refnuc(nx,ifx,nq,a,n3,n,i,j,k)
               if(n.eq.n3) write(6,1001) n,i,j
               if(n.gt.n3) write(6,1002) n,i,j,k
            endif
 12      continue
         go to 90
      elseif(cal.and.nint.eq.1)then
         n=ivec(1)
         if(fls(nx,n)) go to 91
         if(n.lt.n3.or.(.not.ifx(n))) go to 90
         write(6,1010)
         call refnuc(nx,ifx,nq,a,n3,n,i,j,k)
         if(n.eq.n3) write(6,1001) n,i,j
         if(n.gt.n3) write(6,1002) n,i,j,k
         go to 90
      elseif(cal.and.nint.gt.1) then
         go to 91
      endif
      if(nint.eq.1) then
         i=ivec(1)
         if(fls(nx,i)) go to 91
         if(str(1:2).eq.'r ') then
            if(imz(i).ne.0) then
               write(6,1005) i
               imz(i)=0
               nmz=nmz-1
            endif
            go to 90
         elseif(str(1:2).eq.'d ') then
            if(nmz.eq.0) then
               write(6,*) ' No z-matrix information is stored'
               go to 90
            endif
            write(6,1010)
            if(imz(i).ne.0) then
               alt=' '
               if(imz(i).eq.-1) alt='a'
               if(i.eq.3) write(6,1001) i, mz(i,1),mz(i,2)
               if(i.gt.3) write(6,1002) i, mz(i,1),mz(i,2),mz(i,3), alt
            endif
            go to 90
         else
            go to 91
         endif
      endif
      if(nint.ge.3) then
         i=ivec(1)
         n1=ivec(2)
         n2=ivec(3)
         if(fls(nx,i).or.fls(nx,n1).or.fls(nx,n2)) go to 91
         if(n1.ge.i.or.n2.ge.i.or.n1.eq.n2) go to 91
         if(nint.eq.3.and.i.ne.3) go to 91
         mz(i,1)=n1
         mz(i,2)=n2
         if(nint.eq.3) then
            if(imz(i).eq.0) then
               write(6,*) 
     +              ' Z-matrix reference atoms for atom   3  stored'
               nmz=nmz+1
               imz(i)=1
            else
               write(6,*) 
     +              ' Z-matrix reference atoms for atom   3  changed'
            endif
            go to 90
         else
            if(nint.gt.4) go to 91
            if(i.lt.4) go to 91
            n3=ivec(4)
            if(fls(nx,n3).or.n3.ge.i.or.n3.eq.n1.or.n3.eq.n2) go to 90
            mz(i,3)=n3
            if(imz(i).eq.0) then
               write(6,1003) i
               nmz=nmz+1
               imz(i)=1
            else
               write(6,1004) i
            endif
            imz(i)=1
            if(str(last-1:last).eq.' a'.and.i.gt.3) imz(i)=-1
 1003       format(' Z-matrix reference atoms for atom ',i3,' stored')
 1004       format(' Z-matrix reference atoms for atom ',i3,' changed')
 1005       format(' Z-matrix information for atom ',i3,' removed')
            go to 90
         endif
      endif
 91   write(6,*)' Input to ezm (edit z-matrix) wrong.'
      go to 90
 1000 format('                  Z-MATRIX EDIT '/
     +'   Information about the reference atoms (the connectivity) to'/
     +' be used in a display (dzm) or output (ozm) z-matrix can'/
     +' be entered by typing INPUT to this command on the form:'/
     +'   n  i  j  k       where n is the atom number and i, j, and k'/
     +' are the atom numbers of the reference atoms to be used to'/
     +' specify the position of atom number n. For atom number 3 it'/
     +' is enough to give:  n  i  j.  To get the alternative z-matrix'/
     +' specification, place the letter  a  last in the input line.'//
     +'   You can DISPLAY the stored information by typing:  d '/
     +' by typing:  d  n    you get the information for atom n.'/
     +' If you put the character # in the input line you get display'/
     +' of the (nearest neighbour) calculated reference atoms.'//
     +'   Z-matrix information may already exist in the program if it'/
     +' has been entered with this command or if has been read from a'/
     +' file. It can then be REMOVED by typing: r. By typing:  r  n  '/
     +' you can remove the information for atom n. When a z-matrix is'/
     +' displayed or written to file the program will use the stored'/
     +' information; if no information is stored the reference atoms'/
     +' will be calculated according to a nearest neighbour principle.')
      end
*======================================================================
* NUMBERS FROM CHARACTER-STRINGS ROUTINES:
* CHTFLP for double precision reals.
* INTSTR for integers.
*----------------------------------------------------------------------
      subroutine chtflp(str,nun,fnum)
* The nun numbers in str are placed in the nun first elements of fnum.
* Exponential notation cannot be read by chtflp. Everything in str 
* except digits, decimal-points and minus-signs is ignored.
      integer bf(10),af(10),nun,nc,ib,i,j,k,l,iv(20),lt,np,itg,len
      double precision fnum(10)
      character str*(*),pst*70,dig*1 
      common/dgt/dig(10)
      nun=0
      lt=len(str)
      np=1
      do 5 i=1,lt
         pst=str(np:lt)
         k=0
         do 2 j=1,20
            itg=-1
            do 1 l=1,10
               if(pst(j:j).eq.dig(l)) itg=l-1
 1          continue
            if(itg.eq.-1) go to 3
            k=k+1
            iv(k)=itg
 2       continue
         write(6,*) ' WARNING from chtflp, integer string too long.'
 3       continue
         if(k.ne.0) then
            nun=nun+1
            fnum(nun)=0
            do 4 j=1,k
               fnum(nun)=fnum(nun)*10+dble(iv(j))
 4          continue
            bf(nun)=np-1
            af(nun)=np+k
            if(nun.eq.10) go to 6
            np=np+k
         endif
         np=np+1
         if(np.gt.lt) go to 6
 5    continue
 6    continue
      do 7 i=1,nun
         ib=bf(i)
         if(ib.eq.0) go to 7
         if(str(ib:ib).eq.'.') then
            fnum(i)=fnum(i)/10.D0**(af(i)-ib-1)
         endif
 7    continue
         nc=0
      do 8 i=2,nun
         ib=bf(i)
         if(af(i-1).eq.ib.and.str(ib:ib).eq.'.') then
            fnum(i-1)=fnum(i-1)+fnum(i)
            nc=nc+1
            af(nc)=i
         endif
 8    continue
      do 9 i=1,nun
         ib=bf(i)
         if(ib.eq.0) go to 9
         if(str(ib:ib).eq.'-') fnum(i)=-fnum(i)
         if(str(ib:ib).eq.'.') then
            ib=ib-1
            if(ib.eq.0) go to 9
            if(str(ib:ib).eq.'-') fnum(i)=-fnum(i)
         endif
 9    continue
      nun=nun-nc
      j=1
      do 12 i=2,nun
         do 11 k=1,nc
            if(af(k).eq.i) then
               j=j+1
               do 10 l=1,nc
                  af(l)=af(l)-1
 10             continue
            endif
 11      continue
         j=j+1
         fnum(i)=fnum(j)
 12   continue
      return
      end
*----------------------------------------------------------------------
      subroutine intstr(str,nin,ivec)
* This routine is similar to chtflp above but returns nin signed inte-
* gers in the integer vector ivec. Handels only minus signs and digits.
      character str*(*),pst*70,dig*1 
      integer ivec(10),nin,lt,np,int,k,i,j,l,itg,len
      common/dgt/dig(10)
      nin=0
      lt=len(str)
      np=1
      do 4 i=1,lt
         pst=str(np:lt)
         int=0
         k=0
         do 2 j=1,9
            itg=-1
            do 1 l=1,10
               if(pst(j:j).eq.dig(l)) itg=l-1
 1          continue
            if(itg.eq.-1) go to 3
            k=k+1
            int=int*10+itg
 2       continue
         write(6,*) ' WARNING from intstr, integer string too long.'
 3       continue
         if(k.ne.0) then
            j=max(np-1,1)
            if(str(j:j).eq.'-') int=-int
            nin=nin+1
            ivec(nin)=int
            if(nin.eq.10) return
            np=np+k
         endif
         np=np+1
         if(np.gt.lt) return
 4    continue
      return
      end
*===================================================================
* BELOW SOME BASIC GEOMETRIC SUBROUTINES: SCP,ABV,VPR,BIJ AND REFSYS
*------------------------------------------------------------------
      double precision function scp(v1,v2)
* Scalar product of vectors v1 and v2.
      double precision v1(3),v2(3)
      scp=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
      return
      end
*-------------------------------------------------------------------
      double precision function abv(v1)
* The absolute value i.e. length of the vector v1.
      double precision v1(3),scp
      abv=sqrt(scp(v1,v1))
      return
      end
*-------------------------------------------------------------------
      subroutine vpr(n,v1,v2,v3)
* Vector product of v1 and v2 in v3. If n=1 v3 is normalized.
      double precision v1(3),v2(3),v3(3),dd,abv,eps
      integer n,i
      parameter (eps = 1.D-9)
      v3(1)=v1(2)*v2(3)-v1(3)*v2(2)
      v3(2)=v1(3)*v2(1)-v1(1)*v2(3)
      v3(3)=v1(1)*v2(2)-v1(2)*v2(1)
      if(n.ne.1) return
      dd=abv(v3)
      if(dd.lt.eps) return
      do 1 i=1,3
 1       v3(i)=v3(i)/dd
      return
      end
*--------------------------------------------------------------------
      subroutine bij(nx,a,i,j,v1,dd,ierr)
* Calculates the unit vector v1 and the distance dd from i to j .
      double precision a(nx,3),v1(3),dd,abv,eps
      integer nx,i,j,ierr,k
      parameter (eps = 1.D-4)
      v1(1)=a(j,1)-a(i,1)
      v1(2)=a(j,2)-a(i,2)
      v1(3)=a(j,3)-a(i,3)
      dd=abv(v1)
      if(dd.lt.eps) then
         if(ierr.eq.0) write(6,1000) i,j
 1000    format('  Atoms',i3,'  and',i3,'   coincide !!')
         ierr=ierr+1
         return
      endif
      do 1 k=1,3
 1       v1(k)=v1(k)/dd
      return
      end
*--------------------------------------------------------------------
      subroutine refsys(nx,a,i,j,k,v1,v2,v3,sc,d1,d2,ierr)
* Calculates the orthonormal basis v1-v3 defined by atoms i,j,k.
      double precision a(nx,3),v1(3),v2(3),v3(3),sc,d1,d2,scp,one,alm1
      integer nx,i,j,k,ierr,iold
      parameter (one = 1.D0, alm1= 0.999d0)
      iold=ierr
      call bij(nx,a,i,j,v1,d1,ierr)
      if(ierr.gt.iold) return
      call bij(nx,a,i,k,v2,d2,ierr)
      if(ierr.gt.iold) return
      sc=scp(v1,v2)
      sc=dsign(min(abs(sc),one),sc)
      if(abs(sc).ge.alm1) then
         if(ierr.eq.0) write(6,1000) i,j,k
 1000    format(' Atoms ',i3,'  and ',i3,'  and ',i3,'  on line!!')
         ierr=ierr+1
         return
      endif
      call vpr(1,v1,v2,v3)
      call vpr(0,v3,v1,v2)
      return
      end
*======================================================================
*      OTHER GEOMETRICAL SUBROUTINES: 
*   CYLI, SPHERI, ROT, QPR, AND GPRM
*----------------------------------------------------------------------
      subroutine cyli(x,y,r,ph,ic)
      implicit double precision(a-h,o-z)
      common/pp/conv
      if(ic.ne.1) go to 10
      x=r*dcos(conv*ph)
      y=r*dsin(conv*ph)
      return
   10 r=dsqrt(x*x+y*y)
      ph=0.D0
      if(r.lt.1.D-10) return
      if(x.eq.0.D0) ph=dsign(90.D0,y)
      if(x.ne.0.D0) ph=90.D0-dsign(90.D0,x)+datan(y/x)/conv
      if(ph.gt.180.D0) ph=ph-360.D0
      return
      end
*---------------------------------------------------------------------
      subroutine spheri(x,y,z,r,ph,th,ic)
      implicit double precision(a-h,o-z)
      common/pp/conv
      if(ic.ne.1) go to 10
      call cyli(x,y,r,ph,1)
      s=dsin(th*conv)
      x=x*s
      y=y*s
      z=r*dcos(th*conv)
      return
   10 call cyli(x,y,r,ph,-1)
      if(r.ne.0.D0) th=90.D0-datan(z/r)/conv
      if(r.eq.0.D0) th=90.D0-dsign(90.D0,z)
      r=dsqrt(r*r+z*z)
      if(r.lt.1.D-10) th=0.D0
      return
      end
*--------------------------------------------------------------------
      subroutine rot(qr,qv,qw)
*  Rotation is carried out with quaternion multiplication.
      double precision qr(4),qt(4),qv(4),qw(4),qu(4)
      call qpr(qr,qv,qu)
      qt(1)=qr(1)
      do 1 i=1,3
 1       qt(i+1)=-qr(i+1)
      call qpr(qu,qt,qw)
      return
      end
*--------------------------------------------------------------------
      subroutine qpr(a,b,c)
*  c is the product of the quaternions a and b; a*b=c.
      double precision a(4),b(4),c(4),v1(3),v2(3),v3(3)
      c(1)=a(1)*b(1)
      do 1 i=1,3
         v1(i)=a(i+1)
         v2(i)=b(i+1)
 1       c(1)=c(1)-v1(i)*v2(i)
      call vpr(0,v1,v2,v3)
      do 2 i=1,3
 2       c(i+1)=a(1)*v2(i)+b(1)*v1(i)+v3(i)
      return
      end
*--------------------------------------------------------------------
      double precision function gprm(nx,a,i,j,k,l)
      implicit double precision(a-h,o-z)
*  Calculates the distance r(i,j)=gprm(i,j,0,0) or the angle
*  the(i,j,k)=gprm(i,j,k,0) or dihedral angle:gprm(i,j,k,l).
      dimension a(nx,3),v1(3),v2(3),v3(3),v4(3),v5(3),v6(3)
      common/pp/conv
      ierr=1
      if(l.ne.0) go to 3
      if(k.ne.0) go to 2
      call bij(nx,a,i,j,v1,d1,ierr)
      gprm=d1
      return
 2    sc=1.D0
      call refsys(nx,a,j,i,k,v1,v2,v3,sc,d1,d2,ierr)
      gprm=dacos(sc)/conv
      return
 3    call refsys(nx,a,j,k,l,v1,v2,v3,sc,d1,d2,ierr)
      gprm=-360.D0
      if(ierr.gt.1) return
      call refsys(nx,a,j,k,i,v4,v5,v6,sp,d1,d2,ierr)
      if(ierr.gt.1) return
      call vpr(0,v3,v6,v2)
      sc=scp(v3,v6)
      sp=scp(v1,v2)
      if(sp.eq.0.D0) sp=-1.D-18
      if(sc.ne.0.D0) sc=dsign(dmin1(dabs(sc),1.D0),sc)
      gprm=dacos(sc)/conv
      if(gprm.ne.0.D0) gprm=-dsign(gprm,sp)
      return
      end
*===================================================================
* DATA AND SUBROUTINES FOR HANDLING OF CHARGES AND CHEMICAL SYMBOLS:
* BLOCK DATA, ATMNME, SYMB, SYMSTR, SYMSTZ AND NQQ
*--------------------------------------------------------------------
      block data
      character isymv*2,dig*1
      common/atn/isymv(-9:99)
      common/dgt/dig(10)
      common/dfu/lstdha,jt(3)
      save /dfu/
      data isymv/
     +'! ','# ','$ ','% ','& ','* ','? ','@ ','X ','  ',
     +'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne',
     +'Na','Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti',
     +'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se',
     +'Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd',
     +'Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs','Ba','La','Ce',
     +'Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb',
     +'Lu','Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg',
     +'Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th',
     +'Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es'/
      data dig/'0','1','2','3','4','5','6','7','8','9'/
      data lstdha/0/
      end
*-------------------------------------------------------------------
      character*2 function atmnme(nqq)
      character isymv*2
      common/atn/isymv(-9:99)
      if(nqq.le.-10) atmnme=isymv(-9)
      if(nqq.ge.-9.and.nqq.le.99) atmnme=isymv(nqq)
      if(nqq.ge.100) atmnme=isymv(99)
      return
      end
*--------------------------------------------------------------------
      character*5 function symb(xx,nz)
* Takes the chemical symbol xx and n(Z) atom-type count number nz and
* makes a symbol of these which is packed to the left.
      character xx*2,numbst*3
      integer nz
      if(xx(2:2).eq.' ') then
         symb=xx(1:1)//numbst(nz)
      else
         symb=xx//numbst(nz)
      endif
      return
      end
*--------------------------------------------------------------------
      character*3 function numbst(nz)
* Takes the number nz and makes the corresponding
* character string which is packed to the left.
      character dig*1,sdig(3)*1
      integer nz,nnz,ndig(3),i,j
      common/dgt/dig(10)
      numbst='   '
      if(nz.ge.1000.or.nz.lt.0) return
      nnz=nz
      ndig(1)=nnz/100
      nnz=nnz-100*ndig(1)
      ndig(2)=nnz/10
      nnz=nnz-10*ndig(2)
      ndig(3)=nnz
      do 1 i=1,3
         do 2 j=1,10
            if(ndig(i).eq.j-1) sdig(i)=dig(j)
 2       continue
 1    continue
      numbst=sdig(3)
      if(nz.ge.10) numbst=sdig(2)//sdig(3)
      if(nz.ge.100) numbst=sdig(1)//sdig(2)//sdig(3)
      return
      end
*--------------------------------------------------------------------
      subroutine symstr(nx,nq,nnq,n1,n2,n3,n4,sstr)
      integer nq(nx),nnq(nx),nx,n1,n2,n3,n4,last
      character sstr*23, symb*5, atmnme*2
      sstr='  '
      if(n1.le.0) return
      sstr=symb(atmnme(nq(n1)),nnq(n1))
      if(n2.le.0) return
      call strip(sstr,last)
      sstr=sstr(1:last)//'-'//symb(atmnme(nq(n2)),nnq(n2))
      if(n3.le.0) return
      call strip(sstr,last)
      sstr=sstr(1:last)//'-'//symb(atmnme(nq(n3)),nnq(n3))
      if(n4.le.0) return
      call strip(sstr,last)
      sstr=sstr(1:last)//'-'//symb(atmnme(nq(n4)),nnq(n4))
      return
      end
*--------------------------------------------------------------------
      character*20 function symstz(nx,nq,nnq,n1,n2,n3,n4)
      integer nq(nx),nnq(nx),nx,n1,n2,n3,n4,last
      character symb*5, atmnme*2
      symstz='  '
      if(n1.le.0) return
      symstz=symb(atmnme(nq(n1)),nnq(n1))
      if(n2.le.0) return
      call strip(symstz,last)
      symstz=symstz(1:last)//symb(atmnme(nq(n2)),nnq(n2))
      if(n3.le.0) return
      call strip(symstz,last)
      symstz=symstz(1:last)//symb(atmnme(nq(n3)),nnq(n3))
      if(n4.le.0) return
      call strip(symstz,last)
      symstz=symstz(1:last)//symb(atmnme(nq(n4)),nnq(n4))
      return
      end
*--------------------------------------------------------------------
      integer function nqq(atmn)     
      character atmn*2,isymv*2
      common/atn/isymv(-9:99)
      do 1 nqq=-9,99
         if(atmn.eq.isymv(nqq)) return
 1    continue
      return
      end
*====================================================================
* SUBROUTINES TREFST, REFNUC AND NEIGBR FOR Z-MATRIX CONSTRUCTION
*---------------------------------------------------------------------
      subroutine trefst(nx,ifx,n1,n2,n3)
* Finds the numbers of the three first entered atoms.
      logical ifx(nx)
      integer nx,n1,n2,n3,i
      n1=0
      n2=0
      n3=0
      do 1 i=1,nx
         if(ifx(i)) then
            if(n1.eq.0) then
               n1=i
               go to 1
            elseif(n2.eq.0) then
               n2=i
               go to 1
            else
               n3=i
               go to 2
            endif
         endif
 1    continue
 2    continue
      return
      end
*----------------------------------------------------------------------
      subroutine refnuc(nx,ifx,nq,a,n3,n,i,j,k)
* Calculates the three reference atoms i,j,k of atom n. n3 is the 
* third entered atom.  neigbr is used to find the four nearest atoms.
* Preferably the angles the(n,i,j) and the(i,j,k) are not 0 or 180.
* Hydrogens are avoided as middle (i,j) reference atoms.
      logical ifx(nx)
      integer nq(nx),nx,n,i,j,k,ja,ierr,n3,iv(12),nv
      double precision a(nx,3),v1(3),v2(3),v3(3),d1,d2,sc
      neig=4
      call neigbr(nx,ifx,a,neig,n,n,n,iv,nv)      
      if(n.eq.n3) then
         i=iv(1)
         j=iv(2)
         if(nq(i).eq.1.and.nq(j).ne.1) then
            i=iv(2)
            j=iv(1)
         endif
         return
      else
         do 1 ja=1,nv
            i=iv(ja)
            if(nq(i).ne.1) go to 2
 1       continue
         i=iv(1)
 2       continue
      endif
      call neigbr(nx,ifx,a,neig,n,i,i,iv,nv)
      do 3 ja=1,nv
         j=iv(ja)
         ierr=1
         call refsys(nx,a,n,i,j,v1,v2,v3,sc,d1,d2,ierr)
         if(nq(j).ne.1.and.ierr.ne.2) go to 5
 3    continue
      do 4 ja=1,nv
         j=iv(ja)
         ierr=1
         call refsys(nx,a,n,i,j,v1,v2,v3,sc,d1,d2,ierr)
         if(nq(j).ne.1.or.ierr.ne.2) go to 5
 4    continue
      j=iv(1)
 5    continue      
      call neigbr(nx,ifx,a,neig,n,j,i,iv,nv)
      do 6 ja=1,nv
         k=iv(ja)
         ierr=1
         call refsys(nx,a,i,j,k,v1,v2,v3,sc,d1,d2,ierr)
         if(ierr.ne.2) go to 7
 6    continue
      k=iv(1)
 7    continue
      return
      end
*---------------------------------------------------------------------
      subroutine neigbr(nx,ifx,a,neig,mx,n,na,iv,nv)
* Finds the neig (max 12) atoms nearest to atom n and stores their numbers in
* the vector iv. The number of neighbours found and stored is nv. The near-
* est atoms are first in iv. If equally near the lower numbers are first.
* Search among atom numbers less than mx; avoid na.
      logical ifx(nx)
      integer iv(12),nx,mx,n,na,ia,nv,ii,jj,jmax
      double precision a(nx,3),da,d1,d2,d3,dii,dq(12),big,eps
      parameter (big = 1.D10, eps = 1.D-4)
      do 1 ii=1,neig
         dq(ii)=big
 1    continue
      nv=0
      do 6 ii=1,mx-1
         if(ifx(ii).and.ii.ne.n.and.ii.ne.na) then
            d1=a(n,1)-a(ii,1)
            d2=a(n,2)-a(ii,2)
            d3=a(n,3)-a(ii,3)
            dii=d1*d1+d2*d2+d3*d3
            jmax=min0(nv+1,neig)
            do 3 jj=1,jmax
               if(dii.le.dq(jj)) then
                  nv=min0(nv+1,neig)
                  dq(nv)=dii-eps
                  iv(nv)=ii
                  do 2 ia=1,nv-1
                     dq(ia)=dq(ia)-eps
 2                continue
                  go to 4
               endif
 3          continue
 4          continue
            do 5 jj=1,nv-1
               if(dq(jj).gt.dq(nv)) then
                  da=dq(jj)
                  dq(jj)=dq(nv)
                  dq(nv)=da
                  ia=iv(jj)
                  iv(jj)=iv(nv)
                  iv(nv)=ia
               endif
 5          continue
         endif
 6    continue
      return
      end
*=====================================================================
* SUBROUTINES ORDER AND ORDPAR,
* WHICH ORDER ATOMS AFTER DESIRED COORDINATE-TYPE,   AND  CLEAN
*---------------------------------------------------------------------
      subroutine order(nx,ifx,nq,a,mf,ml,it,ifz)
      implicit double precision(a-h,o-z)
      logical ifx(nx),ifz
      dimension w(3),nq(nx),a(nx,3)
      parameter (eps=1.d-5)
      do 50 k=mf,ml-1
         if(.not.ifx(k)) go to 50
         call ordpar(nx,ifx,a,it,k,rk,tk,pk)
         do 10 i=k+1,ml
            if(.not.ifx(i)) go to 10
            call ordpar(nx,ifx,a,it,i,ri,ti,pi)
            if(ri.gt.rk+eps) go to 10
            if(ri.lt.rk-eps) go to 15
            if(ti.gt.tk+eps) go to 10
            if(ti.lt.tk-eps) go to 15
            if(pi.gt.pk+eps) go to 10
 15         continue
            if(ifz) then
               nqq=nq(k)
               nq(k)=nq(i)
               nq(i)=nqq
            endif
            do 20 j=1,3
               w(j)=a(k,j)
               a(k,j)=a(i,j)
               a(i,j)=w(j)
 20         continue
            rk=ri
            pk=pi
            tk=ti
 10      continue
 50   continue
      return
      end
*---------------------------------------------------------------------
      subroutine ordpar(nx,ifx,a,it,n,rn,tn,pn)
      double precision a(nx,3),rn,tn,pn
      integer nx,it,n
      logical ifx(nx)
      if(it.ge.2) go to 2
      rn=a(n,1)
      tn=a(n,2)
      pn=a(n,3)
      return
    2 if(it.eq.3) go to 3
      call cyli(a(n,1),a(n,2),tn,pn,-1)
      rn=-a(n,3)
      if(pn.lt.0.d0) pn=pn+360.d0
      return
    3 call spheri(a(n,1),a(n,2),a(n,3),rn,pn,tn,-1)
      if(pn.lt.0.d0) pn=pn+360.d0
      return
      end
*--------------------------------------------------------------------
      subroutine clean(nx,ifx,a,nf,nl)
* Clean attempts to remove small numerical errors resulting from round 
* off when z-matrices are read.
      implicit double precision(a-h,o-z)
      logical ifx(nx)
      dimension a(nx,3)
      eps=1.d-4
      do 10 i=nf,nl
         if(.not.ifx(i)) go to 10
         do 20 j=1,3
 20         if(dabs(a(i,j)).lt.eps) a(i,j)=0.d0
         call cyli(a(i,1),a(i,2),r,phi,-1)
         ep=-eps
         if(phi.ge.0.d0) ep=eps
         ph=dint(phi+ep)
         if(dabs(ph-phi).lt.2.d0*eps) call cyli(a(i,1),a(i,2),r,ph,1)
 10   continue
      do 30 i=nf,nl-1
         if(.not.ifx(i)) go to 30
         r=dsqrt(a(i,1)**2+a(i,2)**2+a(i,3)**2)
         do 40 j=i+1,nl
            if(.not.ifx(j)) go to 40
            call spheri(a(j,1),a(j,2),a(j,3),rr,ph,th,-1)
            if(dabs(r-rr).lt.eps) 
     +           call spheri(a(j,1),a(j,2),a(j,3),r,ph,th,1)
 40      continue
 30   continue
      return
      end
*=====================================================================
* THE FOLLOWING ROUTINES FIND THE EIGENVALUES AND EIGENVECTORS OF A
* REAL SYMMETRIC 3X3 MATRIX BY DIRECT SOLUTION OF THE CUBIC SECULAR EQ. 
*---------------------------------------------------------------------
      subroutine eigval(x,y,w,ndeg)
*  Calculates the eigenvalues w(3) of a 3*3 real symmetric
*  matrix x(3)=(xx,yy,zz),y(3)=(xy,yz,zx).
      integer ndeg,i
      double precision x(3),y(3),w(3),ww,scp,a,b,c,q,r,arg,cc,ss,dda,ddb
      double precision  zero,onesxt,two,three, eps4,eps6,eps7
      parameter (zero=0.d0, onesxt=(1.d0/6.d0), two=2.d0, three=3.d0)
      parameter (eps4=2.d-4, eps6=1.d-6, eps7=1.d-7)
      do 10 i=1,3
         if(abs(x(i)).lt.eps7) x(i)=zero
 10      if(abs(y(i)).lt.eps6) y(i)=zero
*
      a= x(1)*y(2)*y(2) + x(2)*y(3)*y(3) + x(3)*y(1)*y(1)
     +        - x(1)*x(2)*x(3) - two*y(1)*y(2)*y(3)
      b= x(1)*x(2) + x(2)*x(3) + x(3)*x(1) - scp(y,y)
      c= -(x(1)+x(2)+x(3))/three
*
      q= b/three - c*c
      if(abs(q).lt.eps6) q=zero
      r= (b*c - a)/two - c*c*c
      if(abs(r).lt.eps7) r=zero
*
      a= sqrt(abs(q*q*q + r*r))
      arg= acos(zero)/three
      if(r.ne.zero) arg =atan(a/r)/three + arg - dsign(arg,r)
      q= (r*r + a*a)**onesxt
      cc= cos(arg)
      ss= sin(arg)
*
      w(1)= -c - q*(cc + sqrt(three)*ss)
      w(2)= -c - q*(cc - sqrt(three)*ss)
      w(3)= -c + q*two*cc
*
      ndeg=0
      dda=abs(w(1)-w(2))
      ddb=abs(w(2)-w(3))
      if(dda+ddb.lt.eps4+eps4) ndeg=3
      if(ndeg.eq.3) return
      if(dda.lt.eps4) ndeg=2
      if(ndeg.eq.2) return
      if(ddb.lt.eps4) ndeg=2
      ww=w(3)
      w(3)=w(1)
      w(1)=ww
      return
      end
*---------------------------------------------------------------------
      subroutine eigvec(w,xm,ym,ev,ier)
      double precision xm(3),ym(3),v1(3),v2(3),ev(3),w,scp,alm1
      integer ier,it
      parameter (alm1 = 1.d0 - 1.d-10)
* Finds the eigenvector ev corresponding to the eigenvalue w.
      v1(1)=xm(1)-w
      v1(2)=ym(1)
      v1(3)=ym(3)
      v2(1)=ym(1)
      v2(2)=xm(2)-w
      v2(3)=ym(2)
      it=1
    1 call vpr(1,v1,v2,ev)
      ier=0
      if(scp(ev,ev).gt.alm1) return
      if(it.eq.2) go to 2
      ier=1
      if(it.eq.3) return
      v2(1)=ym(3)
      v2(2)=ym(2)
      v2(3)=xm(3)-w
      it=2
      go to 1
    2 v1(1)=ym(1)
      v1(2)=xm(2)-w
      v1(3)=ym(2)
      it=3
      go to 1
      end
*=====================================================================
* The two routines LUDCMP and LUBKSB below are the basic linear equa-
* tion solvers and are adapted from  W.H.Press et al., Numerical Reci-
* pes (Cambridge, 1986). Ludcmp takes the input coefficient matrix a
* and turns it into LU (lower and upper triangular) form after pivot-
* ing (stored in indx). Lubksb turns the vector of right hand side
* numbers b into the solution vector.
*----------------------------------------------------------------------
      subroutine ludcmp(a,indx,ierr)
      implicit double precision(a-h,o-z)
      parameter (tiny=1.d-20)
      dimension a(3,3),indx(3),vv(3)
      ierr=1
      do 12 i=1,3
         aamax=0.
         do 11 j=1,3
            if(abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
 11      continue
         if(aamax.eq.0.) return
         vv(i)=1./aamax
 12   continue
      do 19 j=1,3
         do 14 i=1,j-1
            sum=a(i,j)
            do 13 k=1,i-1
               sum=sum-a(i,k)*a(k,j)
 13         continue
            a(i,j)=sum
 14      continue
         aamax=0.
         do 16 i=j,3
            sum=a(i,j)
            do 15 k=1,j-1
               sum=sum-a(i,k)*a(k,j)
 15         continue
            a(i,j)=sum
            dum=vv(i)*abs(sum)
            if(dum.ge.aamax) then
               imax=i
               aamax=dum
            endif
 16      continue
         if(j.ne.imax) then
            do 17 k=1,3
               dum=a(imax,k)
               a(imax,k)=a(j,k)
               a(j,k)=dum
 17         continue
            vv(imax)=vv(j)
         endif
         indx(j)=imax
         if(a(j,j).eq.0) a(j,j)=tiny
         if(j.ne.3) then
            dum=1./a(j,j)
            do 18 i=j+1,3
               a(i,j)=a(i,j)*dum
 18         continue
         endif
 19   continue
      ierr=0
      return
      end
*---------------------------------------------------------------------
      subroutine lubksb(a,indx,b)
      implicit double precision(a-h,o-z)
      dimension a(3,3),indx(3),b(3)
      ii=0
      do 12 i=1,3
         ll=indx(i)
         sum=b(ll)
         b(ll)=b(i)
         if(ii.ne.0) then
            do 11 j=ii,i-1
               sum=sum-a(i,j)*b(j)
 11         continue
         elseif(sum.ne.0.) then
            ii=i
         endif
         b(i)=sum
 12   continue
      do 14 i=3,1,-1
         sum=b(i)
         if(i.lt.3) then
            do 13 j=i+1,3
               sum=sum-a(i,j)*b(j)
 13         continue
         endif
         b(i)=sum/a(i,i)
 14   continue
      return
      end
*=====================================================================
* Here follows the routine AMOEBA, also from Press et al. (modified)
* It finds the minimum of a scalar function in 3 dimensional space.
*--------------------------------------------------------------------
      subroutine amoeba(p,y,ftol,iter,ierr)
      implicit double precision(a-h,o-z)
      parameter (alpha=1.d0, beta=0.5d0, gamma=2.0d0, itmax=500)
      dimension p(4,3),y(4),pr(3),prr(3),pbar(3)
      iter=0
 1    ilo=1
      if(y(1).gt.y(2))then
         ihi=1
         inhi=2
      else
         ihi=2
         inhi=1
      endif
      do 11 i=1,4
         if(y(i).lt.y(ilo)) ilo=i
         if(y(i).gt.y(ihi)) then
            inhi=ihi
            ihi=i
         elseif(y(i).gt.y(inhi)) then
            if(i.ne.ihi) inhi=i
         endif
 11   continue
      rtol=abs(y(ihi)-y(ilo))
      if(rtol.lt.ftol) return
      if(iter.eq.itmax) return
      iter=iter+1
      do 12 j=1,3
         pbar(j)=0.d0
 12   continue
      do 14 i=1,4
         if(i.ne.ihi) then
            do 13 j=1,3
               pbar(j)=pbar(j)+p(i,j)
 13         continue
         endif
 14   continue
      do 15 j=1,3
         pbar(j)=pbar(j)/3
         pr(j)=(1.d0+alpha)*pbar(j)-alpha*p(ihi,j)
 15   continue
      ypr=func(pr,ierr)
      if(ierr.ne.0) return
      if(ypr.le.y(ilo)) then
         do 16 j=1,3
            prr(j)=gamma*pr(j)+(1.d0-gamma)*pbar(j)
 16      continue
         yprr=func(prr,ierr)
         if(ierr.ne.0) return
         if(yprr.lt.y(ilo)) then
            do 17 j=1,3
               p(ihi,j)=prr(j)
 17         continue
            y(ihi)=yprr
         else
            do 18 j=1,3
               p(ihi,j)=pr(j)
 18         continue
            y(ihi)=ypr
         endif
      elseif(ypr.ge.y(inhi)) then
         if(ypr.lt.y(ihi)) then
            do 19 j=1,3
               p(ihi,j)=pr(j)
 19         continue
            y(ihi)=ypr
         endif
         do 21 j=1,3
            prr(j)=beta*p(ihi,j)+(1.d0-beta)*pbar(j)
 21      continue
         yprr=func(prr,ierr)
         if(ierr.ne.0) return
         if(yprr.lt.y(ihi)) then
            do 22 j=1,3
               p(ihi,j)=prr(j)
 22         continue
            y(ihi)=yprr
         else
            do 24 i=1,4
               if(i.ne.ilo) then
                  do 23 j=1,3
                     pr(j)=0.5d0*(p(i,j)+p(ilo,j))
                     p(i,j)=pr(j)
 23               continue
                  y(i)=func(pr,ierr)
                  if(ierr.ne.0) return                  
               endif
 24         continue
         endif
      else
         do 25 j=1,3
            p(ihi,j)=pr(j)
 25      continue
         y(ihi)=ypr
      endif
      go to 1
      end
