      subroutine inputbases(natoms,ncenters,nconts,ncontp,ncontd,
     &                      ncontractions,ncds,ncdspd,ncdfuncs,nxcs,
     &                      nxcspd,nxcfuncs,idiffuse,ixcfit,iwkvec1,
     &                      iwkvec2,iwkvec3,iwkvec4,icfunc,ilfunc,
     &                      ngaussians,mtloca,nshels,nshelp,nsheld,
     &                      icdcfunc,icdlfunc,ixccfunc,ixclfunc,nt,
     &                      alpha,coeff,alphacd,coefscd,coefpcd,coefdcd,
     &                      alphaxc,coefsxc,coefpxc,coefdxc,bstring)

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.

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

      character input_string*30,test_string*30,bstring(natoms,2)*30

      dimension iwkvec1(*),iwkvec2(*),iwkvec3(*),iwkvec4(*),icfunc(*),
     &          ilfunc(*),ngaussians(*),mtloca(*),nshels(*),nshelp(*),
     &          nsheld(*),icdcfunc(*),icdlfunc(*),ixccfunc(*),
     &          ixclfunc(*),nt(*)

      dimension alpha(*),coeff(*),alphacd(*),coefscd(*),coefpcd(*),
     &          coefdcd(*),alphaxc(*),coefsxc(*),coefpxc(*),coefdxc(*)

      data pi/3.1415926535898d0/
      data zero,quarter,thrhlf,fivhlf/0.0d0,0.25d0,1.5d0,2.5d0/
      data one,two,three,eight,twelve,twenty/1.0d0,2.0d0,3.0d0,
     &                                       8.0d0,12.0d0,20.0d0/
      data onehund28,twothou48/128.0d0,2048.0d0/

      ncds=0
      nxcs=0

      ncdspd=0
      nxcspd=0

      n1=0
      n2=0

      read(5,2001) test_string

      do 1001 i=1,ncenters

      if(test_string(1:2).eq.'a-') then

        if(i.eq.1) bstring(1,1)=test_string
        if(i.gt.1) read(5,2001) bstring(i,1)

                                   else

        if(test_string(1:5).eq.'aux1a') then

          if(nt(i).eq.00) bstring(i,1)='a-dummy'
          if(nt(i).eq.01) bstring(i,1)='a-hydrogen (3,1;3,1)'
          if(nt(i).eq.02) bstring(i,1)='a-helium (3,1;3,1)'
          if(nt(i).eq.03) bstring(i,1)='a-lithium (4,3;4,3)'
          if(nt(i).eq.04) bstring(i,1)='a-beryllium (4,3;4,3)'
          if(nt(i).eq.05) bstring(i,1)='a-boron (4,4;4,4)'
          if(nt(i).eq.06) bstring(i,1)='a-carbon (4,4;4,4)'
          if(nt(i).eq.07) bstring(i,1)='a-nitrogen (4,4;4,4)'
          if(nt(i).eq.08) bstring(i,1)='a-oxygen (4,4;4,4)'
          if(nt(i).eq.09) bstring(i,1)='a-fluorine (4,4;4,4)'
          if(nt(i).eq.10) bstring(i,1)='a-neon (4,4;4,4)'
          if(nt(i).eq.11) bstring(i,1)='a-sodium (5,4;5,4)'
          if(nt(i).eq.12) bstring(i,1)='a-magnesium (5,4;5,4)'
          if(nt(i).eq.13) bstring(i,1)='a-aluminum (5,4;5,4)'
          if(nt(i).eq.14) bstring(i,1)='a-silicon (5,4;5,4)'
          if(nt(i).eq.15) bstring(i,1)='a-phosphorus (5,4;5,4)'
          if(nt(i).eq.16) bstring(i,1)='a-sulfur (5,4;5,4)'
          if(nt(i).eq.17) bstring(i,1)='a-chlorine (5,4;5,4)'
          if(nt(i).eq.18) bstring(i,1)='a-argon (5,4;5,4)'
          if(nt(i).eq.19) bstring(i,1)='a-potassium (5,5;5,5)'
          if(nt(i).eq.20) bstring(i,1)='a-calcium (5,5;5,5)'
          if(nt(i).eq.21) bstring(i,1)='a-scandium (5,5;5,5)'
          if(nt(i).eq.22) bstring(i,1)='a-titanium (5,5;5,5)'
          if(nt(i).eq.23) bstring(i,1)='a-vanadium (5,5;5,5)'
          if(nt(i).eq.24) bstring(i,1)='a-chromium (5,5;5,5)'
          if(nt(i).eq.25) bstring(i,1)='a-manganese (5,5;5,5)'
          if(nt(i).eq.26) bstring(i,1)='a-iron (5,5;5,5)'
          if(nt(i).eq.27) bstring(i,1)='a-cobalt (5,5;5,5)'
          if(nt(i).eq.28) bstring(i,1)='a-nickel (5,5;5,5)'
          if(nt(i).eq.29) bstring(i,1)='a-copper (5,5;5,5)'
          if(nt(i).eq.30) bstring(i,1)='a-zinc (5,5;5,5)'
          if(nt(i).eq.31) bstring(i,1)='a-gallium (5,5;5,5)'
          if(nt(i).eq.32) bstring(i,1)='a-germanium (5,5;5,5)'
          if(nt(i).eq.33) bstring(i,1)='a-arsenic (5,5;5,5)'
          if(nt(i).eq.34) bstring(i,1)='a-selenium (5,5;5,5)'
          if(nt(i).eq.35) bstring(i,1)='a-bromine (5,5;5,5)'
          if(nt(i).eq.36) bstring(i,1)='a-krypton (5,5;5,5)'
          if(nt(i).eq.37) bstring(i,1)='a-rubidium (5,5;5,5)'
          if(nt(i).eq.38) bstring(i,1)='a-strontium (5,5;5,5)'
          if(nt(i).eq.39) bstring(i,1)='a-yttrium (5,5;5,5)'
          if(nt(i).eq.40) bstring(i,1)='a-zirconium (5,5;5,5)'
          if(nt(i).eq.41) bstring(i,1)='a-niobium (5,5;5,5)'
          if(nt(i).eq.42) bstring(i,1)='a-molybdenum (5,5;5,5)'
          if(nt(i).eq.43) bstring(i,1)='a-technetium (5,5;5,5)'
          if(nt(i).eq.44) bstring(i,1)='a-ruthenium (5,5;5,5)'
          if(nt(i).eq.45) bstring(i,1)='a-rhodium (5,5;5,5)'
          if(nt(i).eq.46) bstring(i,1)='a-palladium (5,5;5,5)'
          if(nt(i).eq.47) bstring(i,1)='a-silver (5,5;5,5)'
          if(nt(i).eq.48) bstring(i,1)='a-cadmium (5,5;5,5)'
          if(nt(i).eq.49) bstring(i,1)='a-indium (5,5;5,5)'
          if(nt(i).eq.50) bstring(i,1)='a-tin (5,5;5,5)'
          if(nt(i).eq.51) bstring(i,1)='a-antimony (5,5;5,5)'
          if(nt(i).eq.52) bstring(i,1)='a-tellurium (5,5;5,5)'
          if(nt(i).eq.53) bstring(i,1)='a-iodine (5,5;5,5)'
          if(nt(i).eq.54) bstring(i,1)='a-xenon (5,5;5,5)'

                                        endif

                                   endif

      rewind(3)

 1002 read(3,2001) input_string
      if(input_string.ne.bstring(i,1)) goto 1002

      read(3,*) n

      if(n.ne.0) then
                       do 1003 j=1,n
                       read(3,*) alphacd(n1+j)
                       icdcfunc(n1+j)=i
                       iwkvec1(ncds+j)=n1+j
 1003                  coefscd(n1+j)=one/sqrt(sqrt(two)*
     &                                  ((pi/alphacd(n1+j))**fivhlf))
                 endif

      ncds=ncds+n
      n1=n1+n

      read(3,*) n

      if(n.ne.0) then
                       do 1004 j=1,n
                       read(3,*) alphacd(n1+j)
                       icdcfunc(n1+j)=i
                       iwkvec2(ncdspd+j)=n1+j
                       factor=sqrt(two)*((pi/alphacd(n1+j))**fivhlf)
                       coefscd(n1+j)=one/sqrt(factor)
                       factor=factor/(twelve*alphacd(n1+j))
                       coefpcd(n1+j)=one/sqrt(factor)
                       factor=factor*three/(twenty*alphacd(n1+j))
 1004                  coefdcd(n1+j)=one/sqrt(factor)
                 endif

      ncdspd=ncdspd+n
      n1=n1+n

      read(3,*) n

      if(n.ne.0) then
                       do 1005 j=1,n
                       read(3,*) alphaxc(n2+j)
                       ixccfunc(n2+j)=i
                       iwkvec3(nxcs+j)=n2+j
 1005                  coefsxc(n2+j)=(eight*alphaxc(n2+j)**3/pi**3)**
     &                                                          quarter
                 endif

      nxcs=nxcs+n
      n2=n2+n

      read(3,*) n

      if(n.ne.0) then
                       do 1006 j=1,n
                       read(3,*) alphaxc(n2+j)
                       ixccfunc(n2+j)=i
                       iwkvec4(nxcspd+j)=n2+j
                       coefsxc(n2+j)=(eight*alphaxc(n2+j)**3/pi**3)
     &                                                       **quarter
                       coefpxc(n2+j)=(onehund28*alphaxc(n2+j)**5/pi**3)
     &                                                       **quarter
 1006                  coefdxc(n2+j)=(twothou48*alphaxc(n2+j)**7/pi**3)
     &                                                       **quarter
                 endif

      nxcspd=nxcspd+n
      n2=n2+n

 1001 continue

      if(ncds.ne.0) then
                          do 1007 i=1,ncds
 1007                     icdlfunc(i)=iwkvec1(i)
                    endif

      if(nxcs.ne.0) then
                          do 1008 i=1,nxcs
 1008                     ixclfunc(i)=iwkvec3(i)
                    endif

      if(ncdspd.ne.0) then
                            do 1009 i=1,ncdspd
 1009                       icdlfunc(ncds+i)=iwkvec2(i)
                      endif

      if(nxcspd.ne.0) then
                            do 1010 i=1,nxcspd
 1010                       ixclfunc(nxcs+i)=iwkvec4(i)
                      endif

      ncdfuncs=ncds+10*ncdspd
      nxcfuncs=nxcs+10*nxcspd

      ncdspd=ncds+ncdspd
      nxcspd=nxcs+nxcspd

      nconts=0
      ncontp=0
      ncontd=0

      m=0
      n=0

      read(5,2001) test_string

      do 1011 i=1,ncenters

      if(test_string(1:2).eq.'o-') then

                          if(i.eq.1) bstring(1,2)=test_string
                          if(i.gt.1) read(5,2001) bstring(i,2)

                                   else

      if(test_string(1:7).eq.'tzp/dzp') then

       if(nt(i).eq.00) bstring(i,2)='o-dummy'
       if(nt(i).eq.01) bstring(i,2)='o-hydrogen (41/1)'
       if(nt(i).eq.02) bstring(i,2)='o-helium (51/1)'
       if(nt(i).eq.03) bstring(i,2)='o-lithium (7111/1/1)'
       if(nt(i).eq.04) bstring(i,2)='o-beryllium (7111/1/1)'
       if(nt(i).eq.05) bstring(i,2)='o-boron (7111/411/1)'
       if(nt(i).eq.06) bstring(i,2)='o-carbon (7111/411/1)'
       if(nt(i).eq.07) bstring(i,2)='o-nitrogen (7111/411/1)'
       if(nt(i).eq.08) bstring(i,2)='o-oxygen (7111/411/1)'
       if(nt(i).eq.09) bstring(i,2)='o-fluorine (7111/411/1)'
       if(nt(i).eq.10) bstring(i,2)='o-neon (7111/411/1)'
       if(nt(i).eq.11) bstring(i,2)='o-sodium (63111/411/1)'
       if(nt(i).eq.12) bstring(i,2)='o-magnesium (63111/411/1)'
       if(nt(i).eq.13) bstring(i,2)='o-aluminum (73111/6111/1)'
       if(nt(i).eq.14) bstring(i,2)='o-silicon (73111/6111/1)'
       if(nt(i).eq.15) bstring(i,2)='o-phosphorus (73111/6111/1)'
       if(nt(i).eq.16) bstring(i,2)='o-sulfur (73111/6111/1)'
       if(nt(i).eq.17) bstring(i,2)='o-chlorine (73111/6111/1)'
       if(nt(i).eq.18) bstring(i,2)='o-argon (73111/6111/1)'
       if(nt(i).eq.19) bstring(i,2)='o-potassium (633111/5211/1)'
       if(nt(i).eq.20) bstring(i,2)='o-calcium (633111/5211/1)'
       if(nt(i).eq.21) bstring(i,2)='o-scandium (633111/5211/311)'
       if(nt(i).eq.22) bstring(i,2)='o-titanium (633111/5211/311)'
       if(nt(i).eq.23) bstring(i,2)='o-vanadium (633111/5211/311)'
       if(nt(i).eq.24) bstring(i,2)='o-chromium (633111/5211/311)'
       if(nt(i).eq.25) bstring(i,2)='o-manganese (633111/5211/311)'
       if(nt(i).eq.26) bstring(i,2)='o-iron (633111/5211/311)'
       if(nt(i).eq.27) bstring(i,2)='o-cobalt (633111/5211/311)'
       if(nt(i).eq.28) bstring(i,2)='o-nickel (633111/5211/311)'
       if(nt(i).eq.29) bstring(i,2)='o-copper (633111/5211/311)'
       if(nt(i).eq.30) bstring(i,2)='o-zinc (633111/5211/311)'
       if(nt(i).eq.31) bstring(i,2)='o-gallium (633111/53111/41)'
       if(nt(i).eq.32) bstring(i,2)='o-germanium (633111/53111/41)'
       if(nt(i).eq.33) bstring(i,2)='o-arsenic (633111/53111/41)'
       if(nt(i).eq.34) bstring(i,2)='o-selenium (633111/53111/41)'
       if(nt(i).eq.35) bstring(i,2)='o-bromine (633111/53111/41)'
       if(nt(i).eq.36) bstring(i,2)='o-krypton (633111/53111/41)'
       if(nt(i).eq.37) bstring(i,2)='o-rubidium (6333111/53211/51)'
       if(nt(i).eq.38) bstring(i,2)='o-strontium (6333111/53211/51)'
       if(nt(i).eq.39) bstring(i,2)='o-yttrium (6333111/53211/5211)'
       if(nt(i).eq.40) bstring(i,2)='o-zirconium (6333111/53211/5211)'
       if(nt(i).eq.41) bstring(i,2)='o-niobium (6333111/53211/5211)'
       if(nt(i).eq.42) bstring(i,2)='o-molybdenum (6333111/53211/5211)'
       if(nt(i).eq.43) bstring(i,2)='o-technetium (6333111/53211/5211)'
       if(nt(i).eq.44) bstring(i,2)='o-ruthenium (6333111/53211/5211)'
       if(nt(i).eq.45) bstring(i,2)='o-rhodium (6333111/53211/5211)'
       if(nt(i).eq.46) bstring(i,2)='o-palladium (6333111/53211/5211)'
       if(nt(i).eq.47) bstring(i,2)='o-silver (6333111/53211/5211)'
       if(nt(i).eq.48) bstring(i,2)='o-cadmium (6333111/53211/5211)'
       if(nt(i).eq.49) bstring(i,2)='o-indium (6333111/533111/531)'
       if(nt(i).eq.50) bstring(i,2)='o-tin (6333111/533111/531)'
       if(nt(i).eq.51) bstring(i,2)='o-antimony (6333111/533111/531)'
       if(nt(i).eq.52) bstring(i,2)='o-tellurium (6333111/533111/531)'
       if(nt(i).eq.53) bstring(i,2)='o-iodine (6333111/533111/531)'
       if(nt(i).eq.54) bstring(i,2)='o-xenon (6333111/533111/531)'

                                        endif

                                 endif

      rewind(3)

 1012 read(3,2001) input_string
      if(input_string.ne.bstring(i,2)) goto 1012

      read(3,*) nfuncs,nfuncp,nfuncd

      do 1013 j=1,nfuncs

      read(3,*) nprimitives
      nshels(nconts+j)=m+j
      ngaussians(m+j)=nprimitives
      icfunc(m+j)=i
      ilfunc(m+j)=n+1

      do 1014 k=1,nprimitives
      read(3,*) alpha(n+k),coeff(n+k)
 1014 coeff(n+k)=coeff(n+k)*(eight*alpha(n+k)**3/pi**3)**quarter

      dn=zero

      do 1015 k=1,nprimitives
      do 1015 l=1,nprimitives
 1015 dn=dn+coeff(n+k)*coeff(n+l)*(pi/(alpha(n+k)+alpha(n+l)))**thrhlf

      do 1016 k=1,nprimitives
 1016 coeff(n+k)=coeff(n+k)/sqrt(dn)

 1013 n=n+nprimitives

      i_add=0
      if((idiffuse.eq.1).and.(nt(i).gt.2)) i_add=1
      if(idiffuse.eq.2) i_add=1

      if(i_add.eq.1) then
                           nshels(nconts+nfuncs+1)=m+nfuncs+1
                           ngaussians(m+nfuncs+1)=1
                           icfunc(m+nfuncs+1)=i
                           ilfunc(m+nfuncs+1)=n+1

                           alpha(n+1)=alpha(n)/three
                           coeff(n+1)=coeff(n)
                           coeff(n+1)=coeff(n+1)
     &                               *(eight*alpha(n+1)**3/pi**3)
     &                               **quarter

                           dn=zero

                           dn=dn+coeff(n+1)*coeff(n+1)
     &                          *(pi/(alpha(n+1)+alpha(n+1)))**thrhlf

                           coeff(n+1)=coeff(n+1)/sqrt(dn)

                           nfuncs=nfuncs+1

                           n=n+1
                     endif

      m=m+nfuncs

      do 1018 j=1,nfuncp

      read(3,*) nprimitives
      nshelp(ncontp+j)=m+j
      ngaussians(m+j)=nprimitives
      icfunc(m+j)=i
      ilfunc(m+j)=n+1

      do 1019 k=1,nprimitives
      read(3,*) alpha(n+k),coeff(n+k)
 1019 coeff(n+k)=coeff(n+k)*(onehund28*alpha(n+k)**5/pi**3)**quarter

      dn=zero

      do 1020 k=1,nprimitives
      do 1020 l=1,nprimitives
      factor=coeff(n+k)*coeff(n+l)*(pi/(alpha(n+k)+alpha(n+l)))**thrhlf
 1020 dn=dn+factor/(two*(alpha(n+k)+alpha(n+l)))

      do 1021 k=1,nprimitives
 1021 coeff(n+k)=coeff(n+k)/sqrt(dn)

 1018 n=n+nprimitives

      i_add=0
      if((idiffuse.eq.1).and.(nt(i).gt.2)) i_add=1
      if(idiffuse.eq.2) i_add=1

      if(i_add.eq.1) then
                           nshelp(ncontp+nfuncp+1)=m+nfuncp+1
                           ngaussians(m+nfuncp+1)=1
                           icfunc(m+nfuncp+1)=i
                           ilfunc(m+nfuncp+1)=n+1

                           alpha(n+1)=alpha(n)/three
                           coeff(n+1)=coeff(n)
                           coeff(n+1)=coeff(n+1)
     &                               *(onehund28*alpha(n+1)**5/pi**3)
     &                               **quarter

                           dn=zero

                           factor=coeff(n+1)*coeff(n+1)
     &                           *(pi/(alpha(n+1)+alpha(n+1)))**thrhlf
                           dn=dn+factor/(two*(alpha(n+1)+alpha(n+1)))

                           coeff(n+1)=coeff(n+1)/sqrt(dn)

                           nfuncp=nfuncp+1

                           n=n+1
                     endif

      m=m+nfuncp

      do 1022 j=1,nfuncd

      read(3,*) nprimitives
      nsheld(ncontd+j)=m+j
      ngaussians(m+j)=nprimitives
      icfunc(m+j)=i
      ilfunc(m+j)=n+1

      do 1023 k=1,nprimitives
      read(3,*) alpha(n+k),coeff(n+k)
 1023 coeff(n+k)=coeff(n+k)*(twothou48*alpha(n+k)**7/pi**3)**quarter

      dn=zero

      do 1024 k=1,nprimitives
      do 1024 l=1,nprimitives
      factor=coeff(n+k)*coeff(n+l)*(pi/(alpha(n+k)+alpha(n+l)))**thrhlf
 1024 dn=dn+factor/((two*(alpha(n+k)+alpha(n+l)))**2)

      do 1025 k=1,nprimitives
 1025 coeff(n+k)=coeff(n+k)/sqrt(dn)

 1022 n=n+nprimitives

      i_add=0
      if((idiffuse.ne.0).and.(nt(i).gt.18)) i_add=1

      if(i_add.eq.1) then
                           nsheld(ncontd+nfuncd+1)=m+nfuncd+1
                           ngaussians(m+nfuncd+1)=1
                           icfunc(m+nfuncd+1)=i
                           ilfunc(m+nfuncd+1)=n+1

                           alpha(n+1)=alpha(n)/three
                           coeff(n+1)=coeff(n)
                           coeff(n+1)=coeff(n+1)
     &                               *(twothou48*alpha(n+1)**7/pi**3)
     &                               **quarter

                           dn=zero

                           factor=coeff(n+1)*coeff(n+1)
     &                           *(pi/(alpha(n+1)+alpha(n+1)))**thrhlf
                           dn=dn+factor
     &                          /((two*(alpha(n+1)+alpha(n+1)))**2)

                           coeff(n+1)=coeff(n+1)/sqrt(dn)

                           nfuncd=nfuncd+1

                           n=n+1
                     endif

      m=m+nfuncd

      nconts=nconts+nfuncs
      ncontp=ncontp+nfuncp
      ncontd=ncontd+nfuncd

 1011 continue

      ncontractions=nconts+3*ncontp+6*ncontd

      ntemp=nxcfuncs
      if(ixcfit.eq.0) ntemp=0

      write(6,2002) ncontractions,ncdfuncs

      do 1026 i=1,nconts
 1026 mtloca(nshels(i))=i

      if(ncontp.ne.0) then
                            do 1027 i=1,ncontp
 1027                       mtloca(nshelp(i))=nconts+3*(i-1)+1
                      endif

      if(ncontd.ne.0) then
                            do 1028 i=1,ncontd
 1028                       mtloca(nsheld(i))=nconts+3*ncontp+6*(i-1)+1
                      endif

 2001 format(a30)
 2002 format(' number of orbital basis functions                =',i5,/,
     &       ' number of charge density fitting functions       =',i5,/)

      return
      end
