*======================================================================
      Program EVCLID
*======================================================================
*     A Geometric Calculator and Editor
*     by Hanno Essen, revision of March 1992 (KTH, Stockholm)
*     originally developed at Mcmaster University,Hamilton,Ont.,Canada,1982
*----------------------------------------------------------------------
*     Declarations:
      logical ifx(+360), ifz,  ifh, ifc
* ifx(n)=.true. if atom number n has been entered. 
* ifz=.true. if atomic numbers Z have been entered.
* ifh=.true. if help command was read by routine readin
* ifc=.true. if a command was found in readin
      integer nq(+360), nnq(+360), mz(+360,3), imz(+360), nxz
      integer nx, lnmx,  lfn,mxn,noa, nnh, nmz,  last, i, maccnt
* nq(n) contains atomic number Z of atom n.
* nnq(n) contains a count of atoms with the same Z.
* mz(n,3) contains the z-matrix reference atoms i,j,k of atom n
* imz(n)= 0 if no z-matrix info for atom n, 1 if normal, -1 if alternative
      double precision a(+360,3), conv, cnzs(+361)
* a(n,3) contains the cartesian coordinates of entered atom number n.
* conv is the conversion factor degrees to radians. 
* cnzs(n) contains numerical z-matrix constants of symbolic z-matrix
      character comm*4,x*1,yz*2, str*80,prmp*79,numbst*3
* symbs contains the symbols in symbolic z-matrices, consts the numbers
      character symbs(+361)*20, consts(+361)*25
*--------------------------------------------------------------------
      common/pp/conv
*---------------------------------------------------------------------
* nx= basic array dimension (max no of atoms). Change above and here :
      nx=+360
* lnmx= max no of output lines on screen.
      lnmx=22
* nxz= max no of symbols in symbolic z-matrix (used in izmcom and ozmcom)
      nxz=+361
*--------------------------------------------------------------------
      prmp=' Evclid, a geometric calculator and editor by H. Essen'
      prmp(68:79)='version:1995'
      write(6,*) prmp
* Initialize conv and arrays nq, ifx and imz
      conv=dacos(-1.D0)/180.D0
      mz(3,3)=0
      do 1 i=1,nx
         nnq(i)=i
         nq(i)=1
         ifx(i)=.false.
         imz(i)=0
 1    continue
* Initialize ifc,  lfn, mxn, noa and nmz, and nnh, ifz
      ifc=.false.
      call naccnt(nx,ifx,lfn,mxn,noa)
      nmz= maccnt(nx,imz)
      call caccnt(nx,ifx,nq,nnq,mxn,nnh,ifz)
*     
      prmp=' Type:command      (h = help,  hh = basic help)'
      go to 93
* Come back here if command was not recognized
 90   continue
*     
* Return to this point after visit to option subroutine or help-writeout
 91   continue
* Command returned from subroutine readin is in  str
      if(ifc) then
         comm=str(1:4)
         ifc=.false.
         go to 80
      endif
      prmp='Type:command'
      if(noa.gt.0) then
         prmp=prmp(1:12)//'   (no of atoms:'//numbst(noa)
         call strip(prmp,last)
         prmp=prmp(1:last)//')'
      endif
      if(nnh.gt.0) then
         call strip(prmp,last)
         prmp=prmp(1:last)//'  (no of Z >1: '//numbst(nnh)
         call strip(prmp,last)
         prmp=prmp(1:last)//')'
      endif
      if(nmz.gt.0) then
         call strip(prmp,last)
         prmp=prmp(1:last)//' (z-mtx rows:'//numbst(nmz)
         call strip(prmp,last)
         prmp=prmp(1:last)//')'
      endif
      prmp=' '//prmp
 93   prmp(75:79)='?xyz?'
      write(6,*) prmp
 930  str(1:4)='Evc>'
      call readin(ifc,str,last,ifh)
      if(str.eq.'file-end') go to 930
      if(str.eq.'assgmnt=') go to 93
      ifc=.false.
      if(last.gt.4) go to 92
      if(last.eq.0) then
         if(noa.eq.0) then
            write(6,*) ' To terminate run type: quit'
            go to 90
         elseif(noa.ne.0) then
            str='dco '
            last=3
         endif
      endif
      comm=str(1:last)
      if(comm.eq.'help') comm='hh '
*     
* Put non-standard abbreviations on standard form
      if(comm.eq.'c   ') comm='ico '
      if(comm.eq.'s   ') then
         if(noa.gt.0) comm='osa '
         if(noa.eq.0.and.(.not.ifz).and.nmz.eq.0) comm='isa '
         if(comm.eq.'s   ') then
            write(6,*) 
     +           ' Zero atoms but not empty.  Specify: isa  or: osa !?'
            go to 91
         endif
      endif
      if(comm.eq.'r   ') comm='erm '
      if(comm(4:4).ne.' ') go to 92
*
 80   x=comm(1:1)
      yz=comm(2:3)
      if(noa.eq.0.and.(x.ne.'i'.and.x.ne.'h'.and.yz.ne.'h ')) 
     +     write(6,*)' When zero atoms option should be of input-type.'
*
* When the command xyz has been given a subroutine with the name
*            subroutine xyzcom
* is called. It issues a prompt for input parameters and when its task
* is finished control is returned to the main program (top level).
      if(x.eq.'i') then
         if(yz.eq.'  '.or.yz.eq.'ge'.or.yz.eq.'g ') then
            call igecom(str,ifc,nx,ifx,a,noa)
         elseif(yz.eq.'co'.or.yz.eq.'c ') then
            call icocom(str,ifc,nx,ifx,a,lfn,mxn,noa,nmz,imz)
         elseif(yz.eq.'po'.or.yz.eq.'p ') then
            call ipocom(str,ifc,nx,ifx,nq,a,mxn,ifz)
         elseif(yz.eq.'an'.or.yz.eq.'a ') then
            call iancom(str,ifc,nx,ifx,nq,noa,ifz)
         elseif(yz.eq.'sa') then
            call isacom(str,ifc,nx,ifx,nq,a,noa,mxn,ifz,nnh,nmz,imz,mz)
         elseif(yz.eq.'cf') then
            call icfcom(str,ifc,nx,ifx,nq,a,ifz,mxn)
         elseif(yz.eq.'zm'.or.yz.eq.'z ') then
            call izmcom(str,ifc,nx,ifx,nq,a,ifz,
     +           imz,mz,mxn,symbs,nxz,consts)
         elseif(ifh) then
            write(6,1001)
            write(6,1002)
         else
            go to 92
         endif
*
      elseif(x.eq.'d')then
         if(yz.eq.'  '.or.yz.eq.'co'.or.yz.eq.'c ') then
            call dcocom(str,ifc,nx,ifx,nq,nnq,a,mxn,ifz,lnmx)
         elseif(yz.eq.'zm'.or.yz.eq.'z ') then
            call dzmcom(str,ifc,nx,ifx,nq,nnq,a,mxn,ifz,lnmx,nmz,imz,mz)
         elseif(yz.eq.'ge'.or.yz.eq.'g ') then
            call dgecom(str,ifc,nx,ifx,nq,nnq,a,mxn,ifz,lnmx)
         elseif(ifh) then
            write(6,1003)
         else
            go to 92
         endif
*
      elseif(x.eq.'o')then
         if(yz.eq.'  '.or.yz.eq.'zm'.or.yz.eq.'z ') then
            if(.not.ifz) write(6,*)' ** Don''t forget atomic charges **'
            call ozmcom(str,ifc,nx,ifx,nq,nnq,a,mxn,
     +           imz,mz,symbs,nxz,cnzs)
         elseif(yz.eq.'cf'.or.yz.eq.'co'.or.yz.eq.'c ') then
            call ocfcom(str,ifc,nx,ifx,nq,nnq,a,mxn,ifz)
         elseif(yz.eq.'sa') then
            call osacom(nx,ifx,nq,a,noa,mxn,ifz,nnh,nmz,imz,mz)
         elseif(ifh) then
            write(6,1004)
            write(6,1024)
         else
            go to 92
         endif
*
      elseif(x.eq.'g')then
         if(yz.eq.'  '.or.yz.eq.'pr'.or.yz.eq.'p ') then
            call gprcom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'tr'.or.yz.eq.'t ') then
            call gtrcom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'ro'.or.yz.eq.'r ') then
            call grocom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'sc') then
            call gsccom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'re') then
            call grecom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'in') then
            call gincom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(ifh) then
            write(6,1051)
            write(6,1052)
         else
            go to 92
         endif
*
      elseif(x.eq.'e')then
         if(yz.eq.'  '.or.yz.eq.'mo'.or.yz.eq.'m ') then
            call emocom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'or'.or.yz.eq.'o ') then
            call eorcom(str,ifc,nx,ifx,nq,a,mxn,ifz)
         elseif(yz.eq.'pe'.or.yz.eq.'p ') then
            call epecom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'rr'.or.yz.eq.'r ') then
            call errcom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'rm'.or.yz.eq.'d ') then
            call ermcom(str,ifc,nx,ifx,nq,mxn,ifz,imz)
         elseif(yz.eq.'cp'.or.yz.eq.'c ') then
            call ecpcom(str,ifc,nx,ifx,nq,a,mxn)
         elseif(yz.eq.'zm'.or.yz.eq.'z ') then
            call ezmcom(str,ifc,nx,ifx,nq,a,mxn,nmz,imz,mz)
         elseif(ifh) then
            write(6,1006)
         else
            go to 92
         endif
