********************************************************************************
*                                                                              * 
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### rinputd ###
c
c     Handles the input to 2DHF.
c
      subroutine rinputd(iform_t,ni_t,mu_t,no_t,nons_t)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      character*8 clabel,char8
      character*8 labellc(30),cmethod(6),pab(3)
      character*8 sigma,sigmag,sigmau,pi,delta,phi,space,angstrom,
     &     endl,endc
      character*20 psiinp,coulinp,exchinp,psiout,coulout,exchout
      include 'commons8.inc'
      dimension id(600)
      data nlabels/28/
      data labellc /'title','method','nuclei','config','initial',
     &     'break','multipol','grid','cutoff','homo',
     &     'scf','sor','debug','help','stop','interp',
     &     'subgrid','fix','xalpha','omega','order',
     &     'plot','fermi','lagra','gauss','damp','canon',
     &     'fefield',2*'        '/
c
      data angstrom/'angstrom'/
      data cmethod/'hf','HF','oed','OED','hfs','HFS'/,
     &     endl/'end'/,endc/'END'/
      data pab/'u','g',' '/
      data sigma/'sigma'/,pi/'pi'/,delta/'delta'/,phi/'phi'/
      data sigmag/'sigma g'/,sigmau/'sigma u'/
      data space/'    '/
      data inpiexit/-99999/

c  no=isum=norb
c  iexd=maxorb
c  bond=orbsym

c  initialization of some variables

      ihomon=0
      ibreak=0
      iinterp=0
     
      cutorb =1.d-6
      cutcoul=1.d-6
      cutexch=1.d-6
      itail  =13    

c     diver  -  if demax(1) in proce is greater than diver
c               divergence in the scf process is detected and
c               the program stops

      diver = 1.d8

c     initial default value
      
      inclorb = 0

c     default canonicalization coefficients

      icanon=1
      
c     multipol default value
      
      impole=8
c     default number of grids and the default ordering of mesh points
      ngrids=1
      iorder(ngrids)=2

c     scf default values
      
      maxscf  = 1000
      nobckup =  100
      ienterm =   10
      inoterm =   10
      iprtlev =    2
      facmul  =    1.15d0

c     fix default values

      iexlorb  = 0
      iexlcoul = 0
      iexlexp  = 2
   
      exlorb  = 0.d0
      exlcoul = 0.d0
      exlexp  = 2.d0

c     xalpha defaul value

      alphaf = 0.0d0
      sltthr = 0.0d0
      
c     sor default value
      
      maxsor2 =10
      maxsor1 = 1
      lsor    = 1
      ipoiss  = 1

c     fermi default value (point nuclei)

      ifermi=0
      
c     plot default value (no data exported for making plots)

      iplot=0
      
c     additional relaxation of the least relaxed orbital is switched off
c     facrel is no longer used	

      facrel=-10.d0

c     search for optimal value of the overrelaxation parameter is
c     switched off
      
      search=1.d-15
      omincr=0.0d0
      
      sflagra=1.d0

      ifefield=0
      ffield =0.d0
      fgrad  =0.d0
      zcutoff=1.d5
      
      cmass =0.d0
      cmass1=0.d0

c     default order of interpolation polynomials when changing grids 
c     for orbitals and potentials

      iord_nu_orb=5
      iord_nu_coul=5
      iord_nu_exch=5
      
      iord_mu_orb=5
      iord_mu_coul=5
      iord_mu_exch=5

      incrni=10
      incrmu=10

