      subroutine input(natoms,ncenters,nptcharges,ilocalize,ibfrozen,
     &                 iafrozen,idfrozen,iifrozen,nfunctional,ngridtype,
     &                 maxiterations,maxgeometries,iptsflag,icharge,
     &                 imultiplicity,idiis,nfdpbt,irestart,imesp,
     &                 ioptimization,ivibrations,intcoordflag,imode,
     &                 nequivalences,nbneutrals,nbsystems,idiffuse,
     &                 ixcfit,iptsweight,izmat1,izmat2,izmat3,natomtype,
     &                 nequivvec,neutral,isys,ibfreeze,iafreeze,
     &                 idfreeze,iifreeze,ibmode,iamode,idmode,iimode,
     &                 econvergence,gconvergence,dmixing,dmaxstep,
     &                 elevelshift,epsilon,cutoff,coord,charge,
     &                 coordptcharge,ptcharge,atomsymbol)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.
c
c  based on steve chin and fiona sim's decode subroutine

      implicit real*8(a-h,o-z)

      parameter(noptions=31,maxitems=40)

      character char_string*79,parsed(maxitems)*32,options(noptions)*4

      dimension nchars(maxitems)

      dimension value(3)

      dimension izmat1(*),izmat2(*),izmat3(*),natomtype(*),
     &          nequivvec(natoms,*),neutral(natoms,*),isys(natoms,*),
     &          ibfreeze(2,*),iafreeze(3,*),idfreeze(4,*),iifreeze(4,*),
     &          ibmode(*),iamode(*),idmode(*),iimode(*)

      dimension cutoff(*),coord(3,*),charge(*),coordptcharge(3,*),
     &          ptcharge(*)

      character atomsymbol(*)*4

      data zero,bohr,one,evolt/0.0d0,0.529177d0,1.0d0,27.2116d0/





      data options/ 'CART', 'ZMAT', 'RUNT', 'POTE', 'MULT',
     &              'CHAR', 'MAXI', 'DCON', 'GCON', 'MAXS',
     &              'MIXI', 'DIIS', 'LEVE', 'MAXG', 'EQUI',
     &              'POIN', 'FDPB', 'MODE', 'ORBI', 'NEUT',
     &              'CUTO', 'SUBS', 'MESP', 'QUAL', 'BFRE',
     &              'AFRE', 'DFRE', 'IFRE', 'DIFF', 'XCFI',
     &              'WEIG'/





 1001 call parser(numwords,nchars,parsed,char_string)





      do 1002 i=1,noptions

      if(parsed(1)(1:1).eq.'#') then
                                      write(6,1003) char_string
 1003                                 format(a79)
                                      goto 1001
                                endif

      if(parsed(1)(1:4).eq.options(i)) goto 1004
      if(parsed(1)(1:3).eq.'END')  goto 1005

 1002 continue

      stop 'keyword error:  please check your input file'





 1004 continue

      goto( 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100,
     &      2110, 2120, 2130, 2140, 2150, 2160, 2170, 2180, 2190, 2200,
     &      2210, 2220, 2230, 2240, 2250, 2260, 2270, 2280, 2290, 2300,
     &      2310) i





 2010 if(parsed(2)(1:3).eq.'ANG') then

           iangstroms=1

           write(6,2011)
 2011      format(/,' input in cartesian coordinates and angstroms')

                                  else

           iangstroms=0

           write(6,2012)
 2012      format(/,' input in cartesian coordinates and bohr radii')

                                  endif

      intcoordflag=0

 2013 call parser(numwords,nchars,parsed,char_string)

      if(parsed(1)(1:3).eq.'END') goto 1001

      ncenters=ncenters+1

      atomsymbol(ncenters)=parsed(1)(1:4)

      call element(ncenters,atomsymbol(ncenters),natomtype)

      if(iangstroms.eq.1) then

           coord(1,ncenters)=getvalue(parsed(2),nchars(2))/bohr
           coord(2,ncenters)=getvalue(parsed(3),nchars(3))/bohr
           coord(3,ncenters)=getvalue(parsed(4),nchars(4))/bohr

                          else

           coord(1,ncenters)=getvalue(parsed(2),nchars(2))
           coord(2,ncenters)=getvalue(parsed(3),nchars(3))
           coord(3,ncenters)=getvalue(parsed(4),nchars(4))

                          endif

      goto 2013





 2020 if(parsed(2)(1:3).eq.'ANG') then

           iangstroms=1

           write(6,2021)
 2021      format(/,' input in zmatrix form and angstroms')

                                  else

           iangstroms=0

           write(6,2022)
 2022      format(/,' input in zmatrix form and bohr radii')

                                  endif

      intcoordflag=1

 2023 call parser(numwords,nchars,parsed,char_string)

      if(parsed(1)(1:3).eq.'END') goto 1001

      ncenters=ncenters+1

      atomsymbol(ncenters)=parsed(1)(1:4)

      call element(ncenters,atomsymbol(ncenters),natomtype)

      if(ncenters.eq.1) call genxyz(1,0,0,0,zero,zero,zero,coord)

      if(ncenters.eq.2) then

           izmat1(ncenters)=int(getvalue(parsed(2),nchars(2)))

           if(iangstroms.eq.1) then

                           value(1)=getvalue(parsed(3),nchars(3))/bohr

                               else

                           value(1)=getvalue(parsed(3),nchars(3))

                               endif

           call genxyz(2,izmat1(2),0,0,value(1),zero,zero,coord)

                        endif

      if(ncenters.eq.3) then

           izmat1(ncenters)=int(getvalue(parsed(2),nchars(2)))
           izmat2(ncenters)=int(getvalue(parsed(4),nchars(4)))

           if(iangstroms.eq.1) then

                           value(1)=getvalue(parsed(3),nchars(3))/bohr

                               else

                           value(1)=getvalue(parsed(3),nchars(3))

                               endif

           value(2)=getvalue(parsed(5),nchars(5))

           call genxyz(3,izmat1(3),izmat2(3),0,value(1),value(2),
     &                 zero,coord)

                        endif

      if(ncenters.gt.3) then

           izmat1(ncenters)=int(getvalue(parsed(2),nchars(2)))
           izmat2(ncenters)=int(getvalue(parsed(4),nchars(4)))
           izmat3(ncenters)=int(getvalue(parsed(6),nchars(6)))

           if(iangstroms.eq.1) then

                           value(1)=getvalue(parsed(3),nchars(3))/bohr

                               else

                           value(1)=getvalue(parsed(3),nchars(3))

                               endif

           value(2)=getvalue(parsed(5),nchars(5))
           value(3)=getvalue(parsed(7),nchars(7))

           call genxyz(ncenters,izmat1(ncenters),izmat2(ncenters),
     &                 izmat3(ncenters),value(1),value(2),value(3),
     &                 coord)

                      endif

      goto 2023




 2030 continue

      do 2031 iword=2,3

      if(parsed(iword)(1:4).eq.'STAR') then

                    irestart=0

                    write(6,2032)
 2032               format(/,' starting up an scf calculation ')

                                       endif

      if(parsed(iword)(1:4).eq.'REST') then

                    irestart=1

                    write(6,2033)
 2033               format(/,' restarting an scf calculation ')

                                       endif

      if(parsed(iword)(1:3).eq.'OPT')  then

                    ioptimization=1

                    write(6,2034)
 2034               format(/,' starting a geometry optimization ')

                                       endif

      if(parsed(iword)(1:4).eq.'CONT') then

                    ioptimization=1

                    irestart=2

                    write(6,2035)
 2035               format(/,' restarting a geometry optimization ')

                                       endif

      if(parsed(iword)(1:4).eq.'VIBR') then

                    ivibrations=1

                    irestart=2

                    write(6,2036)
 2036               format(/,' performing a vibrational analysis ')

                                       endif

      if(parsed(iword)(1:4).eq.'MESP') then

                    imesp=1

                    irestart=3

                    write(6,2037)
 2037               format(/,' going straight to an',
     &                       ' evaluation of the MESP ')

                                       endif

 2031 continue

      goto 1001





 2040 if(parsed(2)(1:4).eq.'BECK') then

        nfunctional=1

        write(6,2041)
 2041   format(/,' using becke exchange/perdew correlation functional ')

        goto 1001

                                   endif

      if(parsed(2)(1:4).eq.'MIXB') then

        nfunctional=2

        write(6,2042)
 2042   format(/,' using becke exchange/perdew correlation energy',
     &           ' with the vwn potential ')

        goto 1001

                                   endif

      nfunctional=0

      write(6,2043)
 2043 format(/,' using the vwn (lsda) functional ')

      goto 1001





 2050 imultiplicity=getvalue(parsed(2),nchars(2))

      write(6,2051) imultiplicity
 2051 format(/,' the multiplicity of the system is ',i3)

      goto 1001





 2060 icharge=int(getvalue(parsed(2),nchars(2)))

      write(6,2061) icharge
 2061 format(/,' the charge on the system is ',i4)

      goto 1001





 2070 maxiterations=int(getvalue(parsed(2),nchars(2)))

      write(6,2071) maxiterations
 2071 format(/,' the maximum number of iterations per geometry is ',i4)

      goto 1001





 2080 econvergence=getvalue(parsed(2),nchars(2))

      write(6,2081) econvergence
 2081 format(/,' the scf convergence criterion is set at ',e14.7)

      goto 1001





 2090 gconvergence=getvalue(parsed(2),nchars(2))

      write(6,2091) gconvergence
 2091 format(/,' the geometry optimization criterion is set at ',e14.7)

      goto 1001





 2100 dmaxstep=getvalue(parsed(2),nchars(2))

      write(6,2101) dmaxstep
 2101 format(/,' maximum stepsize in geometry optimization set at ',
     &       f10.6,' bohr ')

      goto 1001





 2110 dmixing=getvalue(parsed(2),nchars(2))

      write(6,2111) dmixing
 2111 format(/,' using an scf mixing parameter of ',f6.3)

      goto 1001





 2120 if(parsed(2)(1:3).eq.'OFF') then

               idiis=0

               write(6,2121)
 2121          format(/,' the diis scheme has been turned off ')

                                  else

               idiis=1

               write(6,2122)
 2122          format(/,' the diis scheme has been turned on ')

                                  endif

      goto 1001





 2130 elevelshift=getvalue(parsed(2),nchars(2))

      write(6,2131) elevelshift
 2131 format(/,' the levelshift parameter is set to ',f10.3,
     &         ' electron volts ')

      elevelshift=elevelshift/evolt

      goto 1001





 2140 maxgeometries=int(getvalue(parsed(2),nchars(2)))

      write(6,2141) maxgeometries
 2141 format(/,' the maximum number of geometries is ',i4)

      goto 1001





 2150 continue

      do 2151 i=3,numwords

      nequivalences=nequivalences+1

      nequivvec(nequivalences,1)=int(getvalue(parsed(2),nchars(2)))
      nequivvec(nequivalences,2)=int(getvalue(parsed(i),nchars(i)))

      write(6,2152) nequivvec(nequivalences,1),
     &              nequivvec(nequivalences,2)

 2152 format(/,' atoms # ',i4,' and ',i4,' are constrained to have',
     &         ' the same esp fitted charge ')

 2151 continue

      goto 1001





 2160 if(parsed(2)(1:3).eq.'ANG') then

            iangstroms=1

            write(6,2161)
 2161       format(/,' the positions of the point charges are given',
     &               ' in angstroms ')

                                  else

            iangstroms=0

            write(6,2162)
 2162       format(/,' the positions of the point charges are given',
     &               ' in bohr ')

                                  endif

 2163 call parser(numwords,nchars,parsed,char_string)

      if(parsed(1)(1:3).eq.'END') goto 1001

      nptcharges=nptcharges+1

      if(iangstroms.eq.1) then

          coordptcharge(1,nptcharges)=getvalue(parsed(1),nchars(1))/bohr
          coordptcharge(2,nptcharges)=getvalue(parsed(2),nchars(2))/bohr
          coordptcharge(3,nptcharges)=getvalue(parsed(3),nchars(3))/bohr

                          else

          coordptcharge(1,nptcharges)=getvalue(parsed(1),nchars(1))
          coordptcharge(2,nptcharges)=getvalue(parsed(2),nchars(2))
          coordptcharge(3,nptcharges)=getvalue(parsed(3),nchars(3))

                          endif

      ptcharge(nptcharges)=getvalue(parsed(4),nchars(4))

      goto 2163





 2170 if(parsed(2)(1:3).eq.'ON') then

           nfdpbt=1

           epsilon=getvalue(parsed(3),nchars(3))

           write(6,2171) epsilon
 2171      format(/,' a dft-fdpb calculation will be performed,',
     &              ' with a dielectric of ',f8.3)

                                 else

           nfdpbt=0

           epsilon=one

           write(6,2172)
 2172      format(/,' a dft-fdpb calculation will not be performed ')

                                 endif

      goto 1001





 2180 continue

      if(parsed(2)(1:1).eq.'B') then

                imode=1

                ibmode(1)=getvalue(parsed(3),nchars(3))
                ibmode(2)=getvalue(parsed(4),nchars(4))

                write(6,2181) ibmode(1),ibmode(2)
 2181           format(/,' will try to find transition state',
     &                   ' along bond ',2i5)

                                endif

      if(parsed(2)(1:1).eq.'A') then

                imode=2

                iamode(1)=getvalue(parsed(3),nchars(3))
                iamode(2)=getvalue(parsed(4),nchars(4))
                iamode(3)=getvalue(parsed(5),nchars(5))

                write(6,2182) iamode(1),iamode(2),iamode(3)
 2182           format(/,' will try to find transition state',
     &                   ' along bond angle ',3i5)

                                endif

      if(parsed(2)(1:1).eq.'D') then

                imode=3

                idmode(1)=getvalue(parsed(3),nchars(3))
                idmode(2)=getvalue(parsed(4),nchars(4))
                idmode(3)=getvalue(parsed(5),nchars(5))
                idmode(4)=getvalue(parsed(6),nchars(6))

                write(6,2183) idmode(1),idmode(2),idmode(3),idmode(4)
 2183           format(/,' will try to find transition state',
     &                   ' along dihedral angle ',4i5)

                                endif

      if(parsed(2)(1:1).eq.'I') then

                imode=4

                iimode(1)=getvalue(parsed(3),nchars(3))
                iimode(2)=getvalue(parsed(4),nchars(4))
                iimode(3)=getvalue(parsed(5),nchars(5))
                iimode(4)=getvalue(parsed(6),nchars(6))

                write(6,2184) iimode(1),iimode(2),iimode(3),iimode(4)
 2184           format(/,' will try to find transition state',
     &                   ' along improper torsion ',4i5)

                                endif
         
      goto 1001





 2190 continue

      if(parsed(2)(1:3).eq.'CAN') then

               ilocalize=1

               write(6,2191)
 2191          format(/,' canonical molecular orbitals will be',
     &                  ' generated')
 

                                  endif

      if(parsed(2)(1:3).eq.'BOY') then

               ilocalize=3

               write(6,2192)
 2192          format(/,' boys localized molecular orbitals will be',
     &                  ' generated')

                                  endif

      if(parsed(2)(1:3).eq.'BOT') then

               ilocalize=2

               write(6,2193)
 2193          format(/,' canonical molecular orbitals will be',
     &                  ' generated',/,
     &                  ' boys localized molecular orbitals will be',
     &                  ' generated')

                                  endif

      goto 1001





 2200 write (6,2201)
 2201 format(/,' constrained atomic charges (sum is zero)')
 
      nbneutrals=nbneutrals+1

      write(6,2202) nbneutrals
 2202 format(/,'       group # ',i3)

      neutral(nbneutrals,1)=numwords-1
  
      do 2204 i=2,numwords
      neutral(nbneutrals,i)=getvalue(parsed(i),nchars(i))
      write(6,2205) neutral(nbneutrals,i)
 2205 format('             atom #',i4)
 2204 continue

      write(6,2206)
 2206 format(//)

      goto 1001





 2210 continue

      cutoff(1)=getvalue(parsed(2),nchars(2))
      cutoff(2)=getvalue(parsed(3),nchars(3))
      cutoff(3)=getvalue(parsed(4),nchars(4))

      write(6,2211) cutoff(1),cutoff(2),cutoff(3)
 2211 format(/,' buffer atom space (CD fitting) cutoff set at ',
     &       f8.2,' Angstroms',/,
     &         ' buffer atom space (XC fitting) cutoff set at ',
     &       f8.2,' Angstroms',/,
     &         ' overlap/fock matrix element cutoff set at    ',
     &       f8.2,' Angstroms')

      cutoff(1)=cutoff(1)/bohr
      cutoff(2)=cutoff(2)/bohr
      cutoff(3)=cutoff(3)/bohr

      goto 1001





 2220 continue

      write(6,2222) nbsystems
 2222 format(/,' subsystem # ',i3)

      nbsystems=nbsystems+1

      isys(1,nbsystems-1)=numwords-1

      do 2223 i=2,numwords
      index=int(getvalue(parsed(i),nchars(i)))
      isys(i,nbsystems-1)=index
 2223 write(6,2224) index
 2224 format('           atom #',i4)

      write(6,2225)
 2225 format(//)

      goto 1001





 2230 if(parsed(2)(1:4).eq.'READ') then

               imesp=1
               iptsflag=0

               write(6,2231)
 2231          format(/,' the points used to fit the mesp will',
     &                  ' be read in from a user-supplied file ')

                                   else

               imesp=1
               iptsflag=1

               write(6,2232)
 2232          format(/,' the default set of points',
     &                  ' will be used to fit the mesp ')

                                   endif

      goto 1001





 2240 ngridtype=int(getvalue(parsed(2),nchars(2)))

      write(6,2241) ngridtype
 2241 format(/,' the following grid quality has been chosen: ',i4)

      goto 1001





 2250 ibfrozen=ibfrozen+1

      ibfreeze(1,ibfrozen)=int(getvalue(parsed(2),nchars(2)))
      ibfreeze(2,ibfrozen)=int(getvalue(parsed(3),nchars(3)))

      write(6,2251) ibfreeze(1,ibfrozen),ibfreeze(2,ibfrozen)
 2251 format(/,' the following bond will be frozen: ',2i5)

      goto 1001





 2260 iafrozen=iafrozen+1

      iafreeze(1,iafrozen)=int(getvalue(parsed(2),nchars(2)))
      iafreeze(2,iafrozen)=int(getvalue(parsed(3),nchars(3)))
      iafreeze(3,iafrozen)=int(getvalue(parsed(4),nchars(4)))

      write(6,2261) iafreeze(1,iafrozen),iafreeze(2,iafrozen),
     &              iafreeze(3,iafrozen)
 2261 format(/,' the following angle will be frozen: ',3i5)

      goto 1001





 2270 idfrozen=idfrozen+1

      idfreeze(1,idfrozen)=int(getvalue(parsed(2),nchars(2)))
      idfreeze(2,idfrozen)=int(getvalue(parsed(3),nchars(3)))
      idfreeze(3,idfrozen)=int(getvalue(parsed(4),nchars(4)))
      idfreeze(4,idfrozen)=int(getvalue(parsed(5),nchars(5)))

      write(6,2271) idfreeze(1,idfrozen),idfreeze(2,idfrozen),
     &              idfreeze(3,idfrozen),idfreeze(4,idfrozen)
 2271 format(/,' the following dihedral will be frozen: ',4i5)

      goto 1001





 2280 iifrozen=iifrozen+1

      iifreeze(1,iifrozen)=int(getvalue(parsed(2),nchars(2)))
      iifreeze(2,iifrozen)=int(getvalue(parsed(3),nchars(3)))
      iifreeze(3,iifrozen)=int(getvalue(parsed(4),nchars(4)))
      iifreeze(4,iifrozen)=int(getvalue(parsed(5),nchars(5)))

      write(6,2281) iifreeze(1,iifrozen),iifreeze(2,iifrozen),
     &              iifreeze(3,iifrozen),iifreeze(4,idfrozen)
 2281 format(/,' the following improper torsion will be frozen: ',4i5)

      goto 1001





 2290 idiffuse=int(getvalue(parsed(2),nchars(2)))

      if(idiffuse.eq.0) then

           write(6,2291)
 2291      format(/,' diffuse functions will not be added ')

                        endif

      if(idiffuse.eq.1) then

           write(6,2292)
 2292      format(/,' diffuse functions will be added to heavy atoms ')
 
                        endif

      if(idiffuse.eq.2) then

           write(6,2293)
 2293      format(/,' diffuse functions will be added to all atoms ')
 
                        endif

      goto 1001





 2300 if(parsed(2)(1:3).eq.'OFF') then

               ixcfit=0

               write(6,2301)
 2301          format(/,' no xc fitting functions will be used ')

                                  else

               ixcfit=1

               write(6,2302)
 2302          format(/,' xc fitting functions will be used ')

                                  endif

      goto 1001





 2310 if(parsed(2)(1:3).eq.'OFF') then

               iptsweight=0

               write(6,2311)
 2311          format(/,' grid point weight derivatives not',
     &                  ' included in gradient calculation ')

                                  else

               iptsweight=1

               write(6,2312)
 2312          format(/,' grid point weight derivatives',
     &                  ' included in gradient calculation ')

                                  endif

      goto 1001





 1005 continue

      do 1006 i=1,ncenters
 1006 call makelowercase(atomsymbol(i),4)

      do 1007 i=1,ncenters
 1007 charge(i)=real64(natomtype(i))

      return
      end