*
      elseif(x.eq.'h') then
         if(yz.eq.'  ') then
            write(6,901)
            write(6,902) nx
         elseif(yz.eq.'h ') then
            write(6,1901)
            write(6,1902)
         elseif(yz.eq.'i ') then
            write(6,1001)
            write(6,1002)
         elseif(yz.eq.'d ') then
            write(6,1003)
         elseif(yz.eq.'o ') then
            write(6,1004)
            write(6,1024)
         elseif(yz.eq.'g ') then
            write(6,1051)
            write(6,1052)
         elseif(yz.eq.'e ') then
            write(6,1006)
         else
            go to 92
         endif
      else
         go to 92
      endif
*     
       if(x.eq.'i'.or.x.eq.'g'.or.x.eq.'e') then
          call naccnt(nx,ifx,lfn,mxn,noa)
          if(noa.eq.nx) write(6,*) ' Warning, atom storage full !'
          nmz= maccnt(nx,imz)
          call caccnt(nx,ifx,nq,nnq,mxn,nnh,ifz)
          if(noa.eq.0) go to 91
* Coordinates with very small absolut value 
* are assumed to be non-zero by accident, so they are set to zero.
          do 2 i=1,mxn
             if(.not.ifx(i)) go to 2
             if(dabs(a(i,1)).lt.1.D-8) a(i,1)=0.D0
             if(dabs(a(i,2)).lt.1.D-8) a(i,2)=0.D0
             if(dabs(a(i,3)).lt.1.D-8) a(i,3)=0.D0
 2        continue
       endif
       go to 91
* Go back and ask for new command
 92    write(6,*)' Command not recognized.  For help type: h'
       go to 90
* Go back to beginning and remind of the help option
*                                       below only format text  *end*
*---------------------------------------------------------------------
 901  format('              COMMAND MENU              '/
     +'   A *command* has the form:xyz  where  x  specifies the type'/
     +' of option and  yz  the specific option as follows: '/
     +'   x=i (Input)      yz=(ge) co po an h   (from file:) sa cf zm '/
     +'   x=d (Display)    yz=(co) zm ge h'/
     +'   x=o (Output)     yz=(zm) cf sa h              (default:(yz))'/
     +'   x=g (Geometric)  yz=(pr) tr ro sc re in h'/
     +'   x=e (Edit)       yz=(mo) or pe rr rm cp zm h'/
     +'   x=h (Help)       yz=(  ) i  d  o  g  e  h    (for more Help)'/
     +'  ABBREVIATIONS: (ido)  ge =Geometric par',
     +'ameters, co =Coordinates, po =Polygon'/
     +' or Polyhedron, an =Atomic Number, sa =Save',
     +', cf =Coordinate File, zm =Z-Matrix.')
 902  format(
     +' (g)   pr =Principal axes, tr =Translate,',
     +' ro =Rotate, sc =Scale, re =Reflect,  '/
     +' in =Invert.   (e)   mo =Move, or =Order, ',
     +'pe =Permute, rr =Remove & Renumber, '/
     +' rm =Remove, cp =Copy, zm =Z-Matrix.      ',
     +'          **  to stop Type: quit  **'//
     +'     General   NOTATION    in prompts:'/
     +' n  or  m,n  stand for Number(s) (max:',i3,
     +                                       ') assigned (or to be '/
     +'        assigned) to atoms at input by user or automatically.'/
     +' nf,nl  stand for a range of atoms these being the Numbers of'/    
     +'        the First and Last atoms in the range, respectively.'/
     +' ct     denotes Coordinate Type:x =Cartesian, c =Cylindrical '/
     +'        and s =Spherical (polar) coordinates.')
 1901 format('   BASIC HELP ABOUT EVCLID                            '/
     +'  The program works on an array of Cartesian coordinates of'/
     +' points. One can *Input  points to this array in various ways'/
     +' and once this has been done one can *Display  properties of'/
     +' the point set on the screen or *Output  such to file. One can'/
     +' manipulate the point set with *Geometric operations or with'/
     +' *Edit options (copy, remove, renumber). At the top level you'/
     +' can get *Help with one of the following commands:'/
     +'  h    for command menu           hi   for input-commands'/
     +'  hd   for display-commands       ho   for output-commands'/
     +'  hg   for geometric-commands     he   for edit-commands'/
     +'  hh   for the present text')
 1902 format(
     +'  Once you have selected and given a command:xyz you will get'/
     +' a prompt something like:'/
     +' Enter:a1 (a2)... etc.  (& a..)                (h = help) ',
     +'                *xyz*'/
     +' which tells you what arguments the command expects. Now you'/
     +' can enter: h for an explanation of what the argument symbols'/
     +' stand for. This information is maximum one screen long and'/
     +' then the prompt for input,  Enter:. etc., appears again.'/
     +' Note that the number of arguments is important and deter-'/
     +' mines the action taken. Input lines always end with return.'/
     +' Arguments in parenthesis are optional.')
 1001 format('         INPUT FROM TERMINAL      '/
     +'   Input of atoms from the terminal can be done using:'/
     +' Geometric parameters (distances, angles, dihedral angles)'/
     +' referencing already entered atoms;     command:ige  (i) '/
     +' or using Coordinates;                  command:ico  (c) '/
     +'   Make Polygon or Polyhedron with      command:ipo  (ip) '/
     +' Atomic Numbers can be entered with     command:ian  (ia)'/
     +'   When you have entered the command a prompt-line will appear'/
     +' which tells you what input will be accepted. If you still'/
     +' don''t understand you can type:h for an explanation.'/
     +'                           (Short command-form in parenthesis)')
 1002 format('         INPUT FROM FILE            '/
     +'   If the command:s (or:osa) has been given earlier to Save in-'/
     +' put this input will come back with     command:isa  (s)'/
     +' When the program is empty the command:s is taken to mean:isa'/
     +' otherwise it means:osa. Note that the command:isa removes all'/
     +' current contents before it puts back what was stored at save.'/
     +'   To read files containing:'/
     +' charges and Cartesian coordinates;     command:icf '/
     +' charges and Z-matrix;                  command:izm  (iz)'/
     +'   The z-matrix concept is explained if you give command:dzm'/
     +' and then type:h. (See also: ho for information on save.)')
 1003 format('            DISPLAY OPTIONS         '/
     +' Display of Coordinates,          command:dco  ( )'/
     +' Display of Z-matrix,             command:dzm  (dz)'/
     +' Display of Geometric parameters, command:dge  (dg)'/
     +'                          (Short command form in parenthesis) '/
     +'   After the command has been typed a prompt line will appear.'/
     +' No input is needed for the standard (default) display, just'/
     +' press return and the information will appear on the screen.'/
     +' As usual you can always type:h for help.'/
     +'   For coordinates the default is Cartesian coordinates'/
     +' (else you must give coordinate type:ct  as c or s).'/
     +'   A z-matrix is a set of 3N-6 geometric parameters that give'/
     +' the shape of the N-atom set in a coordinate independent way.'/
     +'   Lists of geometric parameters running over all pairs of'/
     +' atoms are generated with:dge. The default parameter type is'/
     +' distance. Give the command and type:h for further information.'/
     +'   One can always give a range: nf,nl to limit display.'/
     +' To cancel display you must type:q (just return will not work).')
 1004 format('         OUTPUT TO FILE               '/
     +'   Output options are     command:ozm  (o)    for Z-matrix'/
     +' and                      command:ocf  (oc)   for Coordinates'/
     +'                      (short command form in parenthesis)'/
     +'   You will be asked for a name for the output file the first'/
     +' time the command is given (it must be a new name).'/
     +'   The lines written to the file are also displayed on the'/
     +' screen. Each time the commands are called and lines appear on'/
     +' the terminal screen corresponding lines are added at the end'/
     +' of the output file (because no rewind is done).'/
     +'   In a subsequent run the output files can be read back in'/
     +' with: izm  or: icf  but the standard restart should be Save.')
 1024 format('         ABOUT SAVE'/
     +'   To save all your input to the program at a given time,'/
     +' just type: s  When there are entered atoms this is short for'/ 
     +'                          command:osa  (s)    for Save'/
     +' The output is written to file: ecldsv.dat and can be read back'/
     +' in later in the same or a subsequent run with command:isa.'/
     +' Each time you give the command s (or osa) output is written to'/
     +' the end of the file ecldsv.dat. When the program is empty  s'/
     +' will mean  isa  and the output of the last save is read in.'/
     +' Otherwise isa issues a prompt which asks which save you wish.')
 1051 format(
     +'   The default GEOMETRIC option is a transformation to a'/
     +' natural molecular coordinate-system: the origin is posi-'/
     +' tioned at the center of charge, and the axes are made Princi-'/
     +' pal axes i.e. eigenvectors of the (charge-)inertia tensor.'/
     +' In a symmetric molecule this system is often a natural sys-'/
     +' tem. A symmetry-axis becomes z-axis.'/
     +'                  Principal axes        command:gpr  (g) '/
     +'   The continuous transformations translation, rotation and'/
     +' scaling are optionally performed on a copy. The commands are:'/
     +'                  Translation           command:gtr  (gt)'/
     +'                  Rotation              command:gro  (gr)'/
     +'                  Scaling               command:gsc')
 1052 format(
     +'   It is also possible to perform reflection or inversion.'/
     +' These are performed on a copy unless otherwise requested.'/
     +' The commands for these operations are:'/
     +'                  Reflection            command:gre'/
     +'                  Inversion             command:gin'/
     +'   All these options accept input of  nf, nl  the range to'/
     +' be affected by the operation. In:gpr  one can also, if'/
     +' desired, specify a (different) range: ndf, ndl  to define'/
     +' the new origin and axes. Further Help available after the '/
     +' appropriate command has been given.')
 1006 format('                EDIT OPTIONS                 '//
     +' The available Edit options can be summarized as follows:'//
     +' EDIT OF ATOMS (and charges):'/
     +'     (0) RENUMBERING (permutational) options'/
     +' Move (with automatic renumbering)    command:emo  (e)'/
     +' Order (ct or an)                     command:eor  (eo)'/
     +' Permute                              command:epe  (ep)'/
     +'     (-) REMOVING options'/
     +' Remove & Renumber (to close gap)     command:err  (er) '/
     +' Remove (i.e. Delete)                 command:erm  (r)'/
     +'     (+) ADDING option'/
     +' Copy (a given range)                 command:ecp  (ec)'//
     +' EDIT OF Z-MATRIX'/
     +' Z-matrix (input, display, remove)    command:ezm  (ez)'//
     +'   When the command has been given a prompt will (as usual)'/
     +' appear. For all these commands one can then type:h for Help'/
     +' to get an explanation of the input and its effect.'/)
      end