c     set maximum allowed unit number
      call setmaxunit()

      write(*,*) '  ...  start of input data  ...'
      write(*,*)'        '

 5    continue
      
      call card
      call inpa(clabel)
      
      do i=1,nlabels
         if (clabel.eq.labellc(i)) go to 12
      enddo
      
      write(iout6,1100) clabel
 1100 format(/10x,'******  ',a8,' not yet in the menu ******'//)
      write(iout6,1102) (labellc(i),i=1,nlabels)
 1102 format(/10x,'try one of the following :'//
     &	     10(10x,a8/))
      stop
      
 12   go to ( 20, 40, 60, 80,100,
     &     120,140,160,180,200,
     &     220,240,260,280,300,
     &     320,162,222,224,244,
     &     242,340,360,370,362,
     &     380,390,392), i
      stop

c     label: title

 20   continue
      read(iinp5,1104,err=1500) header
      write(iout6,1104) header
 1104 format(80a1)
      go to 5

c     label: method
      
 40   continue
c     call card
      call inpa(char8)
      if (char8.eq.cmethod(1).or.char8.eq.cmethod(2)) imethod=1
      if (char8.eq.cmethod(3).or.char8.eq.cmethod(4)) imethod=2
      if (char8.eq.cmethod(5).or.char8.eq.cmethod(6)) imethod=3
     
c     prepare for HFS calculations
      
      if (imethod.eq.3) then
         alphaf=0.70d0
         sltthr=precis
	 exlexp=1.d0
	 islat =1
      endif
      
      if (char8.eq.cmethod(1).or.char8.eq.cmethod(2).or.
     &    char8.eq.cmethod(3).or.char8.eq.cmethod(4).or.
     &    char8.eq.cmethod(5).or.char8.eq.cmethod(6)) go to 5
      write(iout6,*) 'method must be HF, HFS or OED for a time being'
      stop 'rinputd'

c     label: nuclei

 60   continue

      call inpf(z1)
      call inpf(z2)
      call inpf(r)
      call inpa(char8)

c     conversion factor due to Cohen and Taylor (1986)
c     The 1986 Adjustment of the Fundamental Physical Constants

      if (char8.eq.angstrom) r=r/bohr2ang
      r2=r/2.0d0

      if ((abs(z1)+abs(z2)).eq.0.d0.or.r.eq.0.d0) then
         write(iout6,*) 'Error: incomplete input'
         stop 'rinputd'
      endif

      go to 5

c     label: config

 80   continue
      
c     read total charge for the system

      call inpi(itotq)
      totq=z1+z2-dble(itotq)

c     to allow for noninteger nucleus charge needed to force convergence
c     in some open-shell cases read in the nuclei charges again

      z1t=0.d0
      z2t=0.d0
      call inpf(z1t)
      call inpf(z2t)
      if (z1t.gt.precis) z1=z1t
      if (z2t.gt.precis) z2=z2t

c     read in symmetry and occupation and test if
c     the orbital is locked = 1 (open shell case)
c                           = 0 (closed shell case)
      isum=0
 82   call card
      call inpi(nbsym)
      call inpa(char8)

      if(char8.ne.sigma.and.char8.ne.pi.and.
     &     char8.ne.delta.and.char8.ne.phi) then
	 write(iout6,*) 'Error: wrong symmetry or symmetry higher',
     &        ' than phi'
         stop 'rinputd'
      endif

      isum0=isum
      do i=1,nbsym
         isum=isum+1
         orbsym(isum)=char8
	 bond(isum)=orbsym(isum)
         gut(isum)=space
      enddo
      
      if (isum.gt.maxorb) then
         write(iout6,*) 'too many orbitals (max. see maxorb)'
         stop 'rinputd'
      endif
      
c     in homonuclear case read in the parity of each orbital 
c     if BREAK is on (ibreak=1) u/g labels are ignored
c     if |z1-z2|<homolevl then the case is treated as a homonuclear one
c     homolevl is set to 1.0e-6 (see blk_data.inc) 

      if (abs(z1-z2).lt.homolevl.and.ibreak.eq.0) then
         call inpa(gut(isum))
         isum1=isum0
         do i=1,nbsym
            isum1=isum1+1
            gut(isum1)=gut(isum)
         enddo
      endif
      
      cloe=4.d0
      if (orbsym(isum).eq.sigma) cloe=2.d0
      iput=4*(isum-1)
      call inpa(char8)
      if (char8.ne.endl) then
         spin(1+iput)=char8
         call inpa(char8)
         if (char8.ne.endl) then 
            spin(2+iput)=char8
            call inpa(char8)
            if (char8.ne.endl) then
               spin(3+iput)=char8
               call inpa(char8)
               if (char8.ne.endl) then
                  spin(4+iput)=char8
               endif
            endif
         endif
      endif

      if (spin(1+iput).ne.'+'.and.spin(1+iput).ne.'-') then
         do i=1,nbsym
            occ(isum0+i)= cloe
         enddo
      else
         do i=1,4
            if (spin(i+iput).eq.'+'.or.spin(i+iput).eq.'-') then
               do j=1,nbsym 
                  occ(isum0+j)=occ(isum0+j)+1.d0
               enddo
            endif
         enddo
      endif
      norb=isum
      no=norb
      
c     initialize scaling factors for off-diagonal Lagrange multipliers

      do i=1,norb
         do j=i+1,norb
            sflagrat(i,j)=sflagra
         enddo
      enddo

c     initialize damping factors and SOR parameters

      do i=1,norb
         dampf(i)=0.d0
         maxsor(i)=maxsor2
      enddo

      if (char8.eq.endl.or.char8.eq.endc) then
         totchar=0.d0
         do iorb=1,norb
            totchar=totchar+occ(iorb)
         enddo

c        Warning!
c        on Cray Y-MP (F90) the calculated charge was 29. and the one 
c        read in -- 29.00000000000011. 
c        The statement if (totchar.eq.totq) had to be replaced by the following
c        one

         if (abs(totchar-totq).gt.precis*1000.0) then
            write(iout6,*) 'Error: mismatch in given and',
     &           ' calculated total charge'
            stop 'rinputd'
         endif
         goto 5
      endif
      
      goto 82

c     label: initial

 100  continue

c     iform =0  read multiple exchange files &  write multiple exchange files
c     iform =1  read a single exchange file  &  write multiple exchange files
c     iform =2  read multiple exchange files &  write a single exchange file 
c     iform =3  read a single file           &  write a single exchange file 

      call inpi(ini)
      call inpi(iform)
      call inpi(itmp)
      if (itmp.ne.inpiexit) then
         inclorb=itmp
      else
         inclorb=0
      endif
      
      if (imethod.eq.2) then
         if (ini.ne.1) ini=4 
         facmul=10.0**(15)
      endif

      if (imethod.eq.3.and.iform.ne.3) then
         write(iout6,*) 'Error: if method is chosen as HFS ',
     &        'the second parameter has to be set to 3'
         stop 'rinputd'
      endif
      
      if (inclorb.eq.2) then

c        'hydrogen' initialization with screening

         do iorb=1,norb
            call card
            call inpf(co1(iorb))
            call inpi(mgx(1,iorb))
            call inpi(mgx(2,iorb))
            call inpf(eza1(iorb))
            eza1(iorb)=z1-eza1(iorb)
            call inpf(co2(iorb))
            call inpi(mgx(4,iorb))
            call inpi(mgx(5,iorb))
            call inpf(eza2(iorb))
            eza2(iorb)=z2-eza2(iorb)
            call inpi(ifix(iorb))
            call inpi(maxsor(iorb))
            if (maxsor(iorb).eq.0) maxsor(iorb)=maxsor2
            i1ng(iorb)=0
         enddo
      elseif (inclorb.eq.1) then

c        'hydrogen' initialization without screening

         do iorb=1,norb
            call card
            call inpf(co1(iorb))
            call inpi(mgx(1,iorb))
            call inpi(mgx(2,iorb))
            call inpf(eza1(iorb))
            eza1(iorb)=eza1(iorb)
            call inpf(co2(iorb))
            call inpi(mgx(4,iorb))
            call inpi(mgx(5,iorb))
            call inpf(eza2(iorb))
            eza2(iorb)=eza2(iorb)
            call inpi(ifix(iorb))
            call inpi(maxsor(iorb))
            if (maxsor(iorb).eq.0) maxsor(iorb)=maxsor2
            i1ng(iorb)=0
         enddo
      elseif (inclorb.eq.0) then
         
c        orbitals already initialized or taken from Gaussian output

         do iorb=1,norb
            maxsor(iorb)=maxsor2
            co1(iorb)=1.d0
            co2(iorb)=1.d0
         enddo
      else
         write(*,*) 'rinputd: Wrong value of the third parameter'
         stop 'rinputd'
      endif

      do iorb=1,norb
	if (i1ng(iorb).eq.0) then
	  i1ng(iorb)=ngrids
	  i2ng(iorb)=ngrids
	endif
      enddo

      do iorb=1,norb
         if (orbsym(iorb).eq.sigma) mt=0
         if (orbsym(iorb).eq.pi)    mt=1
         if (orbsym(iorb).eq.delta) mt=2
         if (orbsym(iorb).eq.phi)   mt=3

c        asymuthal quantum number of a hydrogen orbital is taken 
c        in accordance with the symmetry of a molecular orbital

         mgx(3,iorb)=mt
         mgx(6,iorb)=mt
         lock(iorb)=0

         clo=4.d0
         if (orbsym(iorb).eq.sigma) clo=2.d0
         clo=occ(iorb)/clo
	 if (abs(clo-1.d0).gt.1.d-06) then
	    lock(iorb)=1
	 endif

         mgi=(-1)**mgx(3,iorb)
      enddo

c     the program works on a set of standard input and output dump files
c     containing orbitals, Coulomb potentials (with extentions .orb and .pot, 
c     respectively), exchange potentials (fort.201, fort.202, ...) and
c     some additional data ( *.ext)

      psiinp ='2dhf_input.orb'
      coulinp='2dhf_input.coul'
      exchinp='2dhf_input.exch'
     
      psiout ='2dhf_output.orb'
      coulout='2dhf_output.coul'
      exchout='2dhf_output.exch'

c     open separate files for reading and writing orbitals and potentials 

      if (iform.eq.0) then
         open(11,file=psiinp, status='unknown',form='unformatted')
         open(12,file=coulinp,status='unknown',form='unformatted')
         iinp13=iout23

         open(21,file=psiout, status='unknown',form='unformatted')
         open(22,file=coulout,status='unknown',form='unformatted')
         
      elseif (iform.eq.1) then
         open(11,file=psiinp, status='unknown',form='unformatted')
         open(12,file=coulinp,status='unknown',form='unformatted')
         open(13,file=exchinp,status='unknown',form='unformatted')
         rewind(13)        

         open(21,file=psiout, status='unknown',form='unformatted')
         open(22,file=coulout,status='unknown',form='unformatted')

      elseif (iform.eq.2) then
         open(11,file=psiinp, status='unknown',form='unformatted')
         open(12,file=coulinp,status='unknown',form='unformatted')

         open(21,file=psiout, status='unknown',form='unformatted')
         open(22,file=coulout,status='unknown',form='unformatted')
         open(23,file=exchout,status='unknown',form='unformatted')
         rewind(23)

      elseif (iform.eq.3) then
         open(11,file=psiinp, status='unknown',form='unformatted')
         open(12,file=coulinp,status='unknown',form='unformatted')
         open(13,file=exchinp,status='unknown',form='unformatted')
         rewind(13)

         open(21,file=psiout, status='unknown',form='unformatted')
         open(22,file=coulout,status='unknown',form='unformatted')
         open(23,file=exchout,status='unknown',form='unformatted')
         rewind(23)
      endif

      rewind(11)
      rewind(12)
      rewind(21)
      rewind(22)

c     ini=5 enables to retrieve from complete/uncomplete dump file
c           if only one (the highest= the first) orbital is missing
c           the rest is retreived and the first is initialized as 
c           a hydrogenic one (ini is changed from 5 to 4)

      go to 5

c     label: break put switch on for breaking the symmetry

 120  continue
      ibreak=1
      go to 5

c     label: multipol -- define number of multipole moments to be calculated

 140  continue

c     facmul  -  if facmul>0  multipole moment expansion
c		 coefficients are recalculated	every time demax(1),
c		 i.e. maximum error in orbital energy, is reduced by
c		 facmul
c		 if facmul<0 these coefficients are not calculated
c
      call inpf(facmul)
      if (facmul.lt.0.d0) facmul=10.0**(15)
      go to 5

c     label: grid

 160  continue

c     call inpi(nni)
c     call inpi(nmu(ngrids))
c     call inpf(rgrid(ngrids))
      
      tmp1=0.d0
      tmp2=0.d0
      call inpi(nni)
      call inpf(tmp1)
      call inpf(tmp2)
      if (tmp2.ne.0.d0) then
c        nnu nmu R_infty
         nmu(ngrids)=tmp1
         rgrid(ngrids)=tmp2
      else
c        nni R_infty
         if (ngrids.ne.1) then
            write(*,*)'rinputd: missing item in GRID card'
            stop 'rinputd'
         endif
         rgrid(ngrids)=tmp1
         nmu(1)=nmucalc()
      endif
     
      rinf=rgrid(ngrids)
      nmutot=nmu(ngrids)

      if (nni.gt.maxni) then
         write(*,*) 'Error: too many grid points in nu variable'
         stop 'rinputd'
      endif
      if (nmutot.gt.maxmu) then
         write(*,*) 'Error: too many grid points in mu variable'
         stop 'rinputd'
      endif

      goto 5

c     label: subgrid

c     up to 10 grids can be formally defined (rgrid(10) though
c     rgrid(9) would sufice nmu(10))
c     the last value of rgrid must be equal to rinf, i.e. the practical infinity

 162  continue

      call inpi(ngrids)
      call card
      call inpi(nni)
      nmutot=0
      do ig=1,ngrids
        call inpi(nmu(ig))
        nmutot=nmutot+nmu(ig)
      enddo
 
      if (nni.gt.maxni) then
         write(*,*) 'Error: too many grid points in nu variable'
         stop 'rinputd'
      endif
      if (nmutot.gt.maxmu) then
         write(*,*) 'Error: too many grid points in mu variable'
         stop 'rinputd'
      endif

      call card
      do ig=1,ngrids
        call inpf(rgrid(ig))
      enddo
      rinf=rgrid(ngrids)

      if (ngrids.gt.maxgrids) then
         write(*,*) 'Error: too many subgrids'
         stop 'rinputd'
      endif

c     set default ordering for subgrids

      do ig=1,ngrids
         iorder(ig)=2
      enddo      
      goto 5

c     label: cutoff -- cutoff values for orbitals and potentials

 180  continue
      call inpi(itmp)
      if (itmp.ge.0) itail=itmp
      goto 5

c     label: homo --  when this label is encountered the g/u symmetry of
c                     orbitals is strictly imposed; in all other respects
c                     a homonuclear case is treated as a heteronuclear one
c     ihomon=0 --  heteronuclear case
c     ihomon=1 --  homonuclear case |z1-z2|<homolevl (set in common block data)
c     ihomon=2 --  homonuclear case with forced symmetry

 200  continue
      ihomon=2
      goto 5

c     label: scf
c
c     maxscf  - max. no. of scf iterations
c     nobckup - no. of iterations between backups 
c     ienterm - if max. error in orbital energy is less than 
c               1/10**ienterm than scf iterations are terminated
c     inoterm - if max. error in orbital norm is less than 
c               1/10**inoterm than scf iterations are terminated
c     iprtlev - level of output during scf process
c               conv. rate, orbital energy, normalization of the worst
c               converged orbital is printed every nobckup iterations
c               if iprtlev=3 or every scf iteration if iprtlev=2 (default).
c               if iprtlev=1 then conv. rate, orbital energy and normalization 
c               of every orbital is printed in every scf iteration;
c               in any case total energy is printed every nobckup iterations;

 220  continue
      call inpi(itmp)
      if (itmp.ne.inpiexit) then
         maxscf = itmp
         call inpi(itmp)
         if (itmp.ne.inpiexit) then
            nobckup = itmp
            call inpi(itmp)
            if (itmp.ne.inpiexit) then
               ienterm = itmp
               call inpi(itmp)
               if (itmp.ne.inpiexit) then
                  inoterm = itmp
                  call inpi(itmp)
                  if (itmp.ne.inpiexit) then
                     iprtlev=itmp
                  endif
               endif
            endif
         endif
      endif
      goto 5

c     label: fix

c     exlwf exlpot- if nonzero wave functions, coulomb potentials or
c     exlexp        exchange potentials are not relaxed; exlpot and 
c                   exlexp cannot be non zero if hydrogen orbitals 
c                   are used to start the scf process
 222  continue

      call inpi(itmp)
      if (itmp.ne.inpiexit) then
         iexlorb=itmp
         call inpi(itmp)
         if (itmp.ne.inpiexit) then
            iexlcoul=itmp
            call inpi(itmp)
            if (itmp.ne.inpiexit) then
               iexlexp=itmp
            endif
         endif
      endif

      exlorb=dble(iexlorb)
      exlcoul=dble(iexlcoul)
      exlexp=dble(iexlexp)
      goto 5

c     label: xalpha
c
c      alphaf  -  if alphaf is not zero the slater exchange approximation
c		 is used as long as demax(1) exceeds the value of sltthr
c      sltthr -  if demax(1) is less then sltthr the local slater
c		 approximation is replaced by the full exchange

 224  continue
      call inpf(ftmp)
      if (ftmp.ne.0.d0) then
         alphaf = ftmp
         call inpf(ftmp)
         if (ftmp.ne.0.d0) then
            sltthr = ftmp
         endif
      endif
      if (alphaf.ne.0.d0) then
         exlexp=1.d0
         islat =1
      endif
      goto 5

c     label: sor

c     maxsor2 - maximal number of (mc)sor iterations (over each subgrid)
c               during relaxation of every orbital and potential in an SCF cycle

c     lsor    - a scaling factor used to change maxsor2 value when relaxing
c               potentials (see relcoul1/2 and relexch1/2); the number of (mc)sor
c               iterations for the potentials is calculalated as maxsor2/lsor

c     ipoiss  -  choice of poisson's equation solving method
c        = 1     sor (modified LPS's code)
c        = 2     mcsor - variables are relaxed in 5 sweeps each relaxing
c                every 5th variable (145 mflops on Cray Y-MP)

 240  continue
      call inpi(itmp)
      if (itmp.ne.inpiexit) then
         maxsor2 = itmp
         call inpi(itmp)
         if (itmp.ne.inpiexit) then
            lsor = itmp
            call inpi(itmp)
            if (itmp.ne.inpiexit) then
               ipoiss = itmp                       
c     maxsor1  -  maximal number of (mc)sor iterations over the whole grid
c                 between which boundary value at infinity are recalculated

c               call inpi(itmp)
c               if (itmp.ne.inpiexit) then
c                  maxsor1 = itmp
c               endif
            endif
         endif
      endif

      goto 5
      
c     label: order

c     iorder  -  type of (mc)sor sweeps (ordering) for each grid
c	= 1   -  natural column-wise ordering for the mesh points
c	= 2   -  'middle' type of sweep - default of 2d-hf program
c	= 3   -  natural row-wise
c	= 4   -  reversed natural (column-wise) ordering (see meshwk for details)
c

 242  continue
      do ig=1,ngrids
         call inpi(iorder(ig))
      enddo
      goto 5

c     label: omega

c     ovforb   -  overelaxation parameter for orbitals and grids
c     ovfcoul  -  overelaxation parameter for coulomb pot. and grids
c     ovfexch  -  overelaxation parameter for exchange pot. and grids

 244  continue
      call card
      call inpf(tmp)
      if (tmp.gt.0.d0) then
         ovforb(1)=tmp
      else
         write(*,*) 'Error: missing or incorrect entry'
         stop 'rinputd'
      endif
      do ig=2,ngrids
        ovforb(ig) = ovforb(1)
      enddo      
      call inpf(tmp)
      if (tmp.gt.0.d0) then
         ovforb(2)=tmp
         call inpf(tmp)
         if (tmp.gt.0.d0) then
            ovforb(3)=tmp
         endif
      endif
      
      call card

      call inpf(tmp)
      if (tmp.ne.0.d0) then
         ovfcoul(1)=tmp
      else
         write(*,*) 'Error: missing or incorrect entry'
         stop 'rinputd'
      endif
      do ig=2,ngrids
        ovfcoul(ig) = ovfcoul(1)
      enddo
      call inpf(tmp)
      if (tmp.gt.0.d0) then
         ovfcoul(2)=tmp
         call inpf(tmp)
         if (tmp.gt.0.d0) then
            ovfcoul(3)=tmp
         endif
      endif
      
      do ig=1,ngrids
        ovfexch(ig) = ovfcoul(ig)
      enddo
      goto 5

c     label: plot
      
 340  continue
      stop 'rinputd: plot disabled'
c      call inpi(itmp)
c      if (itmp.ne.inpiexit) iplot=itmp
c      goto 5
      
c     label: fermi

 360  continue
      call inpf(z1atmass)
      call inpf(z2atmass)
      if (z1atmass.le.0.d0.and.z2atmass.le.0.d0) then
         write(*,*) 'Error: missing or incorrect entry'
         stop 'rinputd'
      endif
      ifermi=1
      goto 5

c     label: gauss

 362  continue
      call inpf(z1atmass)
      call inpf(z2atmass)
      if (z1atmass.le.0.d0.and.z2atmass.le.0.d0) then
         write(*,*) 'Error: missing or incorrect entry'
         stop 'rinputd'
      endif
      ifermi=2
      goto 5

c     label: lagra (not used)
c     scaling factor for off-diagonal Lagrange multipliers

 370  continue
c      call inpf(sflagra)
      stop 'rinputd: lagra disabled'
      goto 5

c     label: damp (not used)
c     damping factor for smoothing out oscillations in orbital values
c     caused by orthogonalization

 380  continue

      stop 'rinputd: damp disabled'

c      do i=1,norb
c         call card
c         call inpf(damp)
c         dampf(i)=damp
c      enddo

      goto 5

c     label: canon
c     choice of canonication coefficients due to
c     Davidson (default)       - 1
c     Guest and Saunders       - 2
c     Roothaan (single matrix) - 3
c     McWeeny and Diercksen    - 4

 390  continue

      call inpi(icanon)
      goto 5
      
c     label: fefield
c     finite field calculations
c     1 field strength
c     2 gradient strength
      
 392  continue
      ifefield=1
c      call inpi(ifefield)
c      if     (ifefield.eq.1) then 
      call inpf(ffield)
      call inpf(zcutoff)
c         write(*,*) 'ifefield,ffield ',ifefield,ffield
c      elseif (ifefield.eq.2) then 
c         call inpf(fgrad)
c         call inpf(zcutoff)
c      else
c         write(*,*) 'Error: incorrect type of finite ',
c     &        'field calculations'
c         stop 'rinputd'
c      endif
      goto 5

c     label: debug

c     idbg    -  additional printout
c     incrni  -  parameters of pmtx routine used to print matrice
c     incrmi	 row-wise. every incrni row and every incrmi column
c		 are printed

c     idbg  -  additional printout for testing purposes

c     set debug flags. if an integer i is encountered debuf flag i
c     is set, i.e. idbg(i)=1. at most 15 integers (1<i<100) can
c     be specified. if fewer integers are in the record terminate
c     it by a slash.

 260  continue

      do i=1,600
         id(i)=0
      enddo

      inzero=0
      do i=1,20
         call inpi(id(i))
         if (id(i).gt.0) then
           inzero=inzero+1
           idbg(id(i))=1
        endif
      enddo
      write(iout6,1110) (id(i),i=1,inzero)
      if (idbg(300).ne.0.or.idbg(301).ne.0) then
         write(*,*) "rinputd: Warning! ",
     &        "Do you realy want to use this flag?!"
      endif
 1110 format(10i7)
      
      if (idbg(600).eq.1) then
         call card
         call inpi(incrni)
         call inpi(incrmu)
      endif
      goto 5
      
      if (idbg(599).eq.1) then
         call card
         call inpi(iorb)
         read(*,*) (sflagrat(iorb,i),i=iorb+1,norb)
         
         do i=iorb+1,norb
            sflagrat(i,iorb)=sflagrat(iorb,i)/2.d0
         enddo
         write(*,'(10d12.2)') (sflagrat(iorb,i),i=iorb+1,norb)
      endif

c  label: help 

 280  continue
      write(*,*) '... no online help is available ...'
      write(*,*) '... see ./doc/users_guide.html or ',
     &     './doc/users_guide.pdf ...'
      stop 'rinputd'
      goto 5

c     label: stop

 300  continue

      write(*,*) '  ... end of input data  ...'
      write(*,*)'        '

      iform_t = iform
      ni_t    = nni
      mu_t    = nmutot
      no_t=norb
      nons_t=0
      do iorb=1,norb
         if (mgx(3,iorb).ne.0) nons_t=nons_t+1
      enddo
      return

c     label: interp

 320  continue

      iflag=0
      call inpi(itmp)
      if (itmp.ne.inpiexit) then
         iord_nu_orb=itmp+1
         call inpi(iord_mu_orb)
         iord_mu_orb=iord_mu_orb+1
         
         call inpi(iord_nu_coul)
         iord_nu_coul=iord_nu_coul+1
         call inpi(iord_mu_coul)
         iord_mu_coul=iord_mu_coul+1

         call inpi(iord_nu_exch)
         iord_nu_exch=iord_nu_exch+1
         call inpi(iord_mu_exch)
         iord_mu_exch=iord_mu_exch+1
      endif
      if ((iord_nu_orb.ne.3).and.(iord_nu_orb.ne.5).and.
     &     (iord_nu_orb.ne.7).and.(iord_nu_orb.ne.9)) goto 1600
      
      if ((iord_mu_orb.ne.3).and.(iord_mu_orb.ne.5).and.
     &     (iord_mu_orb.ne.7).and.(iord_mu_orb.ne.9)) goto 1600
      
      if ((iord_nu_coul.ne.3).and.(iord_nu_coul.ne.5).and.
     &     (iord_nu_coul.ne.7).and.(iord_nu_coul.ne.9)) goto 1600
      
      if ((iord_mu_coul.ne.3).and.(iord_mu_coul.ne.5).and.
     &     (iord_mu_coul.ne.7).and.(iord_mu_coul.ne.9)) goto 1600
      
      if ((iord_nu_exch.ne.3).and.(iord_nu_exch.ne.5).and.
     &     (iord_nu_exch.ne.7).and.(iord_nu_exch.ne.9)) goto 1600
      
      if ((iord_mu_exch.ne.3).and.(iord_mu_exch.ne.5).and.
     &     (iord_mu_exch.ne.7).and.(iord_mu_exch.ne.9)) goto 1600
      
      iinterp=1
      goto 5

01500 write(iout6,*) 'I/O error when reading header'
      stop 'rinputd'
01600 continue
      write(*,1120)
01120 format(1x,'Error: order of interpolation polynomial must be',
     &     ' 2, 4, 6 or 8')
      stop 'rinputd'
      end