*======================================================================
* BELOW THE BASIC ACCOUNTING ROUTINES 
* NACCNT,CACCNT,MACCNT,NOAR,FLS,NCHEK
*----------------------------------------------------------------------
      subroutine naccnt(nx,ifx,lfn,mxn,noa)
*  lfn= Least Free Number,   mxn= MaXimum Number used,
*  noa= Number Of Atoms entered. 
      logical ifx(nx)
      integer nx,lfn,mxn,noa,i
      do 1 i=nx,1,-1
         if(ifx(i)) go to 2
 1    continue
 2    mxn=i
      noa=0
      do 3 i=1,mxn
         if(.not.ifx(i)) go to 4
         noa=noa+1
 3    continue
 4    lfn=i
      do 5 i=lfn,mxn
         if(.not.ifx(i)) go to 5
         noa=noa+1
 5    continue
      if(lfn.gt.nx) lfn=0
      return
      end
*----------------------------------------------------------------------
      subroutine caccnt(nx,ifx,nq,nnq,mxn,nnh,ifz)
* This subroutine counts atoms with the same Z = nq(n) (atomic number)
* and the resulting order number is placed in nnq(n).
      logical ifx(nx),ifz
      integer nq(nx),nnq(nx),nx,mxn,nqn,i,j,nnh
* The number nnh is the number of non hydrogens.
      nnh=0
      do 1 i=1,nx
         if(nq(i).ne.1) nnh=nnh+1
 1    continue
      ifz=(nnh.gt.0)
      if(.not.ifz) return
      do 3 i=1,mxn
         if(.not.ifx(i)) go to 3
         nqn=1
         do 2 j=1,i-1
            if(ifx(j).and.nq(j).eq.nq(i)) nqn=nqn+1
 2       continue
         nnq(i)=nqn
 3    continue
      return
      end
*----------------------------------------------------------------------
      integer function maccnt(nx,imz)
* Finds the number of stored z-matrix rows of reference atoms
      integer nx, imz(nx), i
      maccnt=0
      do 1 i=1,nx
         if(imz(i).ne.0) maccnt=maccnt+1
 1    continue
      return
      end
*----------------------------------------------------------------------
      integer function noar(nx,ifx,nf,nl)
      logical ifx(nx)
      integer nx,nf,nl,i
      noar=0
      do 1 i=nf,nl
         if(ifx(i)) noar=noar+1
 1    continue
      return
      end
*-----------------------------------------------------------------------
      logical function fls(nx,n)
* This function is false if n is the range 1-nx, otherwise true
      integer nx,n
      fls=(n.lt.1.or.n.gt.nx)
* With e.g.:if(fls(nx,n)) return, one avoids n being out of range for array
      return
      end
*------------------------------------------------------------------------
      logical function arenot(nx,ifx,i,j,k)
* This function is true if one of the three atom numbers i,j,k is wrong
* or if one of them have not been entered.
      integer nx,i,j,k
      logical ifx(nx),fls
      arenot=.true.
      if(fls(nx,i)) return
      if(fls(nx,j)) return
      if(fls(nx,k)) return
      if(.not.ifx(i)) return
      if(.not.ifx(j)) return
      if(.not.ifx(k)) return
      arenot=.false.
      return
      end
*-----------------------------------------------------------------------
      subroutine nchek(str,ifc,nx,ifx,nn,io)
      logical ifx(nx),fls,ifh,ifc
      integer ivec(10),nx,nn,io,last,nint,n
      character str*(*)
 1    io=0
      if(.not.ifx(nn)) return
      write(6,1000) nn
 1000 format(
     +'   Atom number: n=',i3,'  already stored; to delete old'/
     +' and store new position, confirm this n, else change it:'/
     +' Enter: n                                       (n=0 to cancel)')
      str='#n?>'
      call readin(ifc,str,last,ifh)
      io=1
      if(last.eq.0.or.ifc) return
      call intstr(str(1:last),nint,ivec)
      if(nint.eq.0.or.nint.gt.1) return
      n=ivec(1)
      if(fls(nx,n)) return
      io=0
      if(nn.ne.n) io=1
      nn=n
      if(io.eq.1) go to 1
      return
      end
*===================================================================
      subroutine readin(ifc,str,last,ifh)
      character str*(*), savestr*80, filnm*20, dolln*2, dig*1
      character*25 consts(0:9)
      integer last, indsc, io, loop, inddl, lec, nocc
      logical ifh, ifc, ifil, noass, ifdln(0:9)
      common/dgt/dig(10)
      save savestr, indsc, ifil, ifdln, consts
      data indsc/0/, ifil/.false./, ifdln/10*.false./
 90   ifh=.false.
      ifc=.false.
      if(indsc.eq.0) then
* Write a prompt on the same line as the input is read
         if(.not.ifil) write(6,1000) str(1:4)//'> '
 1000    format(tr1,a6,$)
*
         savestr= '          '
         if(.not.ifil) read(5,1,iostat=io) savestr
         if(ifil) read(7,1,iostat=io) savestr
 1       format(80a)
         call strip(savestr,last)
         if(last.eq.0) then
            str='  '
            if(ifil.and.io.ne.0) then
               close (7)
               ifil=.false.
               str='file-end'
            endif
            return
         endif
      endif
      indsc=index(savestr,';')
      if(indsc.eq.0) then
         str=savestr
      else
         str= '  '
         if(indsc.gt.1) str=savestr(1:indsc-1)
         savestr=savestr(indsc+1:80)
      endif
      call strip(str,last)
 91   if(last.eq.0) return
      if(last.eq.4.and.str(1:4).eq.'quit') then
         write(6,*) ' EVCLID-run terminated.'
         stop
      endif
      if(index(str(1:last),'h').ne.0.and.last.le.2) ifh=.true.
      if(ifh) return
      if(last.eq.3) then
         ifc=lge(str(1:1),'a').and.lle(str(1:1),'z')
         ifc=ifc.and.lge(str(2:2),'a').and.lle(str(2:2),'z')
         ifc=ifc.and.lge(str(3:3),'a').and.lle(str(3:3),'z')
      endif
      if(str(1:2).eq.'<<') then
         if(last.eq.2) last=0
         if(last.eq.0) return
         filnm=str(3:last)
         call strip(filnm,last)
         inquire(file=filnm,exist=ifil)
         if(.not.ifil) then
            write(6,*) ' Command file:', filnm,'  not found.'
            return
         endif
         open(7,file=filnm,err=2,status='old')
         rewind (7)
         ifil=.true.
         indsc=0
         go to 90
      endif
      inddl=index(str,'$')
      if(inddl.eq.1) then
         if(str(3:3).eq.'=') then
            do 3 loop=0,9
               if(str(2:2).eq.dig(loop+1))  then
                  ifdln(loop)=.true.
                  consts(loop)=str(4:last)
                  write(6,*)
     +                 ' $'//dig(loop+1)//' assigned:'//consts(loop)
               endif
 3          continue
            str='assgmnt='
            last=0
            return
         elseif(last.eq.2.and.str(2:2).eq.'*') then
            noass=.true.
            do 4 loop=0,9
               if(ifdln(loop)) write(6,*)
     +              ' $'//dig(loop+1)//' assigned:'//consts(loop)
               noass=noass.and.(.not.ifdln(loop))
 4          continue
            if(noass) write(6,*) ' No assignments recorded.'
            str='assgmnt='
            last=0
            return
         else
            inddl=99
         endif
      endif
      if(inddl.gt.1) then
         do 5 loop=0,9
            dolln='$'//dig(loop+1)
            if(ifdln(loop)) then
               call strip(consts(loop),lec)
               call replac(str,dolln,consts(loop)(1:lec),nocc)
            else
               call replac(str,dolln,' ',nocc)
               if(nocc.ne.0) write(6,*)
     +              ' Warning, '//dolln//' not assigned.'
            endif
 5       continue
         if(index(str,'$').ne.0) then
            write(6,*) ' Warning $-sign occurs erroneously in line.'
            return
         endif
         call strip(str,last)
         go to 91
      endif
      return
 2    write(6,*)  ' Could not open command file:', filnm
      return
      end
*---------------------------------------------------------------------
      subroutine strip(str,last)
      character str*(*)
      integer last,first,len
      do 2 last= len(str), 1, -1
         if(str(last:last).ne.' ') go to 3
 2    continue
 3    if(last.eq.0) return
      do 4 first=1,last
         if(str(first:first).ne.' ') go to 5
 4    continue
 5    str=str(first:last)
      last=last-(first-1)
      return
      end
*---------------------------------------------------------------------
      subroutine replac(str,old,new,nocc)
      integer ls, lo, i, nocc
      character str*(*), old*(*), new*(*), STR1*100, STR2*100
      ls=len(str)
      lo=len(old)
      nocc=0
 1    continue
      STR1=' '//str
      STR2=' '//old//' '
      i=index(STR1(1:ls+1),STR2(1:1+lo+1))
      if(i.eq.0) then
         STR2(1:1)='-'
         i=index(STR1(1:ls+1),STR2(1:1+lo+1))
      endif
      if(i.eq.0) return
      if(i.eq.1) then
         str=new//str(lo+1:ls)
      else
         str=str(1:i-1)//new//str(i+lo:ls)
      endif
      nocc=nocc+1
      go to 1
      end
      
*===================================================================
* MAJOR GEOMETRIC PARAMETER INPUT SUBROUTINES (WITH ROOT-SEARCH)
*===================================================================
      subroutine igecom(str,ifc,nx,ifx,a,noa)
      implicit double precision(a-h,o-z)
      logical ifx(nx),fls,linear,ifh,ifc
      character prmp*79,str*(*),numbst*3
      dimension a(nx,3),fnum(10),v1(3),v2(3),v3(3),iv(6),rld(3),rr(3)
      common/dfu/lstdha,jt(3)/dfv/pa(3),ed(3,3),av(3,3),bv(3,3),cv(3,3)
      common/chk/jtr(3),im(3,3) /chp/gev(3) /pp/conv
      save /dfu/, /chk/, mem,nold, /dfv/, /chp/
      data mem/0/, nold/0/
*--------------------------------------------------------------------
* Initialize some case flags.
* Return here to prompt for more input: 
 90   nol=3
      linear=.false.
      nflg=0
      ori=1.
      if(noa.eq.0.or.noa.eq.1) nol=1
      if(noa.eq.2) nol=2
      prmp='Enter: ge i j ( k ( l))'//'    '//numbst(nol)
      call strip(prmp,last)
      prmp=' '//prmp(1:last)//' line'
      if(nol.gt.1) prmp=prmp(1:last+6)//'s'
      prmp(60:69)='(h = help)'
      prmp(75:79)='*ige*'
* Prompt for triple of geometric parameters.
*
      write(6,*) prmp
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*  Read in loop over number of required parameters +++++++++++++++++++
*
      do 19 mc=1,nol
*---------------------------------------------------------------------
 80   str='ige'//numbst(mc)
      call readin(ifc,str,last,ifh)
      if(last.eq.0.or.ifc) return
      if(ifh) then
         if(str(1:last).eq.'ha') then
            write(6,1012)
            go to 90
         elseif(str(1:last).eq.'hs') then
            write(6,1013)
            write(6,1014)
            go to 90
         elseif(str(1:last).eq.'hx') then
            write(6,1015)
            go to 90
         elseif(str(1:last).eq.'he') then
            write(6,1016)
            go to 90
         elseif(str(1:last).eq.'hr') then
            write(6,1018)
            go to 90
         elseif(str(1:last).eq.'hh') then
            write(6,1017)
            go to 90
         else
            write(6,1010)
            write(6,1011)
            go to 90
         endif
      endif
      call chtflp(str(1:last),nuno,fnum)
      ierr=0

      if(nuno.eq.0) then
         if(str(1:2).eq.'a ') nflg=1
         if(str(1:2).eq.'b ') nflg=2
         if(str(1:2).eq.'c ') nflg=3
         if(str(1:2).eq.'d ') nflg=4
         if(nflg.gt.0) go to 91
         write(6,*) ' Geometric parameter input wrong. Re-enter line.'
         go to 80
      elseif(nuno.ge.3) then
         ge=fnum(1)
         i=int(fnum(2))
         j=int(fnum(3))
         if(fnum(2)-float(i).ne.0..or.fnum(3)-float(j).ne.0.) ierr=1
         if(fls(nx,i).or.fls(nx,j)) ierr=1
         k=0
         l=0
         if(nuno.ge.4) then
            k=int(fnum(4))
            if(fnum(4)-float(k).ne.0..or.fls(nx,k)) ierr=1
         endif
         if(nuno.eq.5) then
            l=int(fnum(5))
            if(fnum(5)-float(l).ne.0..or.fls(nx,l)) ierr=1
         endif
         if(nuno.gt.5) then
            ierr=1
            write(6,*)' To many input numbers in this line !' 
         endif
      else
         ierr=1
         write(6,*)' To few input numbers in this line !' 
      endif
*--------------------------------------------------------------------
* If this point has been reached with ierr=0 then a line of input with
* correct ge,i,j,k,l format has been read in. If ierr=1 something was
* wrong with the line and it should be re-read, at 80.
*
      if(ierr.ne.0) then
         write(6,1030) mc
         go to 80
      endif
* 
* Handle the case of equal distance input
      if(str(1:2).eq.'r ') then
         n=int(ge)
         if(ge-float(n).ne.0..or.fls(nx,n)) then
            write(6,*) ' First number not acceptable atom number n.'
            go to 80
         elseif(ifx(n)) then
            go to 82
         elseif((.not.ifx(i)).or.(.not.ifx(j))) then
            go to 83
         endif
         if(nuno.eq.4) then
            if(.not.ifx(k)) go to 83
         elseif(nuno.eq.5) then
            if((.not.ifx(k)).or.(.not.ifx(l))) go to 83
         endif
         if(i.eq.j.or.i.eq.k.or.i.eq.l) go to 84
         if(j.eq.k.or.j.eq.l) go to 84
         if(k.ne.0.and.k.eq.l) go to 84
         ierr=0
         call eqdist(nx,a,nuno,i,j,k,l,rr,ierr)
         if(ierr.ne.0) go to 90
         nold=n
         go to 200
      endif
* Now determine syntactic correctness of input parameter. Find the
* number of the input atom n, and determine the type of the
* parameter jt(mc)=1,2,3,4, or 5. Find reference atom numbers i,j,k.
*---------------------------------------------------------------------
      n=0
*CASE: R,  DISTANCE JT=1 
* Check first if distance: r,i,j. Make: r,n,i and put jt=1.
      if(nuno.eq.3) then
         if(i.eq.j) go to 84
         if(ifx(j).and.(.not.ifx(i))) then
            n=i
            i=j
         endif
         if(ifx(i).and.(.not.ifx(j))) n=j
         if(noa.eq.0) then
            n=max(i,j)
            i=min(i,j)
         endif
         if(n.eq.0) then
            if(ifx(i).and.ifx(j)) go to 82
            go to 83
         endif
         j=0
         jt(mc)=1
*CASE: THT,  ANGLE JT=2 OR 4
* Check then if angle: tht,i,j,k. Make: tht,n,i,j and put jt=2
* or: tht,i,n,j and put jt=4.
      elseif(nuno.eq.4) then
         if(i.eq.j.or.i.eq.k.or.j.eq.k) go to 84
         if(.not.ifx(i)) n=i
         if(.not.ifx(j)) n=j
         if(.not.ifx(k)) n=k
         if(n.eq.0) go to 82
         if(n.eq.i) then
            if(ifx(j).and.ifx(k)) then
               i=j
               j=k
               jt(mc)=2
            else
               go to 83
            endif
         elseif(n.eq.k) then
            if(ifx(i).and.ifx(j)) then
               k=i
               i=j
               j=k
               jt(mc)=2
            else
               go to 83
            endif
         elseif(n.eq.j) then
            if(ifx(i).and.ifx(k)) then
               j=k
               jt(mc)=4
            else
               go to 83
            endif
         endif
         k=0
*CASE: PHI,  DIHEDRAL ANGLE JT=3 OR 5
* Check then if dihedral angle: phi,i,j,k,l. Make: phi,n,i,j,k and put 
* jt=3 or: phi,k,n,i,j and put jt=5.
      elseif(nuno.eq.5) then
         iv(1)=i
         iv(2)=j
         iv(3)=k
         iv(4)=l
         do 11 md=2,4
            do 11 me=1,md-1
               if(iv(md).eq.iv(me)) go to 84
 11      continue
         nnum=0
         do 12 md=1,4
            if(.not.ifx(iv(md))) then
               n=iv(md)
               nnum=nnum+1
            endif
 12      continue
         if(nnum.gt.1) go to 83
         if(n.eq.0) go to 82
*            We now know that all of i,j,k,l except that which is =n 
*            have been entered.
         if(n.eq.l) then
            l=k
            k=i
            i=j
            j=l
            ge=-ge
            jt(mc)=3
         endif
         if(n.eq.i) then
            i=j
            j=k
            k=l
            jt(mc)=3
         endif
         if(n.eq.j) then
            j=l
            l=k
            k=i
            i=l
            jt(mc)=5
         endif
         if(n.eq.k) then
            k=i
            i=j
            j=l
            ge=-ge
            jt(mc)=5
         endif
      endif
*ALL CASES HAVE BEEN ACCOUNTED 
* Now put dihedral angles and angles jt=4 on standard form.
      if(jt(mc).eq.3.or.jt(mc).eq.4) then
         if(i.gt.j) then
            l=i
            i=j
            j=l
            if(jt(mc).eq.3) ge=-ge
         endif
      elseif(jt(mc).eq.5) then
         if(j.gt.k) then
            l=j
            j=k
            k=l
            ge=-ge
         endif
      endif
* All parameters are now on a standard form so that if they are geome-
* trically equivalent they are identical (except the value ge, of 
* course, which is handeled below).
*--------------------------------------------------------------------
*
**** Handle distance R here. Case R   11111   
      if(jt(mc).eq.1) then
         pa(mc)=abs(ge)
         if(pa(mc).lt.1.D-3) go to 85
      endif
****  -----------------------------   11111
* Skip rest of do-loop here if noa less than or equal to 1. *EXIT*
      if(noa.le.1) go to 19
*-------------------------------------------------------------------
*  Now check if n is consistent, handle orientation sign
      if(mc.gt.1.and.n.ne.nold) write(6,1040) nold,n
      nold=n
      ori=ori*dsign(1.d0,ge)
* Put value of angles on standard form.
      if(jt(mc).ne.3.and.jt(mc).ne.5) ge=abs(ge)
      if( (jt(mc).ne.1) .and. (ge.le.-180. .or. ge.gt.180.) ) then
         ge=dmod(ge,360.d0)
         if(ge.le.-180.) ge=ge+360.
         if(ge.gt.+180.) ge=ge-360.
      endif
* Temporary print-out (now commented).
*      write(6,3000) jt(mc),ge,n,i,j,k,mc
* 3000 format(
*     +' Parm-type jt=',i1,'  ge=',f8.3,tr1,
*     +     '(  n=',i3,' i=',i3,' j=',i3,' k=',i3,' )  parm-no:',i1)
* Store standard form parameter: type, value and reference atoms:
      jtr(mc)=jt(mc)
      gev(mc)=ge
      im(mc,1)=i
      im(mc,2)=j
      im(mc,3)=k
* Store position of reference atom i in array av.
      do 13 md=1,3
         av(mc,md)=a(i,md)
 13   continue
*-------------------------------------------------------------------
**** Handle angle THT here. Case THT  22222   44444
      if(jt(mc).eq.2.or.jt(mc).eq.4) then
         call bij(nx,a,i,j,v1,d1,ierr)
         if(ierr.eq.1) go to 85
         pa(mc)=cos(conv*ge)
         if(jt(mc).eq.2) then
            do 14 md=1,3
               bv(mc,md)=v1(md)
 14         continue
         endif
         if(jt(mc).eq.4) then
            do 15 md=1,3
               bv(mc,md)=a(j,md)
 15         continue
         endif
*  Test if 0 or 180 and handle that case.
         if(abs(sin(conv*ge)).lt.1.D-4) then
            if(linear) go to 85
            write(6,*) ' When angle 0 or 180 two parameters suffice.'
            linear=.true.
            ml=mc
            mc2=3
            if(mc.eq.3) then
               write(6,*)' Your first parameter is ignored !!'
               mc2=1
               jtr(1)=0
            endif
            v2(1)=0.
            v2(2)=1.
            v2(3)=0.
            if(v1(1)**2+v1(3)**2.lt.1.D-5) v2(1)=1.
            call vpr(1,v1,v2,v3)
            call vpr(0,v1,v3,v2)
            jt(mc)=3
            jt(mc2)=3
            do 16 md=1,3
               bv(mc,md)=v2(md)
               bv(mc2,md)=v3(md)
               av(mc2,md)=a(i,md)
               ed(1,md)=a(i,md)
               ed(2,md)=0.1d0*pa(mc)*v1(md)
 16         continue
            pa(mc)=0.
            pa(mc2)=0.
            lstdha=0
            go to 81
         endif
      endif
***** ----------------------------   22222    44444 
* Skip rest of read-in do-loop from here if noa equal to 2.* EXIT*
      if(noa.eq.2) go to 19
*--------------------------------------------------------------------
**** Handle dihedral angle PHI here. Case PHI  33333   55555 
      if(jt(mc).eq.3.or.jt(mc).eq.5) then
         call refsys(nx,a,i,j,k,v1,v2,v3,sc,d1,d2,ierr)
         if(ierr.gt.0) go to 85
         if(jt(mc).eq.3) then
            pa(mc)=0.
            cc=cos(conv*ge)
            sc=sin(conv*ge)
            do 17 md=1,3
               bv(mc,md)=sc*v2(md)+cc*v3(md)
               ed(mc,md)=cc*v2(md)-sc*v3(md)
 17         continue
         elseif(jt(mc).eq.5) then
            pa(mc)=cos(conv*ge)
            call bij(nx,a,i,k,v2,d2,ierr)
            do 18 md=1,3
               bv(mc,md)=v1(md)
               cv(mc,md)=v2(md)
               ed(mc,md)=dsign(1.d0,ge)*v3(md)
 18         continue
         endif
      endif
**** ------------------------------------  33333    55555
* Leave read-in do-loop here if linear and mc at least two. Go 
* to check the two parameters before calling root-search.      *EXIT*
 81   if(linear.and.mc.ge.2) go to 91
*--------------------------------------------------------------------
      go to 19
*-------------------------Error-printouts for errors in loop 19.------
 82   write(6,1030) mc
      write(6,*) ' No new free atom number n given in this line.'
      go to 80
 83   write(6,1030) mc
      write(6,*) ' Not entered atom number referenced in this line.'
      go to 80
 84   write(6,1030) mc
      write(6,*) ' Atom numbers must be different in input line.'
      go to 80
 85   write(6,1030) mc
      write(6,*) ' Some geometric problem with this input line.'
      go to 80
*------------------------------------------------------------------
 19    continue
* Read-in do-loop number 19 ends here.**  ++++++++++++++++++++++++++
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* Jump-outs from loop land on next line.
 91   continue
*
* First 1 handle Case: SPECIAL OPTION ---------
      if(nflg.ne.0) then
         if(noa.le.2) go to 95
         n=-nflg
         if(n.gt.-4.and.mem.eq.0) go to 95
         if(n.eq.-3.and.mc.ne.1) go to 95
         if(n.eq.-3) go to 94
         if(n.eq.-1) then
*  If n=-1 change sign of last entered dihedral angle,recalculate.
            if(lstdha.eq.0) go to 95
            if(jtr(lstdha).ne.3.and.jtr(lstdha).ne.5) go to 95
            i=im(lstdha,1)
            j=im(lstdha,2)
            k=im(lstdha,3)
            ge=-gev(lstdha)
            gev(lstdha)=ge
            call refsys(nx,a,i,j,k,v1,v2,v3,sc,d1,d2,ierr)
            if(ierr.gt.0) go to 95
            cc=cos(conv*ge)
            sc=sin(conv*ge)
            do 21 md=1,3
               if(jt(lstdha).eq.3) then
                  bv(lstdha,md)=sc*v2(md)+cc*v3(md)
                  ed(lstdha,md)=cc*v2(md)-sc*v3(md)
               else
                  ed(lstdha,md)=-ed(lstdha,md)
               endif
 21         continue
            go to 94
         elseif(n.eq.-2) then
*  If n=-2 change sign of alternative sign-convention.
            if(lstdha.ne.0) go to 95
            ed(2,1)=-ed(2,1)
            ed(2,2)=-ed(2,2)
            ed(2,3)=-ed(2,3)
            go to 94
         elseif(n.eq.-4.and.mc.eq.3.and.nol.eq.3) then
            n=nold
            nflg=4
         else
            go to 95
         endif
      endif 
*
* Second 2 handle Case: NOA=0 OR 1 ------------
      if(noa.eq.0) then
         ifx(n)=.true.
         ifx(i)=.true.
         noa=noa+2
         do 22 md=1,3
            a(i,md)=0.
            a(n,md)=a(i,md)
 22      continue
         a(n,1)=a(i,1)+pa(1)
         write(6,1050)i,a(i,1),a(i,2),a(i,3)
         write(6,1050)n,a(n,1),a(n,2),a(n,3)
         go to 90
      elseif(noa.eq.1) then
         ifx(n)=.true.
         noa=noa+1
         do 23 md=1,3
            a(n,md)=a(i,md)
 23      continue
         a(n,1)=a(i,1)+pa(1)
         write(6,1050)n,a(n,1),a(n,2),a(n,3)
         go to 90
      endif
*
* Third 3 handle Case: NOA=2  -------------------
* Jump in here for special option d (nflg=4, treat as third atom).
      if(noa.eq.2.or.nflg.eq.4.and.(.not.linear)) then
         jtr(3)=0
         if(im(1,1).eq.im(2,1).and.im(1,2).eq.im(2,2).and.
     +        jt(1).eq.jt(2)) go to 97
         if(im(1,3).ne.0.or.im(2,3).ne.0) go to 96
         i=min0(im(1,1),im(2,1))
         if(im(1,2).ne.0) i=min0(i,im(1,2))
         if(im(2,2).ne.0) i=min0(i,im(2,2))
         j=max0(im(1,1),max0(im(1,2),max0(im(2,1),im(2,2))))
         call bij(nx,a,i,j,v1,d1,ierr)
         if(ierr.eq.1) go to 96
         v2(1)=1.
         v2(2)=0.
         v2(3)=0.
         if(v1(2)**2+v1(3)**2.lt.1.D-5) v2(2)=1.
         call vpr(1,v1,v2,v3)
         call vpr(0,v3,v1,v2)
         if(v2(2).lt.0.) v2(2)=-v2(2)
         do 24 md=1,3
            av(3,md)=a(i,md)
            bv(3,md)=v3(md)
            ed(3,md)=v2(md)
 24      continue
         jt(3)=3
         pa(3)=0.
         lstdha=3
         go to 94
      endif
*
* Fourth 4 now handle Case: LINEAR  -----------
      if(linear) then
         if(jtr(1).ne.0) jtr(3)=0
         ma=0
         do 25 md=1,3
            if(jt(md).ne.3) ma=md
 25      continue
         if(ma.ne.0) then
            if( ( jt(ma).eq.2 ) .and.
     +      (im(ma,1).eq.im(ml,1).and.im(ma,2).eq.im(ml,2)) ) go to 97
         endif
         go to 94
      endif
*
* Fifth 5 now handle Case: NORMAL -----------
* When 3 parameters have been entered check them:
      do 27 i=1,2
         do 271 j=i+1,3
            lch=iabs(jt(i)-jt(j)) + iabs(im(i,1)-im(j,1))+
     +           iabs(im(i,2)-im(j,2))+iabs(im(i,3)-im(j,3))
            if(lch.eq.0) go to 97
            if(jt(j).lt.2.or.jt(j).gt.3) go to 271
            if(jt(i).lt.2.or.jt(i).gt.3) go to 271
            if(pa(i).lt.1.D-11.and.pa(j).lt.1.D-11) then
               sc=0.
               do 26 k=1,3
                  sc=sc+bv(i,k)*bv(j,k)
 26            continue
               if(abs(sc).ge.0.999d0) then
                  write(6,*)
     +                 ' No solution; input specifies parallel planes !'
                  go to 96
               endif
            endif
 271     continue
 27   continue
*
* Find out if dihedral angle and if so the number of the last one lstdha.
      lstdha=0
      do 28 md=1,3
         if(jt(md).eq.3.or.jt(md).eq.5) lstdha=md
 28   continue
* Handle Case LSTDHA=0, no dihedral angle among the three parameters.
      if(lstdha.eq.0) then
         ncount=0
         do 30 i=1,3
            do 301 m=1,3
               if(im(m,i).eq.0) go to 301
               do 29 j=1,ncount
                  if(im(m,i).eq.iv(j)) go to 301
 29            continue
               ncount=ncount+1
               iv(ncount)=im(m,i)
 301        continue
 30      continue
         if(ncount.lt.3) then
            write(6,*) ' To few reference atoms.'
            go to 96
         endif
         do 31 i=1,ncount-1
            do 311 j=i+1,ncount
               if(iv(j).gt.iv(i)) go to 311
               isv=iv(j)
               iv(j)=iv(i)
               iv(i)=isv
 311        continue
 31      continue
         n1=iv(1)
         n2=iv(2)
         n3=iv(3)
         call refsys(nx,a,n1,n2,n3,v1,v2,v3,sc,d1,d2,ierr)
         if(ierr.gt.0) go to 96
         do 33 md=1,3
            ed(1,md)=(a(n1,md)+a(n2,md)+a(n3,md))/3.
            ed(2,md)=ori*v3(md)
 33      continue
      endif
*
 94   continue
*---------------------------------------------------------------------
*  Find root of f=0 in correct half-space, display it, store it.
      call rotsea(nx,a,n,ierr,rld,rr,f)
      if(ierr.gt.0) then
        write(6,*) ' No solution found!!!!'
        go to 90
      endif
      if(n.lt.0.and.ifx(nold)) then
         n=0
      else
         n=nold
      endif
 200  write(6,1050) n,rr(1),rr(2),rr(3)
      if(n.eq.0) then
        call nchek(str,ifc,nx,ifx,nold,ierr)
        if(ifc) return
        if(ierr.eq.1) go to 90
        write(6,1050) nold,rr(1),rr(2),rr(3)
      endif
      do 34 md=1,3
         a(nold,md)=rr(md)
 34   continue
      ifx(nold)=.true.
      noa=noa+1
      mem=mem+1
      go to 90
*
*  Go and prompt for new triple of parameters. Main part of
*  this routine ends. 
*---------- Error printouts for errors after do-loop.----------------
 95   write(6,*) ' Input wrong!!!! Special option not applicable.'
      go to 90
 96   write(6,1060) n
      go to 90
 97   write(6,*)' Input wrong, same parameter given twice !'
      write(6,1060) n
      go to 90
*********************below only format-text******               *end****
 1010 format('        GENERAL INFORMATION           '/
     +'   Notation in geometric parameter (ge i j k l) input:'/
     +' ge= a distance (r), an angle (tht), or a dihedral angle (phi).'/
     +'   Among the numbers: i j k l  one must find:'/
     +' n= the number of the atom you wish to specify;   '/
     +' i j k= numbers of atoms that already have been entered.'/
     +'   Thus, enter:                                   '/
     +' r n i        if distance=r from n to i is known,'/
     +' tht n i j    if angle=tht at i between i-n and i-j is known'/
     +' phi n i j k  if angle=phi at line i-j between planes',
     +' n-i-j and i-j-k is known.'/
     +' tht i n j    if angle=tht at n between n-i and n-j is known'/
     +' phi i n j k  if angle=phi at line n-j between planes',
     +' i-n-j and n-j-k is known.'/
     +'  (The relative order of the numbers is not important. E.g. r i',
     +' n also works.)'/
     +'   phi  are signed angles (-180 to 180 degrees) but',
     +' tht  un-signed (0 to 180).'/
     +'   Normally 3 lines (parameters) of input are needed.')
 1011 format('        MENU FOR FURTHER HELP           '/
     +' Type:hh   for example of input with explanations.'/
     +' Type:ha   for help about Starting when fewer than three atoms.'/
     +' Type:hs   for help about Sign conventions.     '/
     +' Type:hx   for help about Extra features and special cases.'/
     +' Type:he   for help about Error codes.'/
     +' Type:hr   for info on equal (but unknown) distances input [r n',
     +' i j (k (l))].')
 1012 format('    THE FIRST THREE ATOMS                    '/
     +'   To specify the positions of the first two atoms one line of'/
     +' input is enough: r m n (no reference atom needed in this case)'/
     +' This places the atom with the lower number at the origin '/
     +' (0,0,0) and the other r along the x-axis (r,0,0).'/
     +'   Should one atom have been entered (e.g. using coordinates)'/
     +' it must be referenced i.e. the input must be: r n i. Atom '/
     +' number n is then placed at (x+r,y,z) if atom i is at (x,y,z).'/
     +'   For the third atom two input lines are needed specifying two'/
     +' distances and/or angles (of type tht). It is placed in the xy-'/
     +' plane with positive y-coordinate if the first two are along'/
     +' the x-axis as a result of input described above. Else it is'/
     +' placed in a suitably choosen plane containing the first two'/
     +' atoms. Input for the third atom should thus be two parameters'/
     +' of the following types:  r n i  or   tht n i j  or   tht i n j')
 1013 format('    SIGN-CONVENTIONS                   '/
     +'   The specification of a position with three geometric para-'/
     +' meters is often not unique. Unless some sign convention is'/
     +' used there will normally be two solutions. The sign-conven-'/
     +' tion for dihedral (or torsion) angles is one of these. If '/
     +' there is no dihedral angle in the specification an alterna-'/
     +' tive described further below is used.               '/
     +'   The dihedral angle phi(i,j,k,l) increases in the mathemati-'/
     +' cally positive sense from 0 to 360 around the directed line-'/
     +' segment from j to k (j-k) from the plane i-j-k to j-k-l.'/
     +'   Equivalently one can say that phi is between 0 and 180 when'/
     +' i is on the negative side of the oriented triangle j-k-l or,'/
     +' alternatively, when l is on the positive side of the oriented'/
     +' triangle i-j-k. Orientation given by the right hand rule.')
 1014 format(
     +'   If no dihedral angle is among the 3 parameters, normally two'/
     +' different solutions (positions) are compatible with them. In '/
     +' this case the (numerically first) three reference atoms used '/
     +' in the specification are considered to define an oriented tri-'/
     +' angle (oriented by the numerical order) and a solution on the'/
     +' positive side is choosen. By giving a negative sign to one of'/
     +' the three parameters the other solution is obtained.'/
     +'   Note: neither convention guarantees uniqueness in all cases.')
 1015 format('    EXTRA FEATURES               '/
     +'   Since the sign-conventions are difficult to get right every'/
     +' time the following possibility exists: After a specification'/
     +' has successfully produced a solution one can type:'/
     +'  a  to find a solution with the opposite sign of the last'/
     +'     entered dihedral angle in specification just given. Or'/
     +'  b  to find new solution with opposite sign-convention for'/
     +'     the case of no dihedral angle. Or, finally,'/
     +'  c  to make the program look for a new solution to your'/
     +'     preceding specification (test for non-uniqueness).'/
     +'   If the shape is linear it is not possible to give a third'/
     +' parameter. In that case it is possible to just type:'/
     +'  d  instead of the third line of input. A solution is then'/
     +'     found according to the conventions for the two atom case.'/
     +'   Note that if an angle tht=0 or =180 the atom is specified'/
     +' to lie along a line (not a surface as otherwise) and only one'/
     +' further parameter is needed to specify it. The program can '/
     +' handle this case and will only prompt for two parameters.')
 1016 format(
     +'          ERROR CODES (for ige root search)'/
     +' -1 means that solution found by converged amoeba routine.'/
     +'      Results in no message.'/
     +' -2 means that amoeba did not converge in 500 iterations.'/
     +'      Results in message: Possibly reduced accuracy.'/
     +' 1  means that solution approached branchpoint of surface.'/
     +' 2  means that determinant of linear equations was zero.'/
     +' 3  means that sum of squares function was not zero at minimum.'/
     +' 10 means that dihedral angle had wrong sign.'/
     +' 11 means that alternative sign convention was violated.'/
     +' 20 means that new root was within 1e-10 of old root.'/
     +' 21 means that new amoeba root within 1e-4 of old root.'/
     +'      Three errors with code 3 result in message: Suspect non-'/
     +'    zero minimum; solution might not exist.')
 1017 format('       EXAMPLE OF INPUT TO ige'/
     +'   The following input lines result in the corners of a tetra-'/
     +' hedron with side 1.5 (there must be no previous input so the'/
     +' prompt asks for 1 line):'/
     +'  1st prompt, 1st line:  1.5  2  1'/
     +' Now you get two points stored and their coordinates are shown.'/
     +' The prompt is now for two lines and possible input is:'/
     +'  2nd prompt, 1st line:  1.5  3  1'/
     +'  2nd prompt, 2nd line:  60.0  1  3  2'/
     +' This results in storage of atom 3 at a distance 1.5 from atom'/
     +' 1 and the angle at 3 between 3-1 and 3-2 is 60.0 degrees. The'/
     +' resulting Cartesian coordinates are shown on the screen. Now'/
     +' the prompt will be for 3 lines of input and these might be:'/
     +'  3rd prompt, 1st line: 60.0  4  1  2'/
     +'  3rd prompt, 2nd line: 60.0  4  1  3'/
     +'  3rd prompt, 3rd line: 1.5  1  4'/
     +' The calculated Cartesian coordinates are again shown, and the'/
     +' regular tetrahedron is finished. Note that the program keeps'/
     +' track of which atoms are entered and which atoms are new.')
 1018 format('   EQUAL (UNKNOWN) DISTANCES INPUT'/
     +'   The position of an atom can be specified by demanding that'/
     +' its distances to other (entered) atoms are equal. This is'/
     +' done by typing: r  n  i j (k (l))'/
     +' i.e. you put the letter r  first in the input line, then the'/
     +' number, n, of the atom to be specified, and then the numbers'/
     +' of the entered atoms to which the distances must be equal.'/
     +'   If two reference atoms, i and j, are given the atom n will'/
     +' be placed at the mid-point between them.'/
     +'   If three reference atoms (i j k) are given the atom n will'/
     +' be placed in the same plane as these three at the center of'/
     +' the circumscribed circle of the triangle defined by:  i j k.'/
     +' These may not thus not lie on the same line.'/
     +'   If four reference atoms are given the point at the center'/
     +' of the circumscribed sphere of the tetrahedron defined by:'/
     +' i j k l  is specified. These may thus not all lie in a plane.')
 1030   format(' Something wrong with this input line no:',i1,
     +  ' try again!')
 1040 format(' Warning, you changed input atom no from:',i3,' to:',i3)
 1050 format(1x,i3,3x,3(2x,f15.9))
 1060 format(' Input wrong for atom number:',i3,' ,try again !!!')
      end
*=====================================================================
      subroutine rotsea(nx,a,nn,ierr,rld,r,f)
      implicit double precision(a-h,o-z)
*  A root of f=0 is sought in the desired half-space(s). If zrgrh
*  finds a zero of the gradient in wrong half-space a new initial
*  r is tried up to 16 times. nn=-3 means find another root!
      dimension aa(3),ee(3),r(3),rld(3),rm(3),ierrvc(0:16),a(nx,3)
      common/dfu/lstdha,jt(3)/dfv/pa(3),ed(3,3),av(3,3),bv(3,3),cv(3,3)
      common/chk/jtr(3),im(3,3) /chp/gev(3)
      common/pp/conv
      save /dfu/, /chk/, /dfv/, /chp/
      ic=0
* First determine a start guess in the correct half-space.
      do 11 i=1,3
         if(lstdha.eq.0) then
            aa(i)=ed(1,i)
            ee(i)=ed(2,i)
         else
            aa(i)=av(lstdha,i)
            ee(i)=ed(lstdha,i)
         endif
 11   continue
* Determine scale for root search.
      if(jtr(1).gt.1)fac=gprm(nx,a,im(1,1),im(1,2),0,0)
      if(jtr(1).eq.1) then
         if(jtr(2).gt.1)fac=gprm(nx,a,im(2,1),im(2,2),0,0)
         if(jtr(2).eq.1)fac=gprm(nx,a,im(1,1),im(2,1),0,0)
      endif
      fac=fac*2.03461947208
*
      do 12 i=1,3
         r(i)=aa(i)+fac*ee(i)
 12      rm(i)=r(i)
      nerrcn=0
*
* Now seek a minimum of the sum of the squares of the surface functions
* which is also a zero (i.e. a root). 
* Return to this point with new starting guesses if necessary.
 90   ierr=0
      call zrgrh(r,ierr)
 91   continue
      ierrvc(ic)=ierr
      if(ierr.eq.3) nerrcn=nerrcn+1
      if(nerrcn.eq.3) go to 93
      if(ierr.gt.0) then
* Make a new start-vector r here. Then go back to 90
         ic=ic+1
         if(ic.gt.8) fac=2.58324813*fac
         if(ic.gt.16) go to 92
         nc=0
         do 13 i=-1,1,2
            do 13 j=-1,1,2
               do 13 k=-1,1,2
                  if(nc.eq.mod(ic,8)) go to 14
                  nc=nc+1
 13      continue
 14      r(1)=rm(1)-i*fac
         r(2)=rm(2)-j*fac
         r(3)=rm(3)-k*fac
         go to 90
      endif
*
* Now check that signs of dihedral angles are correct:
      if(lstdha.gt.0) then
         do 18 i=1,lstdha
            if(jt(i).ne.3.and.jt(i).ne.5) go to 18
            coll=0.
            do 17 j=1,3
 17            coll=coll+(r(j)-av(i,j))*ed(i,j)
            if(coll.lt.-1.D-10) ierr=10
 18      continue
* If no dihedral angle check the alternative sign-convention.
      elseif(lstdha.eq.0) then
         coll=0.
         do 19 i=1,3
            coll=coll+(r(i)-aa(i))*ee(i)
 19      continue
         if(coll.lt.-1.D-10) ierr=11
      endif
* If nn=-3 then the task is to find a root different from the previous.
      if(nn.eq.-3) then
         dis=abs(r(1)-rld(1))+abs(r(2)-rld(2))+abs(r(3)-rld(3))
         if(ierr.le.-1.and.dis.lt.1.D-4) ierr=21
         if(dis.le.1.D-10) ierr=20
      endif
      if(ierr.gt.0) go to 91
*
* If success this point is reached. Save result and return.
      do 20 i=1,3
 20      rld(i)=r(i)
         if(ierr.eq.-2) write(6,*) ' Possibly reduced accuracy.'
      return
*--------------------------------------------------------------------
 92   ierr=1
      write(6,*) ' The following errors occured:'
      write(6,*) (ierrvc(i), i=0,ic-1)
      return
 93   ierr=1
      write(6,*) ' Suspect non-zero minimum; solution might not exist.'
      return
      end
*====================================================================
      subroutine zrgrh(r,ierr)
      implicit double precision(a-h,o-z)
      dimension r(3),f(3),gr(3,3),gv(3),hs(3,3),indx(3)
      dimension p(4,3),y(4),rold(3)
*--------------------------------------------------------------
      ita=0
      itl=0
      ith=1
      fold=1.D35
 1    call fgrhs(r,fun,gv,hs,ierr)
      if(ierr.gt.0) go to 100
      if(fun.gt.fold) go to 6 
      fold=fun
      call ludcmp(hs,indx,ierr)
      if(ierr.gt.0) go to 200
      call lubksb(hs,indx,gv)
      do 2 i=1,3
         rold(i)=r(i)
 2       r(i)=r(i)-gv(i)
      if(fun.lt.1.D-22) go to 90
      if(fun.lt.1.D-8) go to 4
      if(ith.gt.17) go to 6
      ith=ith+1
      go to 1
*------------------------------------------------------------------
 4    itl=itl+1
      call fncgr(r,fun,f,gr,ierr)
      if(ierr.gt.0) go to 100
      if(fun.gt.fold) then
         ftol=1.D-38
         go to 8
      endif
      fold=fun
      call ludcmp(gr,indx,ierr)
      if(ierr.gt.0) go to 200
      call lubksb(gr,indx,f)
      do 5 i=1,3
         rold(i)=r(i)
 5       r(i)=r(i)-f(i)
      if(fun.lt.1.D-22) go to 90
      go to 4
*--------------------------------------------------------------------
 6    ftol=1.D-10
 8    continue
      do 9 i=1,4
         do 9 j=1,3
 9          p(i,j)=rold(j)
      p(2,1)=p(2,1)+0.1234901
      p(3,2)=p(3,2)+0.0967324
      p(4,3)=p(4,3)+0.0828011
      do 11 i=1,4
         do 10 j=1,3
 10         r(j)=p(i,j)
 11      y(i)=func(r,ierr)
      call amoeba(p,y,ftol,ita,ierr)
      if(ierr.gt.0) go to 100
      do 12 i=1,3
 12      r(i)=p(1,i)
      if(ftol.gt.1.D-35) go to 4
      ierr=-1
      if(ita.ge.500) ierr=-2
*--------------------------------------------------------------------
 90   if(func(r,ier).gt.1.D-22) go to 300
      if(ier.gt.0) go to 100
*      write(6,1000) ith,itl,ita
* 1000 format('Convergence in zrgrh after ith=',i2,' itl=',i2,' ita=',i3)
      return
*
 100  ierr=1
      go to 500
 200  ierr=2
      go to 500
 300  ierr=3
 500  continue
*      write(6,2000) ierr,ith,itl,ita
* 2000 format(
*     +'Error type:',i1,' in zrgrh after ith=',i2,' itl=',i2,' ita=',i3)
      return
*
      end
*====================================================================
      function func(r,ierr)
      implicit double precision(a-h,o-z)
      dimension r(3),dr(3),er(3),fr(3)
      common/dfu/lstdha,jt(3)/dfv/pa(3),ed(3,3),av(3,3),bv(3,3),cv(3,3)
      save /dfu/, /dfv/
*
      ierr=1
      func=0.
      do 4 mc=1,3
         do 10 i=1,3
 10         dr(i)=r(i)-av(mc,i)
         d2=scp(dr,dr)   
         dd=dsqrt(d2)
         if(dd.lt.1.D-4) return
*--------------------------------------------------------------------
         if(jt(mc).eq.1) then
            fm=dd-pa(mc)
*--------------------------------------------------------------------
         elseif(jt(mc).eq.2) then
            do 12 i=1,3
 12            er(i)=bv(mc,i)
            fm=scp(dr,er)-dd*pa(mc)
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.3) then
            do 14 i=1,3
 14            er(i)=bv(mc,i)
            fm=scp(dr,er)
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.4) then
            do 16 i=1,3
 16            er(i)=r(i)-bv(mc,i)
            de=dsqrt(scp(er,er))
            if(de.lt.1.D-4) return
            fm=scp(dr,er)-pa(mc)*dd*de
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.5) then
            do 18 i=1,3
               er(i)=bv(mc,i)
 18            fr(i)=cv(mc,i)
            bb=scp(er,fr)
            sc1=scp(dr,er)
            sc2=scp(dr,fr)
            dd1=dsqrt(abs(d2-sc1*sc1))
            dd2=dsqrt(abs(d2-sc2*sc2))
            if(dd1.lt.1.D-4.or.dd2.lt.1.D-4) return
            fm=d2*bb-sc1*sc2-pa(mc)*dd1*dd2
         endif
*------------------------------------------------------------------------
         func=func+fm*fm
 4    continue
      ierr=0
      return
      end
*====================================================================
      subroutine fncgr(r,fun,f,gr,ierr)
      implicit double precision(a-h,o-z)
      dimension r(3),f(3),gr(3,3),dr(3),er(3),fr(3)
      common/dfu/lstdha,jt(3)/dfv/pa(3),ed(3,3),av(3,3),bv(3,3),cv(3,3)
      save /dfu/, /dfv/
*
      ierr=1
*  Do-loop 4 sums over contributions from the three parameters.
*  The values at r of the three functions are returned in f, the
*  matrix constructed from their gradients in gr.
      do 4 mc=1,3
         do 10 i=1,3
 10         dr(i)=r(i)-av(mc,i)
         d2=scp(dr,dr)   
         dd=dsqrt(d2)
         if(dd.lt.1.D-4) return
*--------------------------------------------------------------------
         if(jt(mc).eq.1) then
            f(mc)=dd-pa(mc)
            do 11 i=1,3
 11            gr(mc,i)=dr(i)/dd
*--------------------------------------------------------------------
         elseif(jt(mc).eq.2) then
            do 12 i=1,3
 12            er(i)=bv(mc,i)
            f(mc)=scp(dr,er)-dd*pa(mc)
            do 13 i=1,3
 13            gr(mc,i)=er(i)-dr(i)*pa(mc)/dd
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.3) then
            do 14 i=1,3
 14            er(i)=bv(mc,i)
            f(mc)=scp(dr,er)
            do 15 i=1,3
 15            gr(mc,i)=er(i)
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.4) then
            do 16 i=1,3
 16            er(i)=r(i)-bv(mc,i)
            de=dsqrt(scp(er,er))
            if(de.lt.1.D-4) return
            f(mc)=scp(dr,er)-pa(mc)*dd*de
            do 17 i=1,3
 17            gr(mc,i)=dr(i)+er(i) - pa(mc) * (dr(i)*de/dd+er(i)*dd/de)
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.5) then
            do 18 i=1,3
               er(i)=bv(mc,i)
 18            fr(i)=cv(mc,i)
            bb=scp(er,fr)
            sc1=scp(dr,er)
            sc2=scp(dr,fr)
            dd1=dsqrt(abs(d2-sc1*sc1))
            dd2=dsqrt(abs(d2-sc2*sc2))
            if(dd1.lt.1.D-4.or.dd2.lt.1.D-4) return
            f(mc)=d2*bb-sc1*sc2-pa(mc)*dd1*dd2
            do 19 i=1,3
               gr(mc,i)= 2.*bb*dr(i) - er(i)*sc2 - fr(i)*sc1 -
     -         ((dr(i)-er(i)*sc1)*dd2/dd1 + (dr(i)-fr(i)*sc2)*dd1/dd2)
 19         continue
         endif
*------------------------------------------------------------------------
 4    continue
      fun=scp(f,f)
      ierr=0
      return
      end
*====================================================================
      subroutine fgrhs(r,fun,gv,hs,ierr)
      implicit double precision(a-h,o-z)
      dimension r(3),gr(3),gv(3),hs(3,3)
      dimension dr(3),er(3),fr(3),der(3),dfr(3),h(3,3)
      common/dfu/lstdha,jt(3)/dfv/pa(3),ed(3,3),av(3,3),bv(3,3),cv(3,3)
      save /dfu/, /dfv/
      dlt(i,j)=float(max0(1-iabs(i-j),0))
* Statement function defines the Kronecker delta.
*
      ierr=1
      fun=0.
      do 9 i=1,3
         gv(i)=0.
         do 9 j=1,i
            hs(i,j)=0.
 9    continue
*
*  Do-loop 22 sums over contributions from the three parameters.
*  The gradient and the hessian of the sum of the squares of the
*  three surface functions are returned in gv and hs resp.
*
      do 22 mc=1,3
         do 10 i=1,3
 10         dr(i)= r(i) - av(mc,i)
         dsq= scp(dr,dr)
         dd= dsqrt(dsq)
         if(dd.lt.1.D-4) return
*--------------------------------------------------------------------
         if(jt(mc).eq.1) then
            f= dd - pa(mc)
            do 11 i=1,3
               gr(i)= dr(i)/dd
               do 11 j=1,i
 11               h(i,j)=  ( dlt(i,j) - dr(i) * dr(j)/dsq ) /dd
*--------------------------------------------------------------------
         elseif(jt(mc).eq.2) then
            do 12 i=1,3
 12            er(i)= bv(mc,i)
            f= scp(dr,er) - dd*pa(mc)
            do 13 i=1,3
               gr(i)= er(i) - dr(i)*pa(mc)/dd
               do 13 j=1,i
 13               h(i,j)= - pa(mc) * ( dlt(i,j) - dr(i)*dr(j)/dsq ) /dd
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.3) then
            do 14 i=1,3
 14            er(i)=bv(mc,i)
            f= scp(dr,er)
            do 15 i=1,3
               gr(i)= er(i)
               do 15 j=1,i
 15               h(i,j)= 0.
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.4) then
            do 16 i=1,3
 16            er(i)= r(i) - bv(mc,i)
            esq= scp(er,er)
            de= dsqrt(esq)
            if(de.lt.1.D-4) return
            ddd=dd*de
            d1=dd/de
            d2=de/dd
            f= scp(dr,er) - pa(mc)*ddd
            do 17 i=1,3
               gr(i)= dr(i) + er(i) - pa(mc) * ( dr(i)*d2 + er(i)*d1 )
               do 17 j=1,i
 17               h(i,j)=  dlt(i,j) * ( 2.- pa(mc) * (d2+d1) ) +
     +            pa(mc) * ( dr(i)*dr(j)*d2/dsq + er(i)*er(j)*d1/esq
     -            - ( dr(i)*er(j) + dr(j)*er(i) )/ddd )
*-----------------------------------------------------------------------
         elseif(jt(mc).eq.5) then
            do 18 i=1,3
               er(i)=bv(mc,i)
 18            fr(i)=cv(mc,i)
            bb=scp(er,fr)
            sc1=scp(dr,er)
            sc2=scp(dr,fr)
            d1sq=dsq-sc1*sc1
            d2sq=dsq-sc2*sc2
            dd1=dsqrt(abs(d1sq))
            dd2=dsqrt(abs(d2sq))
            if(dd1.lt.1.D-4.or.dd2.lt.1.D-4) return
            ddd=dd1*dd2
            d12=dd1/dd2
            d21=dd2/dd1
            f= dsq*bb - sc1*sc2 - pa(mc)*ddd
            do 19 i=1,3
               der(i)=dr(i)-er(i)*sc1
 19            dfr(i)=dr(i)-fr(i)*sc2
            do 20 i=1,3
               gr(i)= 2.*bb*dr(i) - er(i)*sc2 - fr(i)*sc1 -
     -         pa(mc) * ( der(i)*d21 + dfr(i)*d12 )
               do 20 j=1,i
 20               h(i,j)= dlt(i,j) * ( 2.*bb-pa(mc)*(d21+d12) )-
     -            ( er(i)*fr(j) + er(j)*fr(i) ) + 
     +            pa(mc)* (  er(i)*er(j)*d21 + fr(i)*fr(j)*d12 +
     +            der(i)*der(j)*d21/d1sq + dfr(i)*dfr(j)*d12/d2sq -
     -            ( der(i)*dfr(j) + der(j)*dfr(i) )/ddd )
         endif
*------------------------------------------------------------------------
         fun=fun + f*f
         do 21 i=1,3
            gv(i)= gv(i) + 2. * f * gr(i)
            do 21 j=1,i
               hs(i,j)= hs(i,j) + 2. * ( gr(i) * gr(j) + f * h(i,j) )
               hs(j,i)= hs(i,j)
 21      continue
 22   continue
      ierr=0
      return
      end
*------------------------------------------------------------------------
      subroutine eqdist(nx,a,nuno,i,j,k,l,rr,ierr)
      double precision a(nx,3),rr(3),b1(3),b2(3),b3(3),bm(3,3),scp
      double precision v1(3),v2(3),v3(3),sc,d1,d2,vv1(3),vv2(3)
      integer i,j,k,l,ierr,nuno,ii,indx(3)
      ierr=0
      do 1 ii=1,3
         v1(ii)=0.5*(a(i,ii)+a(j,ii))
         rr(ii)=v1(ii)
 1    continue
      if(nuno.eq.3) return
      do 2 ii=1,3
         v2(ii)=0.5*(a(i,ii)+a(k,ii))
         if(nuno.eq.4) v3(ii)=a(i,ii)
         if(nuno.eq.5) v3(ii)=0.5*(a(i,ii)+a(l,ii))
 2    continue
      call bij(nx,a,i,j,b1,d1,ierr)
      if(ierr.ne.0) return
      call bij(nx,a,i,k,b2,d1,ierr)
      if(ierr.ne.0) return
      if(nuno.eq.5) then
         call bij(nx,a,i,l,b3,d1,ierr)
      elseif(nuno.eq.4) then
         call refsys(nx,a,i,j,k,vv1,vv2,b3,sc,d1,d2,ierr)
      endif         
      if(ierr.ne.0) return
      do 4 ii=1,3
         bm(1,ii)=b1(ii)
         bm(2,ii)=b2(ii)
         bm(3,ii)=b3(ii)
 4    continue
      call ludcmp(bm,indx,ierr)
      if(ierr.ne.0) return
      rr(1)=scp(v1,b1)
      rr(2)=scp(v2,b2)
      rr(3)=scp(v3,b3)
      call lubksb(bm,indx,rr)
      return
      end
