c
c     pa60 - polyatom properties program for s,p,d and f-orbitals
c
c      250 basis functions allowed
c
c     input modified according to the karlsruhe version of pitzer's
c     scf-program by c.m.ehrhardt
c     for this purpose the subroutines inpsym,getran,scfinp and comtrn
c     are taken from h.f. scfaefer's gradient-program
c
c     the subroutines rhftce,propa and prod were made faster (4.3.1984)
c
c     the program was modified to handle f-functions by p.maier and
c     c.m.ehrhardt     ( 1.3.1984 )
c
c     mulliken and roby population analysis included ( 12.7.1984 )
c
c     modified to calculate the kinetic energy and its relativistic
c     correction by s. brode in april,1984 : a subroutine opap4 was
c     added and the main program and the subroutines dintfm, pint,
c     propa and prnt were slightly modified
c
c     boys localization included by r.ahlrichs and h.schiffer
c     (1.3.1985)
c
c     ft02,ft04,ft07,ft09       scratch files
c     ft22                      localized mos for plopr
c     ft21                      electrostatic potential
c     ft23                      file for plot
c
c     input description:
c
c     input files:
c     ft24   integral information from the karlsruhe integral input
c            program
c     ft08   eigenvectors and occupation as written out by subroutine
c            outgrd of scf-program
c     ft05   input unit
c     ft28   optional file for atomic coordinates (compatible to ints)
c
c     input on unit ft05:
c
c     1.    title                         (format 12a6)
c
c     2.    isl1 , isl2 , isl4 , nsupp    (format 4i5)
c           isl1 = 1   full property calculation (default)
c           isl1 = 2   property only in ao-basis
c           isl2 = 1   no output of property matrix in ao-basis(default)
c           isl2 = 2   output of these integrals
c           isl2 = 3   output, and integrals are written on unit 3
c           isl4 = 1   renormalisation of wavefunctions and output
c           isl4 = 2   no renormalisation
c           isl4 = all other values   renormalisation but no output
c                      (default)
c           nsupp= 1   output of diagonalelements of property matrix
c                      in mo-basis
c           nsupp= 2   controll-output of integral and scf input data
c           nsupp= 3   output of mo vectors in ao basis
c           nsupp= 4   output of overlap population per mo
c           nsupp= 0   no output (default=0)
c
c     3.    m, ilnmc, icntr, (c(k),k=1,3) (format i5,6x,a4,i5,3f10.5)
c           step 3 can be repeated several times
c           population analysis or localization must always be the
c           last option
c           m specifies the property :
c           potential(1), electric field(2), field gradient(3),
c           dipole moment(4), quadrupole moment(5), third moment(6),
c           planar density(7), linear density(8), charge density(9),
c           overlap(10), second moment(11), diamagnetic shielding(12),
c           relativistic correction(13), mulliken and roby population
c           analysis(14), boys localization(15), plot(16)
c
c           ilnmc      name of the point at which property shall be
c                      calculated
c           icntr      number of atom at which property shall be
c                      calculated, if icntr.ne.0 you need not read in
c                      coordinates
c           c(k)       cartesian coordiantes of the point
c
c     4. special input for population analysis (m=14)
c
c           nmao, (mor(im),im=1,nmao)  (format i10, 20 i3)
c           nmao       number of maos for atom 1,2...,n
c           mor        numbers of selected maos
c           mor(1)= 0  automatic selection
c           nm, (iatnr(i),i=1,nm)  (format 20i3)
c           nm.ne.0    special multicenter term of order nm
c           iatnr(i)   atoms that form the multicenter term
c           (a blank line is required if nm=0)
c
c     5. special input for localization (m=15)
c           nmo  (format 5x,i5)
c           nmo =  -1  all occupied mo's shall be localized
c           nmo =   0  now the program wants to read a threshold
c                      for the selection of mo's, wich shall be
c                      localized
c           nmo = nmo  number of mo's, wich shall be localized
c                      now the program wants to read the corresponding
c                      mo numbers in ascending order
c either thr (5x,f10.4) for nmo = 0
c or     mon(i),i=1,nmo (5x,(5i5)) for nmo > 0
c        nop (5x,i5)
c           nop =   0  no further input needed for localization
c           nop = nop  number of operators ( 1,2 or 3 )
c                      now the program wants to read the corresponding
c                      operator numbers ( x=1, y=2, z=3 ) in canonical
c                      order
c           opn(i),i=1,nop (5x,3i5) for nop > 0
c
c
c                 important dimensions
c
c   number of atoms:                        nocmx  =      50
c   number of primitive basis functions  :  nbmx   =     720
c   number of contracted basis functions
c   and number of mos                    :  nomx   =     250
c   dimension of buff                    :  lbuff  = nomx*(nomx+1)/2
c   dimensions in common /big/           :  maxuv  = nomx**2
c   max number of orbitals to be local.  :  mloc   =     250
c   (in common eig and scfinf)
c   dimension of rnabcd                  :  n*(n-3)/2+m
c   with:            n=nocmx*(nocmx-3)/2+nocmx-1
c                    m=(nocmx-2)*(nocmx-5)/2+nocmx-3
c   max number of grids for 'x' direction:  nmaxx  =        50
c   max number of grids for 'y' direction:  nmaxy  =        50
c   max number of orbitals to be ploted  :  norbmx =        10
c
      implicit real*8(a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375,nrabcd=748379)
      parameter(nmaxx=50,nmaxy=50,norbmx=10,icomx=12)
      real*8 icod(12),ilab(12),ilbl(12)
      dimension eta(nbmx,25),fcntr(nbmx),
     1         ftype(nbmx),frocc(nomx),c(3),dint(3,10),
     2         type(20),val(10),tot(10),nstart(nomx),nstop(nomx),
     3         ncontc(nomx),buff(lbuff),a(39),title(150),
     4         cint(nomx,10),izahl(nbmx),itrnao(nbmx),idxxt(nbmx),
     5         rnabcd(nrabcd),contrc(nomx,icomx)
      dimension nocc(norbmx),cc(nmaxx*nmaxy,3),
     1         dd(nmaxx*nmaxy,3),dinx(nmaxx*nmaxy),
     2         cin(nmaxx*nmaxy,norbmx+1),
     3         v(nmaxx*nmaxy),ddsq(nmaxx*nmaxy)
      equivalence (buff(nocmx+1),rnabcd(1))
      external opaa0,opaa1,opaa2,opab1,opab2,opab3,opac1,opac2, opac3,
     1opad1,opab4,opaa3,opap4
      external potnuc,fldnuc,grdnuc,dipnuc,qudnuc,octnuc,nonuc,magnuc,
     1sldnuc
      integer centr,fcntr,ftype,type
      common/ioind/isl1, isl2
      common/muli/vlist(nocmx,4),ianr(nomx),centr(nocmx),itnr(nomx)
      common/holl/blabel(19),alabel(19),fmt(4),jl,lomao,iopen,nsym
      common/supp/nsupp
      common/ncnt/ilnmcm
      common/prptyp/m
      common /contrc/ chomo,ncont
      common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
c
      data nromx/3/,nfcol/25/,ntcol/20/,iwarn/0/,ir/8/
      data icod/6h1pa60 ,6h      ,6h      ,6hnyu pr,6hoperti,6hes pac,
     1          6hkage  ,6h      ,6h      ,6h      ,6h12/66 ,6hspdf  /
      data title/6h   1/r,6h    ex,6h   fxx,6h     x,6h   qxx,6h  oxxx,
     1  6h    dx,6h   dxy,6h   den,6h   ovl,6h    xx,6h fx*x*,6h qx*x*,
     2  6h  x*x*,6h<p**2>,                  6h      ,6h    ey,6h   fyy,
     3  6h     y,6h   qyy,6h  oyyy,6h    dy,6h   dxz,6h      ,6h      ,
     4  6h    yy,6h fy*y*,6h qy*y*,6h  y*y*,6h<p**4>,
     5  6h      ,6h    ez,6h   fzz,6h     z,6h   qzz,6h  ozzz,6h    dz,
     6  6h   dyz,6h      ,6h      ,6h    zz,6h fz*z*,6h qz*z*,6h  z*z*,
     7  6h sum  ,                  6h      ,6h      ,6h   fxy,6h      ,
     8  6h   qxy,6h  oxxy,6h      ,6h      ,6h      ,6h      ,6h    xy,
     9  6h fx*y*,6h qx*y*,6h  x*y*,6h      ,                  6h      ,
     1  6h      ,6h   fxz,6h      ,6h   qxz,6h  oxxz,6h      ,6h      ,
     2  6h      ,6h      ,6h    xz,6h fx*z*,6h qx*z*,6h  x*z*,6h      ,
     3                    6h      ,6h      ,6h   fyz,6h      ,6h   qyz,
     4  6h  oxyy,6h      ,6h      ,6h      ,6h      ,6h    yz,6h fy*z*,
     5  6h qy*z*,6h  y*z*,6h      ,
     6  4*6h      ,6h   rsq,6h  oyyz,9*6h      ,
     7           5*6h      ,6h  oxzz,9*6h      ,
     8           5*6h      ,6h  oyzz,9*6h      ,
     9           5*6h      ,6h  oxyz,9*6h      /
      data a/6hpotent,6helectr,6hfield ,6hdipole,6hquadru,6hthird ,
     1       6hplanar,6hline d,6hcharge,6hoverla,6hsecond,6hdiamag,
     &       6hrelati,
     2       6hial   ,6hic fie,6hgradie,6h momen,6hpole m,6hmoment,
     3       6h densi,6hensity,6h densi,6hp     ,6h momen,6hnetic ,
     &       6hvistic,
     4       6h      ,6hld    ,6hnt    ,6ht     ,6homent ,6h      ,
     5       6hty    ,6h      ,6hty    ,6h      ,6ht     ,6hshldg ,
     &       6h corr./
      data type/4hs   ,4hx   ,4hy   ,4hz   ,4hxx  ,4hyy  ,4hzz  ,
     1     4hxy  ,4hxz  ,4hyz  ,4hxxx ,4hyyy ,4hzzz ,4hxxy ,4hxxz ,
     2     4hyyx ,4hyyz ,4hzzx ,4hzzy ,4hxyz /
c
c**      call xuflow(0)
c        file shout short output
      open(unit=5,file='propin')
      open(6,file='propls')
      nshout=13
      open(unit=nshout,file='shout',status='unknown'
     &,access='sequential',form='formatted')
      rewind nshout
c        file vectgrd scf vectors
c        ir=8
      open(unit=ir,file='vectgrd',status='old',
     1access='sequential',form='formatted')
c        file plot: plot output
      iplot=23
      open(unit=iplot,file='plot',status='unknown',form='formatted')
      rewind iplot
c
      maxuv=nomx*nomx
      lomao=22
   76 do 20 i=1,12
   20 ilab(i) = icod(i)
      do 22 i=1,nomx
   22 frocc(i) = 0.0d0
      read(5,101,end=999) (ilbl(i), i=1,12)
      write(6,101) (ilab(i), i=1,12)
      write(6,103)
      write(6,104) (ilbl(i), i=1,12)
      read(5,100) isl1,isl2,isl4,nsupp
      if(isl1.eq.0) isl1=1
      if(isl2.eq.0) isl2=1
      if(isl4.eq.0) isl4=3
      if(isl2.eq.3) rewind 3
      call input(eta,nfcol,fcntr,ftype,n,nbmx,vlist,centr,non,nocmx,
     1 ntcol,type,ilab,ilbl,nomx,nstart,
     2 nstop,contrc,ncontc,buff,izahl,nbfao,icomx,ianr,itnr)
      call nrml(eta,nbmx,nfcol,nomx,contrc,nstart,nstop,fcntr,ftype,
     1isl4,c,opad1)
c
c     ivmax.ne.0 : transformation of mo-vectors in ao basis has already
c                  been done
      rewind ir
      read (ir,940)(alabel(i1),i1=1,19)
      read(ir,940)(blabel(i1),i1=1,19)
      read(ir,950)iopen,nsym,ivmax
      read(ir,940)(fmt(i1),i1=1,4)
      if(ivmax.eq.0) goto 170
      if(ivmax.gt.maxuv) goto 1500
      read(ir,fmt) (y(i),i=1,ivmax)
      norb=iopen
      do 180 nn=1,norb
      frocc(nn)=2.0d0
  180 chomo=chomo-2.0d0
      write(6,970) chomo
      goto 25
c
  170 call getran(nbfao,frocc,norb,itrnao,idxxt)
   25 read(5,102,end=999)m,ilnmcm,icntr,(c(k),k=1,3)
      if(iwarn.eq.1) goto 998
      if(m.eq.14.or.m.eq.15) iwarn=1
      if(icntr.eq.0) go to 800
      do 801 i=1,3
  801 c(i)=vlist(icntr,i)
  800 if(m.eq.0) go to 76
      go to (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),m
    1 call dintfm(opaa0, potnuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
    2 call dintfm(opaa1, fldnuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
    3 call dintfm(opaa2, grdnuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
    4 call dintfm(opab1, dipnuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
    5 call dintfm(opab2, qudnuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
    6 call dintfm(opab3, octnuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
    7 call dintfm(opac1,  nonuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
    8 call dintfm(opac2,  nonuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
    9 call dintfm(opac3,  nonuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
   10 call dintfm(opad1,  nonuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
   11 call dintfm(opab4, magnuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
   12 call dintfm(opaa3,sldnuc, dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
   13 call dintfm(opap4, nonuc, dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
   14 call dintfm(opad1,  nonuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
   15 call dintfm(opab1, nonuc,dint,    vlist,c,n,eta,val,nt
     1,norb,frocc,tot,non,nromx,ilbl,a,title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      go to 25
   16 call plotin(nmaxx,nmaxy,kpmax,cc,norbpl,mdens,nocc,norbmx,iplot,
     1 iat)
      if(iat.gt.0) call molcor(iplot)
      call orbval(eta,nbmx,norbmx,nomx,nstart,nstop,contrc,ftype,
     1frocc,cc,dd,dinx,cin,v,ddsq,kpmax,norbpl,mdens,nocc,y,iplot)
      go to 25
  999 stop
 1500 write(6,960) ivmax,maxuv
      stop
  100 format(6i5)
  101 format(12a6)
  102 format(i5,6x,a4,i5,3f10.0)
  103 format(/)
  104 format(1x,12a6)
  105 format(2x,10a4)
  940 format(19a4)
  950 format(20i4)
  960 format(//' dimension of y in main program  too small !  ivmax= ',
     &i6,'  maxuv= ',i6//)
  970 format(/' charge of the molecule ',f10.5/)
c
  998 write(6,106)
  106 format(1x,'option 14 (population analysis) or option 15 (localizat
     1ion) must be the last one'/,' ***** calculation stopped *****')
      stop
      end
c
c
      subroutine input( eta, nfcol, fcntr, ftype, n, nbmx, vlist, centr,
     1 non,nocmx,ntcol,type,ilab,ilbl,nomx,
     2 nstart,nstop,contrc,ncontc,buff,izahl,nbfao,icomx,ianr,
     3 itnr)
      implicit real*8(a-h,o-z)
c
c  parameter statement as in integral program
      parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
     &          lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
     &          lpru=18000,lsu=50,lstu=8,lconu=12,
     &          lcu=8,lccu=182)
c
c     this subroutine prepares the integral input from the karlsruhe
c     input program for polyatom.
c     one part is taken from routine convrt of schaefer's gradient
c     program
c
      real*8    ilab(12), ilbl(12)
      integer fcntr, ftype, centr, type
      dimension eta(nbmx,nfcol), vlist(nocmx,4),ianr(1),itnr(1)
      dimension fcntr(1), ftype(1), centr(1), type(1),contrc(nomx,8)
      dimension nstart(1),nstop(1),ncontc(1),buff(1),frocc(1),izahl(1)
      common /contrc/ chomo,ncont
      common /ioind/isl1,isl2
      common/infoa/natoms,mul,nbf,nij,nelec,nalpha,nbeta
      common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
     1mcxu,ngcs,mdum(45),mccu,mblu,ns,ng,nsf,nst,ncons
      common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
     1  ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
     2  ccoef(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
     3  z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
     4  la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
     5  mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu)
c
      msu=lsu
      mconsu=lconsu
      mconu=lconu
      mcu=lcu
      msfu=lsfu
      mgu=lgu
      mnsfu=lnsfu
      mpru=lpru
      mccu=lccu
      kaords=laords
      mgcsu=lgcsu
      mstu=lstu
      mru=lru
      mcsu=lcsu
      mctu=lctu
      msfru=lsfru
      call inpsym(mcu,msu,mru,mctu,mcsu,icomx,mgcsu)
      n=0
      nbf=0
      ispf=0
      natoms=0
      nr=0
      do 500 is=1,ns
      icu=nc(is)
      ifu=nf(is)
      do 400 ic=1,icu
      isf=ispf
      natoms=natoms+1
      if(natoms.le.nocmx) goto 504
      write(6,1200)
 1200 format(1x,'too many atoms in subroutine input')
      stop
  504 centr(natoms)=mtype(is)
      vlist(natoms,1)=x(ic,is)
      vlist(natoms,2)=yy(ic,is)
      vlist(natoms,3)=z(ic,is)
      vlist(natoms,4)=chg(is)
      chomo=chomo+chg(is)
      nr=nr+1
      do 400 if=1,ifu
      isf=isf+1
      icons=mcons(isf)
      ncns=ncon(icons)
      jtyp=lmnp1(icons)
      if(jtyp.le.4) goto 11
      write(6,1000) jtyp
 1000 format(1x,' jtyp = ',i3,4x,' only s,p,d and f orbitals are allowed
     1 ',/,'****** calculation stopped ******')
      stop
   11 continue
      jtypm=jtyp*(jtyp+1)/2
      jtyp = (jtypm*(jtyp-1))/3 +1
      do 400 i=1,jtypm
      nbf=nbf+1
      if(nbf.lt.nomx) goto 503
      write(6,1100) nbf,nomx
 1100 format(1x,' dimension too small in input, nbf =',i5,3x,
     1' nomx =',i5)
      stop
  503 ncontc(nbf)=ncns
      itnr(nbf)=jtyp-1+i
      ianr(nbf)=nr
      do 400 icon=1,ncns
      n=n+1
      if(n.le.nbmx) goto 505
      write(6,1300)
 1300 format(1x,'too many primitive functions in subroutine input')
      stop
  505 izahl(n)=nr
      fcntr(n)=mtype(is)
      ftype(n)=type(jtyp-1+i)
      eta(n,4)=zet(icon,icons)
      contrc(nbf,icon)=ccoef(icon,icons)
  400 continue
      ispf=ispf+ifu
  500 continue
      non=natoms
      ncont=nbf
      nbfao=nbf
      call rdetac(eta,nfcol,fcntr,ftype,n,nbmx,vlist,centr,non,
     1 nocmx,ntcol,type,nstart,nstop,contrc,ncontc,nomx,
     2 izahl)
      return
      end
c
c
      subroutine rdetac(eta,nfcol,cntr,typ,n,nfmx,vlist,centre,non,
     1ncmx,ntcol,type,nstart,nstop,contrc,ncontc,nomx,izahl)
c
c     polyatom inputroutine for basis functions
c
      implicit real*8(a-h,o-z)
      common /contrc/ chomo,ncont
      dimension ncontc(1),nstart(1),nstop(1),contrc(nomx,8)
      dimension eta(nfmx,nfcol),vlist(ncmx,4)
      dimension cntr(1),typ(1),centre(1),type(1),izahl(1)
      integer cntr,typ,centre,blank,type
      write(6,105)
c.....center specifications
      write(6,104)
      write(6,106)
      do 1 i=1,non
      write(6,107) centre(i), (vlist(i,j), j=1,4)
    1  continue
      if(ncont.ne.n) go to 202
      do 201 i=1,n
  201 ncontc(i) = 1
  202 nstart(1)=1
      nstop(1)=ncontc(1)
      do 31 i=2,ncont
      nstart(i)=nstart(i-1)+ncontc(i-1)
   31 nstop(i)=nstop(i-1)+ncontc(i)
c.....set up coordinates of function in eta
      icon=1
      jcon=0
      do 3 i=1,n
      if(i-nstop(icon)) 60,60,61
   60 jcon=jcon+1
      go to 62
   61 icon=icon+1
      jcon=1
   62 j=izahl(i)
      do 6 m=1,3
    6 eta(i,m)=vlist(j,m)
c.....set up unnormalized hybridization coeffs of function in eta
      do 7 j=1,ntcol
      if(noteqv(type(j),typ(i))) 7, 8, 7
    7 continue
      write(6,2000) type(j),typ(i)
 2000 format('0***** types not equal  ',2a4)
      stop
    8 continue
      do 9 m=1,ntcol
    9 eta(i,m+5)=0.0d0
      do 2 m=1,ntcol
      if(typ(i).ne.type(m)) goto 2
      eta(i,5)=m
      eta(i,m+5)=1.0
      goto 3
    2 continue
    3 continue
      return
  101 format(2i3)
 1020 format(a4,6x,4f12.0)
 1030 format(36i2)
  103 format(a4,6x,a4,i3,3x,2f12.0)
  104 format(/ ,1x,15hnuclear centres)
  105 format(///,1x,32hgaussian function specifications)
  106 format(/ ,1x,6hcentre,18x,12hco-ordinates,17x,6hcharge,//)
  107 format(1x,a6,6x,3f12.8,6x,f4.1)
  110 format(///,1x,8hfunction,3x,11hcontraction,3x,6hcentre,3x,4htype,
     111x,8hexponent,3x,11hcoefficient)
  111 format(1x,i5,8x,i5,8x,i3,a3,3x,a6,3x,2f13.7)
  114 format(//,1x,18hadditional centres,//)
  115 format(//,1x,6hcentre,18x,12hco-ordinates//)
      end
c
c
      subroutine nrml(eta,nfmx,nfcol,nomx,c,nstart,nstop,nc,nt,isl4,d,
     1ovl)
      implicit real*8 (a-h,o-z)
      dimension a(20),eta(nfmx,nfcol),c(nomx,8),nstart(1),nstop(1),
     1nc(1),nt(1),d(3)
      data  cons/15.7496099d0/,c1/1.0d0/,zero/0.d0/,c3/3.0d0/,c4/4.0d0/,
     1 delt/1.0d-6/,c15/15.0d0/
      common /contrc/chomo,ncont
c     normalize primitive functions.
      do 5 i=1,ncont
      is=nstart(i)
      if=nstop(i)
      do 4 ii=is,if
      z=eta(ii,4)*c4
      do 1 m=1,20
    1 a(m)=eta(ii,m+5)
      s=cons*(a(1)+(a(2)+a(3)+a(4)+(a(8)+a(9)+a(10)+c3*(a(5)+a(6)+a(7))+
     1  (a(20)+c3*(a(14)+a(15)+a(16)+a(17)+a(18)+a(19))+c15*(a(11)+a(12)
     2 +a(13))) /z)/z)/z)/(z*sqrt(z))
      s=c1/sqrt(s)
      do 3 m=6,nfcol
    3 eta(ii,m)=eta(ii,m)*s
    4 continue
    5 continue
      if(isl4.eq.2) return
c     renormalize contraction coefficients
      d(1)=zero
      d(2)=zero
      d(3)=zero
      if(isl4.ne.1) go to 2
      write(6,101)
      write(6,110)
    2 do 9 i=1,ncont
      is=nstart(i)
      if=nstop(i)
      ip=0
      s=zero
      do 7 ii=is,if
      ip=ip+1
      jp=0
      do 7 jj=is,if
      call propa(ovl,ii,jj,eta,a,d,1,nfmx)
      jp=jp+1
      s=s+c(i,ip)*c(i,jp)*a(1)
    7 continue
      s=c1/sqrt(s)
      ip=0
      do 8 ii=is,if
      ip=ip+1
      c(i,ip)=c(i,ip)*s
      if(abs(c(i,ip)+c1).le.delt) c(i,ip)=c1
      if(isl4.ne.1) go to 8
      write(6,111)ii,i,nc(ii),nt(ii),eta(ii,4),c(i,ip)
    8 continue
    9 continue
  110 format(//,1x,8hfunction,3x,11hcontraction,3x,6hcentre,3x,4htype,
     111x,8hexponent,3x,11hcoefficient)
  111 format(1x,i5,8x,i5,8x,a6,3x,a6,3x,2f13.7)
  101 format(///,10x,15hrenormalization)
      return
      end
c
c
      subroutine dintfm(opa, nuc, dint, vlist, c, n, eta,
     1 val, nt, norb, frocc, tot, noc, nro, ilbl, a, title,
     2 nstart,nstop,contrc,cint,buff,rnabcd)
      implicit real*8(a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
      common /contrc/ chomo,ncont
      common /prptyp/ m
      common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
      dimension dint(nro,*), vlist(nocmx,4), c(3), eta(nbmx,16), val(1),
     1          frocc(1), tot(1)
      dimension nstart(1),nstop(1),contrc(nomx,8),cint(nomx,1),
     1 buff(lbuff),rnabcd(1)
      real*8 ilbl
      common /ioind/isl1,isl2
      dimension ilbl(1), a(1), title(1), dinx(310)
      dimension lin(84),min(84),nin(84)
      data pi/3.14159265359d0/, cl/137.03604d0/
      data lin/0,1,0,0,2,0,0,1,1,0,3,0,0,2,2,1,0,1,0,1,4,0,0,3,3,1,0,1,0
     1 ,2,2,0,2,1,1,5,0,0,3,3,2,2,0,0,4,4,1,0,0,1,1,3,1,2,2,1,6,0,0,3,
     2 3,0,5,5,1,0,0,1,4,4,2,0,2,0,3,3,1,2,2,1,4,1,1,2/,min/0,0,1,0,0,
     3 2,0,1,0,1,0,3,0,1,0,2,2,0,1,1,0,4,0,1,0,3,3,0,1,2,0,2,1,2,1,0,5,
     4 0,2,0,3,0,3,2,1,0,4,4,1,0,1,1,3,2,1,2,0,6,0,3,0,3,1,0,0,1,5,5,2,0
     5,0,2,4,4,2,1,3,1,3,2,1,4,1,2/,nin/0,0,0,1,0,0,2,0,1,1,0,0,3,0,1,0,
     6 1,2,2,1,0,0,4,0,1,0,1,3,3,0,2,2,1,1,2,0,0,5,0,2,0,3,2,3,0,1,0,1,
     7 4,4,3,1,1,1,2,2,0,0,6,0,3,3,0,1,5,5,1,0,0,2,4,4,0,2,1,2,2,3,1,3,
     8 1,1,4,2/
c
c
      call dumdum(nuc,noc,vlist,c,dint,nro,nocmx,nt)
      if(m.eq.15) nt=3
      if(m.eq.13) nt=3
      call pint(opa,       n, eta, c, val, nt, nro, nbmx, a, title,
     1 ilbl,nstart,nstop,contrc,nomx,buff,lbuff)
      if (m.ne.13) go to 7
c
c     calculation of the darwin-term
c
      valval=0.0d0
c     + loop over the centers
      do 80 ic=1,noc
      c(1)=vlist(ic,1)
      c(2)=vlist(ic,2)
      c(3)=vlist(ic,3)
      ipp=0
c     + + loop over the c-gto's
      do 40 iz=1,ncont
      dinx(iz)=0.d0
      istart=nstart(iz)
      istop=nstop(iz)
      iprim=0
c     + + + loop over the gto's
      do 40 ii=istart,istop
      ipp=ipp+1
      iprim=iprim+1
      ieta=eta(ii,5)+0.00001
      li=lin(ieta)
      mi=min(ieta)
      ni=nin(ieta)
      xc=c(1)-eta(ii,1)
      yc=c(2)-eta(ii,2)
      zc=c(3)-eta(ii,3)
      expo=-eta(ii,4)*(xc*xc+yc*yc+zc*zc)
      din=1.d0
      if(li.ne.0) din=xc**li
      if(mi.ne.0) din=din*yc**mi
      if(ni.ne.0) din=din*zc**ni
      din=din*exp(expo)
      ieta=ieta+5
      dinx(iz)=dinx(iz)+din*contrc(iz,iprim)*eta(ipp,ieta)
   40 continue
      vval=0.d0
      ii=0
c     + + loop over the orbitals
      do 70 mi=1,norb
      dens=0.d0
c     + + + loop over the c-gtos
      do 75 mj=1,ncont
      ii=ii+1
      dens=dens+dinx(mj)*y(ii)
   75 continue
      vval=vval+frocc(mi)*dens*dens
   70 continue
      valval=valval+vlist(ic,4)*vval
   80 continue
      dint(1,3)=valval*0.5d0*pi/cl**2
c
c
    7 if (m.ne.14) go to 6
c
c     mulliken  and roby population analysis
c
      call cntrcm(cint,norb,noc,buff,frocc,rnabcd)
      return
c
    6 go to (4,5), isl1
    4 call output(dint,norb,n,nt,nro,frocc,tot,c,ilbl,a,title,
     1 cint,buff)
    5 return
      end
c
c
      subroutine pint(opa,n,eta,c,val,nt,nro,nbmx,a,title,
     1 ilbl,nstart,nstop,contrc,nomx,buff,lbuff)
c
c     calculation of integrals in ao-basis. data are written on
c     file ft002
c
      implicit real*8(a-h,o-z)
      common /contrc/ chomo,ncont
      dimension nstart(1),nstop(1),contrc(nomx,8)
      dimension a(13,3), title(14,10)
      dimension c(1),eta(nbmx,16),val(1),dinx(10),buff(lbuff)
      real*8 ilbl(1)
      common/ioind/isl1,isl2
      common/prptyp/mm
      ipr=ncont*(ncont+1)/2
      if(ipr.gt.lbuff.and.mm.eq.14) goto 99
      nu = 3
      rewind 2
      ix=0
      do 50 i=1,ncont
      do 50 j=1,i
      nu = nu + 1
      istart = nstart(i)
      istop = nstop(i)
      jstart = nstart(j)
      jstop = nstop(j)
      iprim = 0
      do 2  m=1,nt
    2 dinx(m)=0.0d0
      do 40 ii=istart,istop
      iprim = iprim + 1
      jprim = 0
      do 40 jj=jstart,jstop
      jprim = jprim + 1
      call propa(opa,ii,jj,eta,val,c,nt,nbmx)
      do 30 m=1,nt
   30 dinx(m)    = contrc(i,iprim)*contrc(j,jprim)*val(m) + dinx(m)
   40 continue
      if(mm.ne.14) goto 55
      ix=ix+1
      buff(ix)=dinx(1)
      goto 50
   55 do 60 m=1,nt
      ix=ix+1
      buff(ix)=dinx(m)
      if(ix.lt.lbuff)go to 60
      ix=0
      write(2)buff
   60 continue
   50 continue
      if(mm.eq.14) return
      write(2)buff
      if(isl2.eq.1)return
      if(isl2.eq.3) go to 1
      write(6,19) (ilbl(i), i=1,12)
      write(6,20) (a(mm,i), i=1,3)
      write(6,21)
      rewind 2
      ix=lbuff
      do 6 i=1,ncont
      do 7 j=1,i
      write(6,21)
      do 70 m=1,nt
      ix=ix+1
      if(ix.le.lbuff)go to 70
      ix=1
      read(2)buff
   70 write(6,10)i,j,buff(ix)
    7 continue
    6 continue
      return
    1 rewind 2
      ix=lbuff
      do 66 i=1,ncont
      do 77 j=1,i
      do 770 m=1,nt
      ix=ix+1
      if(ix.le.lbuff) go to 770
      ix=1
      read(2) buff
      write(3) buff
  770 continue
   77 continue
   66 continue
      return
c
   99 write(6,7000) ipr,lbuff
 7000 format(1x,'dimension too small in subroutine pint  ipr =',i8,
     1 ' lbuff = ',i8,/,' population analysis is not possible')
      stop
c
   19 format(1h0,/ ,1x,12a6)
   20 format(    //,1x,34hintegrals over basis functions for,5x,3a6//)
   21 format(/)
   10 format(1x,2i5,d20.5)
      end
c
c
      subroutine cntrac(cint,norb,n,nt,nro,buff)
c
c     transformation of integrals in mo-basis. only diagonalelements
c     are calculated
c
      implicit real*8(a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
      common /contrc/ chomo,ncont
      common/ioind/isl1, isl2
      common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
      dimension cint(nomx,1),buff(lbuff)
      rewind 2
      do 10 m=1,nt
      do 10 i=1,norb
   10 cint(i,m)=0.0d0
      nend=ncont*norb
      ix=lbuff
      do 50 k=1,ncont
      do 50 l=1,k
      do 40 m=1,nt
      ix=ix+1
      if(ix.le.lbuff)go to 20
      ix=1
      read(2)buff
   20 do 30 i=1,norb
      kk=ncont*(i-1)+k
      ll=ncont*(i-1)+l
      tmp=y(kk)*y(ll)
      if(k.ne.l) tmp=tmp+tmp
   30 cint(i,m)=buff(ix)*tmp+cint(i,m)
   40 continue
   50 continue
      return
      end
c
c
      subroutine cntrcm(cint,norb,non,buff,frocc,rnabcd)
c
c     mulliken and roby population analysis
c
      implicit real*8(a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
      integer centr
      common/scfinf/dmh(nomx),roc(nomx)
      common /big/ y(nomx*nomx),su(nomx*nomx),uv(nomx*nomx)
      common/holl/blabel(19),alabel(19),fmt(4),jl,lomao,iopen,nsym
      common/muli/vlist(nocmx,4),ianr(nomx),centr(nocmx),itnr(nomx)
      common /contrc/ chomo,ncont
      common/ioind/isl1, isl2
      common/supp/nsupp
      data ifile/4/,jfile/7/,kfile/9/,tol/0.d0/
      dimension cint(2),buff(lbuff),rnmy(20),frocc(2),mor(20)
c        number of aos                : 250
c        number of atoms              :  16
c        dimension of ovpop           : 16*(16+1)/2
c        dimension of nat and jatnr   : 16+1
      dimension ovpop(136),mos(250),iatnr(16),nat(17),rmopo(250,250)
      dimension rnab(15,16),jatnr(17),rnabc(14,15,16),rnamo(250,16)
      dimension rnabcd(1),katnr(16)
c        rnabcd(1) equivalenced to buff(16+1) in main program
      equivalence (uv(1),rmopo(1,1)), (rnabc(1,1,1),rnamo(1,1))
      tri=1.d0/3.d0
      four=1.d0/4.d0
      ipr=ncont*(ncont+1)/2
      ncmax=ncont*ncont
      maxuv=nomx*nomx
      idmax=maxuv/2
      nocm1=nocmx-1
      nocm2=nocmx-2
      nocp1=nocmx+1
c
c         mulliken population analysis
c
c                overlap * mos, written on kfile (for roby analysis)
      call hm(buff,y,su,norb,ncont)
      momax=norb*ncont
      rewind kfile
      write(kfile) (su(ir),ir=1,momax)
      rewind kfile
c
c               calculation of density , written on jfile
      do 205 i=1,ipr
  205 su(i)=0.d0
      do 210 k=1,norb
      ij=0
      ist=(k-1)*ncont
      do 220 i=1,ncont
      iy=ist+i
      do 220 j=1,i
      ij=ij+1
      jy=ist+j
  220 su(ij)=su(ij)+y(iy)*y(jy)*frocc(k)
  210 continue
      rewind jfile
      write(jfile) (su(ir),ir=1,ipr)
      rewind jfile
c
c
      do 10 i=1,ncont
   10 cint(i)=0.0d0
      ix=0
      iend=non*(non+1)/2
      ij=0
      do 15 i=1,iend
      ovpop(i)=0.d0
      do 15 j=1,norb
   15 rmopo(i,j)=0.d0
c                         calculation of overlap charges
      do 50 k=1,ncont
      do 50 l=1,k
      ij=ij+1
      ix=ix+1
      tmp=0.d0
      do 30 i=1,norb
      kk=ncont*(i-1)+k
      ll=ncont*(i-1)+l
      ka=ianr(k)
      la=ianr(l)
      kla=(ka*(ka-1))/2+la
      tmpb=y(kk)*y(ll)*frocc(i)*buff(ix)
      if(ka.ne.la) tmpb=tmpb*2.0d0
      rmopo(kla,i)=rmopo(kla,i)+tmpb
      if(k.eq.l) goto 30
      if(ka.eq.la) rmopo(kla,i)=rmopo(kla,i)+tmpb
   30 continue
      tmpb=su(ij)*buff(ix)
      cint(k)=cint(k)+tmpb
      if(k.ne.l) cint(l)=cint(l)+tmpb
   50 continue
      ren=0.d0
      do 150 ii=1,ncont
  150 ren=ren+cint(ii)
      write(6,1000)
c                             calculation of occupation numbers
      do 70 j=1,non
      do 71 k=1,20
   71 rnmy(k)=0.d0
      rnmy2=0.d0
      rnmy3=0.d0
      rnmy4=0.d0
      rnm=0.d0
      do 60 i=1,ncont
      if(ianr(i).ne.j) goto 60
      itp=itnr(i)
      rnmy(itp)=rnmy(itp)+cint(i)
   60 continue
      rnmy2=rnmy2+rnmy(2)+rnmy(3)+rnmy(4)
      do 80 k=5,10
   80 rnmy3=rnmy3+rnmy(k)
      do 90 k=11,20
   90 rnmy4=rnmy4+rnmy(k)
      rnm=rnmy(1)+rnmy2+rnmy3+rnmy4
      chg=vlist(j,4)-rnm
      write(6,2000)
      write(6,3000) j,centr(j),rnmy(1),rnmy2,rnmy3,rnmy4,rnm,chg
      if(rnmy2.ne.0.d0) write(6,4000) (rnmy(ii),ii=2,4)
      if(rnmy3.ne.0.d0) write(6,4001) (rnmy(ii),ii=5,10)
      if(rnmy4.ne.0.d0) write(6,4002) (rnmy(ii),ii=11,20)
   70 continue
      write(6,5000) ren
      write(6,6000)
      fakt=1.0d0
      do 115 i=1,iend
      do 115 j=1,norb
  115 ovpop(i)=ovpop(i)+rmopo(i,j)
      call outpak(ovpop,non,1,1)
      if(nsupp.ne.5) goto 218
      ii=0
      write(6,6400)
      do 110 i=1,non
      do 110 j=1,i
      write(6,6500) i,centr(i),j,centr(j)
      ii=ii+1
      write(6,1001) (rmopo(ii,k),k=1,norb)
  110 continue
c
c        roby population analysis
c
  218 continue
      write(6,8000)
      nnc=ncont
c
c                 calculate blocks of atoms
      itmp=1
      ij=0
      do 400 i=2,ncont
      if(ianr(i).ne.ianr(i-1)) goto 405
      itmp=itmp+1
      goto 400
  405 ij=ij+1
      iatnr(ij)=itmp
      itmp=1
  400 continue
      iatnr(ij+1)=itmp
c
      ndim=0
      do 733 i=1,non
      ia=iatnr(i)
  733 ndim=ndim+(ia*ia+ia)/2
c
c                 atom blocked overlap matrix on su
      call sort(buff,su,iatnr,non,0,ncmax)
c
c               calculation of s**1/2 and s**-1/2 (diagonal)
      ist=1
      jst=1
      ij=0
      ll=0
      do 530 i=1,non
      ia=iatnr(i)
      do 540 j=1,ia
      do 540 k=1,ia
      ij=ij+1
      uv(ij)=0.d0
      if(j.eq.k) uv(ij)=1.d0
  540 continue
      call erduw(su(jst),uv(ist),ia,1.d-12)
      ii=jst-1
      do 560 l=1,ia
      ii=ii+l
      ll=ll+1
      cint(ll)=sqrt(su(ii))
  560 dmh(ll)=1.d0/cint(ll)
      ist=ist+ia*ia
      jst=jst+(ia*ia+ia)/2
  530 continue
      itest=ist-1
      if(itest.le.idmax) goto 94
      write(6,7100) itest,idmax
 7100 format(1x,'itest = ',i8,'  greater than idmax = ',i8,
     1 ' in subroutine cntrcm',/,' enlarge nomx in main program')
      stop
   94 continue
c
c                       calculation of s**1/2 (atom blocked)
      call uduk(uv,cint,su,iatnr,non,ndim)
      rewind jfile
      read(jfile) (y(ii),ii=1,ipr)
      rewind jfile
c
c                       atom blocked density matrix on uv
      call sort(y,uv(idmax+1),iatnr,non,0,ncmax)
c                        s**1/2  *  density  *  s**1/2
      call h1h2h1(su,uv(idmax+1),cint,y,iatnr,non,ipr)
c
c                     calculation of s**-1/2 (atom blocked)
      call uduk(uv,dmh,su,iatnr,non,ndim)
c
c             calculation of maos   (written on uv)
      ij=0
      ist=1
      jst=1
      do 570 i=1,non
      ia=iatnr(i)
      do 580 j=1,ia
      iuv1=jst-1+(j-1)*ia
      iuv2=jst-1+j
      do 580 k=1,j
      ij=ij+1
      iuv=iuv1+k
      uv(iuv)=su(ij)
      iuv=iuv2+(k-1)*ia
  580 uv(iuv)=su(ij)
      call erduw(y(ist),uv(jst),ia,1.d-12)
      ist=ist+(ia*ia+ia)/2
      jst=jst+ia*ia
  570 continue
c           full mao matrix on y
      call sort(uv,y,iatnr,non,1,ncmax)
c           calculation of  s * mao
      call hm(buff,y,uv,ncont,ncont)
c           read in density and calculate  (mao*s)# * d * (mao*s)
      rewind jfile
      read(jfile) (buff(ii),ii=1,ipr)
      rewind jfile
      call ukhu(uv,buff,cint,su,ncont,ncont)
      rewind ifile
      write(ifile) (su(ir),ir=1,ncmax)
      rewind ifile
c
      ii=ncont+1
      ij=1
      do 585 i=1,ncont
      roc(i)=su(ij)
  585 ij=ij+ii
c                  mao * (s*mao), written on jfile
      call mkm(y,uv,su,ncont,ncont,ncont)
      rewind jfile
      write(jfile) (su(ir),ir=1,ncmax)
      rewind jfile
c
      do 649 i=1,ncont
  649 mos(i)=0
c
c            roby occupation numbers
      nnc=0
      ii=0
      ist=0
      do 600 i=1,non
      jatnr(i)=iatnr(i)
      ia=iatnr(i)
      rat=0.d0
      do 610 j=1,ia
      ii=ii+1
  610 rat=rat+roc(ii)
      write(6,6700) i,centr(i),rat
      write(6,6800) (roc(ir),ir=ist+1,ii)
c
c          selection of maos by user
      rat=0.d0
      read(5,1500) nmao,(mor(im),im=1,nmao)
      if(mor(1).eq.0) goto 851
      do 854 im=1,nmao
      kk=mor(im)
      mos(kk)=kk
      cint(im)=roc(kk)
  854 rat=rat+roc(kk)
      iatnr(i)=nmao
      nnc=nnc+nmao
      goto 855
  851 iatnr(i)=nmao
c          automatic selection of maos according to occupation numbers
      if(nmao.eq.0) write(6,6900) nmao,i
      if(nmao.eq.0) goto 600
      nnc=nnc+nmao
      iij=0
      do 850 k=1,nmao
      rpr=0.d0
      do 860 l=1,ia
      if(roc(ist+l).le.rpr) goto 860
      kk=ist+l
      rpr=roc(kk)
  860 continue
      iij=iij+1
      mos(kk)=kk
      rat=rat+rpr
      cint(iij)=rpr
      roc(kk)=-rpr
  850 continue
  855 rnmy(i)=rat
      write(6,6900) nmao,i,rat
      write(6,7500) (cint(ir),ir=1,nmao)
  600 ist=ist+ia
c
      do 871 i=1,non
  871 ianr(i)=jatnr(i)-iatnr(i)
      min=1000
      do 872 i=1,non
      if(ianr(i).lt.min) min=ianr(i)
  872 continue
c               hypervalence contributions
      call hypval(jatnr,mos,tol,ipr,ifile,non,ncmax,ren,
     1 jfile,rtot)
c
c
      nnc2=nnc*nnc
      mm=0
      do 870 m=1,ncont
      if(mos(m).eq.0) goto 870
      mm=mm+1
      mos(mm)=mos(m)
      roc(mm)=abs(roc(m))
  870 continue
c
c            calculate trace (mao*b*mao * s**-1)
c            calculation of unassigned charge
      rewind ifile
      read(ifile) (uv(iw),iw=1,ncmax)
      rewind ifile
      rewind jfile
      read(jfile) (su(iw),iw=1,ncmax)
      rewind jfile
c         order according to selected maos
      call order(uv,su,mos,nnc,ncont)
c         ordered overlap(mao basis) written on jfile
      rewind jfile
      write(jfile) (su(iw),iw=1,nnc2)
      rewind jfile
c
      if(nsupp.ne.4) goto 655
      write(6,9100)
      call prtblk(su,nnc,nnc,'mao ','mao ',mos)
      write(6,9200)
  655 call osinv1(su,nnc,d,tol,itnr,ianr)
c
      rna=0.d0
      do 910 i=1,nnc
      a=0.d0
      i1=i
      i20=(i-1)*nnc
      do 920 j=1,nnc
      a=a+su(i1)*uv(i20+j)
  920 i1=i1+nnc
      rna=rna+a
  910 continue
      uach=ren-rna
      write(6,7700) uach
      rtot=(uach-rtot)/non
      do 500 i=1,non
  500 dmh(i)=dmh(i)+rtot
      write(6,7800) rtot
c
c               selection and reordering of mao vector y
      imo=0
      ie=0
      ij=0
      nnc=0
      nst=1
      nat(1)=1
      do 650 l=1,non
      nmao=iatnr(l)
      nat(l+1)=nat(l)+nmao
      if(nmao.eq.0) goto 650
      do 660 i=1,nmao
      imo=imo+1
      imos=mos(imo)-1
      ist=imos*ncont+1
      ie=ist+ncont-1
      do 680 k=ist,ie
      ij=ij+1
  680 y(ij)=y(k)
  660 continue
      write(6,7600) l,centr(l)
      call prtblk(y(nst),ncont,nmao,' ao ','mao ',mos(nnc+1))
      nst=nst+nmao*ncont
      nnc=nnc+nmao
  650 continue
c                       write out maos for plot
      ndimw=nnc*ncont
      write(lomao,8900) (alabel(i1),i1=1,19)
      write(lomao,8900) (blabel(i1),i1=1,19)
      write(lomao,9000) nnc,ncont,ndimw
      write(lomao,8900) (fmt(i1),i1=1,4)
      write(lomao,fmt) (y(i),i=1,ndimw)
c
c            calculation of occupation numbers (mo-contributions)
      rewind kfile
      read(kfile) (uv(ir),ir=1,momax)
      rewind kfile
c                      (maos)# * (sc) = su
      call mkm(y,uv,su,nnc,ncont,norb)
      if(nsupp.eq.5) write(6,8650)
      do 960 m=1,non
      rat=0.d0
      mm=0
      do 970 j=1,norb
      ind=(j-1)*nnc
      do 980 k=nat(m),nat(m+1)-1
      mm=mm+1
  980 uv(mm)=su(ind+k)
  970 continue
      ndim=iatnr(m)
      i0=-ndim
      do 950 i=1,norb
      v=0.d0
      i0=i0+ndim
      do 955 j=1,ndim
  955 v=v+uv(i0+j)*uv(i0+j)
      rnamo(i,m)=v*frocc(i)
      rat=rat+rnamo(i,m)
  950 continue
      if(nsupp.eq.5) write(6,8655) m,centr(m),(rnamo(ir,m),ir=1,norb)
  960 continue
c
c            two center shared electron numbers
      write(6,8350)
      ncen=1
      ishare=2
      rewind jfile
      read(jfile) (uv(ir),ir=1,nnc2)
      rewind jfile
      do 715 i=1,non
  715 buff(i)=0.d0
      nat(non+1)=nnc+1
c
      do 700 i=1,non-1
      ij=0
      if(iatnr(i).eq.0) goto 700
      do 710 ii=nat(i),nat(i+1)-1
      ij=ij+1
  710 mos(ij)=ii
      do 720 j=i+1,non
      if(iatnr(j).eq.0) goto 720
      ik=ij
      do 730 jj=nat(j),nat(j+1)-1
      ik=ik+1
  730 mos(ik)=jj
      itest=ik*ik
      if(itest.gt.idmax) goto 95
c
      call sen(su,uv,y(1),y(idmax),mos,cint,rna,ik,d,tol,itnr,ianr,
     1 norb,nnc,roc,frocc)
      rnab(i,j)=rna
      write(6,7900) i,centr(i),j,centr(j)
      rna=0.d0
      do 810 im=1,norb
      roc(im)=rnamo(im,i)+rnamo(im,j)-roc(im)
  810 rna=rna+roc(im)
      write(6,8100) (roc(ir),ir=1,norb)
      buff(i)=buff(i)+rna*0.5d0
      buff(j)=buff(j)+rna*0.5d0
      write(6,8200) rna
  720 continue
  700 continue
c
c                three center shared electron numbers
c
      if(non.le.3) goto 515
      ncen=-1
      ishare=3
      write(6,8400)
      do 740 i=1,non-2
      ij=0
      if(iatnr(i).eq.0) goto 740
      do 745 ii=nat(i),nat(i+1)-1
      ij=ij+1
  745 mos(ij)=ii
      do 750 j=i+1,non-1
      if(iatnr(j).eq.0) goto 750
      ik=ij
      do 755 jj=nat(j),nat(j+1)-1
      ik=ik+1
  755 mos(ik)=jj
      do 760 k=j+1,non
      if(iatnr(k).eq.0) goto 760
      il=ik
      do 765 kk=nat(k),nat(k+1)-1
      il=il+1
  765 mos(il)=kk
      itest=il*il
      if(itest.gt.idmax) goto 95
c
      call sen(su,uv,y(1),y(idmax),mos,cint,rna,il,d,tol,itnr,ianr,
     1 norb,nnc,roc,frocc)
      rnabc(i,j,k)=rna
      rna=rnmy(i)+rnmy(j)+rnmy(k)-rnab(i,j)-rnab(i,k)-rnab(j,k)+rna
      buff(i)=buff(i)-tri*rna
      buff(j)=buff(j)-tri*rna
      buff(k)=buff(k)-tri*rna
      if(abs(rna).ge.0.01d0) write(6,8500) i,j,k,rna
  760 continue
  750 continue
  740 continue
c
c                four center shared electron numbers
c
      if(non.le.4) goto 515
      ncen=1
      ishare=4
      write(6,8450)
      do 770 i=1,non-3
      ij=0
      if(iatnr(i).eq.0) goto 770
      do 775 ii=nat(i),nat(i+1)-1
      ij=ij+1
  775 mos(ij)=ii
      do 780 j=i+1,non-2
      if(iatnr(j).eq.0) goto 780
      ik=ij
      do 785 jj=nat(j),nat(j+1)-1
      ik=ik+1
  785 mos(ik)=jj
      do 790 k=j+1,non-1
      if(iatnr(k).eq.0) goto 790
      il=ik
      do 795 kk=nat(k),nat(k+1)-1
      il=il+1
  795 mos(il)=kk
      do 800 l=k+1,non
      if(iatnr(l).eq.0) goto 800
      im=il
      do 805 ll=nat(l),nat(l+1)-1
      im=im+1
  805 mos(im)=ll
      itest=im*im
      if(itest.gt.idmax) goto 95
c
      call sen(su,uv,y(1),y(idmax),mos,cint,rna,im,d,tol,itnr,ianr,
     1 norb,nnc,roc,frocc)
      index1=j*(j-3)/2+i
      index2=l*(l-3)/2+k
      index=index2*(index2-3)/2+index1
      rnabcd(index)=rna
      rna=rnmy(i)+rnmy(j)+rnmy(k)+rnmy(l)-rnab(i,j)-rnab(i,k)-rnab(i,l)
     1    -rnab(j,k)-rnab(j,l)-rnab(k,l)+rnabc(i,j,k)+rnabc(i,j,l)
     2    +rnabc(i,k,l)+rnabc(j,k,l)-rna
      buff(i)=buff(i)+four*rna
      buff(j)=buff(j)+four*rna
      buff(k)=buff(k)+four*rna
      buff(l)=buff(l)+four*rna
      if(abs(rna).ge.0.002d0) write(6,8550) i,j,k,l,rna
  800 continue
  790 continue
  780 continue
  770 continue

  515 read(5,6100) nm,(iatnr(i),i=1,nm)
      if(nm.eq.0) goto 933
      call spemu(su,uv,y(1),y(idmax),mos,cint,rna,d,tol,itnr,ianr,
     1 norb,nnc,roc,frocc,iatnr,rnmy,rnab,rnabc,rnabcd,nm,nat,idmax,
     2 jatnr,katnr,nocm1,nocm2,nocp1)
c
  933 continue
c
c         calculation of atomic charges
c           chomo: charge of molecule, cmult: multicenter contribution
      sum=0.d0
      do 786 i=1,non
      rnmy(i)=vlist(i,4)-rnmy(i)+buff(i)-dmh(i)
      sum=sum+rnmy(i)
  786 continue
      sum=sum-chomo
      cmult=sum*ncen
      if(non.eq.3) write(6,8840) cmult
      if(non.eq.4) write(6,8850) cmult
      if(non.ge.5) write(6,8800) cmult
      write(6,8600)
      sum=sum/non
      do 787 i=1,non
      q=rnmy(i)-sum
      write(6,8700) i,centr(i),q
  787 continue

c
c
 1000 format(///,' mulliken population analysis',//)
 1001 format(1x,12f10.6)
 1500 format(i10,20i3)
 2000 format(//,' atom nr. symbol  s-occupation  p-occupation  ',
     1 'd-occupation  f-occupation      total        charge  ')
 3000 format(1x,i5,4x,a4,2x,6(f12.5,2x))
 4000 format(/,1x,'n(px)  =',f9.5,2x,'n(py)  =',f9.5,2x,'n(pz)  =',f9.5)
 4001 format(1x,'n(dxx) =',f9.5,2x,'n(dyy) =',f9.5,2x,'n(dzz) =',f9.5,
     1 2x,'n(dxy) =',f9.5,2x,'n(dxz) =',f9.5,2x,'n(dyz) =',f9.5)
 4002 format(1x,'n(fxxx)=',f9.5,2x,'n(fyyy)=',f9.5,2x,'n(fzzz)=',f9.5,
     1 2x,'n(fxxy)=',f9.5,2x,'n(fxxz)=',f9.5,2x,'n(fyyx)=',f9.5,/,
     2 1x,'n(fyyz)=',f9.5,2x,'n(fzzx)=',f9.5,2x,'n(fzzy)=',f9.5,
     3 2x,'n(fxyz)=',f9.5)
 5000 format(//,' total number of electrons = ',f15.8)
 6000 format(///,' mulliken overlap population ',/)
 6100 format(20i3)
 6400 format(//,'      overlap population per mo',/)
 6500 format(/,' atom ',i5,a4,'    /    ','atom ',i5,a4)
 6700 format(//,' atom nr ',i5,2x,a4,'  occupation',f10.6)
 6800 format(1x,'occupation of maos',/,(10f10.6))
 6900 format(1x,i5,'  maos selected for atom ',i4,'  occupation ',f10.6)
 7500 format(1x,'occupation of maos ',/,(10f10.6))
 7600 format(///,' selected maos for atom ',i3,a4)
 7700 format(/,10x,'unassigned charge = ',f10.4,/)
 7800 format(1x,'rest charge per atom = ',f10.6)
 7900 format(//,' atom ',i3,a4,'  /  atom ',i3,a4)
 8000 format(///,' roby population analysis with modified atomic',
     1 ' orbitals   ')
 8100 format(1x,'contributions of mos to shared electron number',
     1 /,(10f10.5))
 8200 format(1x,' shared electron number = ',f10.4)
 8300 format(1x,'n( 1 2 3 )  = ',f10.5)
 8350 format(///,' two center shared electron numbers '/)
 8400 format(///,' three center shared electron numbers ge 0.01'/)
 8450 format(///,' four center shared electron numbers ge 0.002'/)
 8500 format(1x,'n(',3i3,') = ',f10.4,f10.6)
 8550 format(1x,'n(',4i3,') = ',f10.4)
 8600 format(//,' atomic charges with multicenter corrections',
     1 /,'    atom     charge')
 8650 format(//,' contributions of mos to occupation numbers')
 8655 format(/,' atom nr ',i5,2x,a4,/,(12f10.6))
 8700 format(1x,i5,a4,f10.4)
 8800 format(//,' multicenter rest contribution = ',f10.6,
     1 '  (correct sign for five center terms)')
 8840 format(//,' three center shared electron number = ',f10.6)
 8850 format(//,' four center shared electron number = ',f10.6)
 8900 format(19a4)
 9000 format(20i4)
 9100 format(///' s matrix in mao basis'/)
 9200 format(//)
 9500 format(a1,i4,14i5)
 9550 format(1x,2i5)
      return
c
   95 write(6,7200) itest,idmax,ishare
 7200 format(1x,'itest = ',i8,'  greater than idmax = ',i8,
     1 ' in subroutine cntrcm for',i3,'-center sens',/,
     2 ' enlarge all arrays of dimension nomx')
      stop
      end
c
c
      subroutine spemu(su,uv,ssu,suv,mos,cint,rna,d,tol,itnr,ianr,
     1 norb,nnc,roc,frocc,iatnr,rnmy,rnab,rnabc,rnabcd,nm,nat,idmax,
     2 jatnr,katnr,nocm1,nocm2,nocp1)
      implicit real*8(a-h,o-z)
      dimension su(1),uv(1),ssu(1),suv(1),mos(1),cint(1),itnr(1),ianr(1)
      dimension roc(1),frocc(1),iatnr(1),rnmy(1),rnab(nocm1,1),nat(1)
      dimension rnabc(nocm2,nocm1,1),rnabcd(1),jatnr(nocp1),katnr(1)
c
      iend=nm
      do 310 j=1,nm-1
      max=iatnr(1)
      ind=1
      do 320 i=2,iend
      if(iatnr(i).lt.max) goto 320
      max=iatnr(i)
      ind=i
  320 continue
      jj=nm-j+1
      iend=jj-1
      iatnr(ind)=iatnr(jj)
  310 iatnr(jj)=max
      write(6,2000) (iatnr(i),i=1,nm)
c
      ij=0
      do 50 i=1,nm
      i1=iatnr(i)
      do 60  ii=nat(i1),nat(i1+1)-1
      ij=ij+1
   60 mos(ij)=ii
   50 continue
      itest=ij*ij
      if(itest.gt.idmax) goto 99
      call sen(su,uv,ssu,suv,mos,cint,rna,ij,d,tol,itnr,ianr,
     1 norb,nnc,roc,frocc)
      tmult=0.d0
      do 100 i=1,nm
      i1=iatnr(i)
      tmult=tmult+rnmy(i1)
  100 continue
      tmult=tmult-rna
      write(6,3500) tmult
      if(nm.le.4) return
      ipha=(-1)**(nm-1)
      sum=ipha*rna
c
      lmax=int(4.e+09)
      do 10 i=1,nm
      i1=iatnr(i)
      sum=sum+rnmy(i1)
      if(i.eq.1) goto 10
      do 20 j=1,i-1
      j1=iatnr(j)
      sum=sum-rnab(j1,i1)
      if(j.eq.1) goto 20
      do 30 k=1,j-1
      k1=iatnr(k)
      sum=sum+rnabc(k1,j1,i1)
      if(k.eq.1) goto 30
      do 40 l=1,k-1
      l1=iatnr(l)
      index1=k1*(k1-3)/2+l1
      index2=i1*(i1-3)/2+j1
      index=index2*(index2-3)/2+index1
   40 sum=sum-rnabcd(index)
   30 continue
   20 continue
   10 continue
      if(nm.eq.5) goto 410
c
      do 110 i=5,nm-1
      do 115 ii=1,i
      katnr(ii)=ii
  115 jatnr(ii)=iatnr(ii)
      ij=0
      do 81 k=1,i
      i1=jatnr(k)
      do 91  ii=nat(i1),nat(i1+1)-1
      ij=ij+1
   91 mos(ij)=ii
   81 continue
      itest=ij*ij
      if(itest.gt.idmax) goto 99
      call sen(su,uv,ssu,suv,mos,cint,rna,ij,d,tol,itnr,ianr,
     1 norb,nnc,roc,frocc)
      ipha=(-1)**(i-1)
      sum=sum+ipha*rna
c
      jatnr(i+1)=1000
      jst=i
      do 140 l=1,lmax
      do 120 j=jst,i
      do 130 k=jst+1,nm
      if(iatnr(k).le.jatnr(j)) goto 130
      if(iatnr(k).ge.jatnr(j+1)) goto 130
      jatnr(j)=iatnr(k)
      katnr(j)=k
      ij=0
      do 80 ll=1,i
      i1=jatnr(ll)
      do 90  ii=nat(i1),nat(i1+1)-1
      ij=ij+1
   90 mos(ij)=ii
   80 continue
      itest=ij*ij
      if(itest.gt.idmax) goto 99
      call sen(su,uv,ssu,suv,mos,cint,rna,ij,d,tol,itnr,ianr,
     1 norb,nnc,roc,frocc)
      ipha=(-1)**(i-1)
      sum=sum+ipha*rna
  130 continue
  120 continue
      ki=0
      do 145 kk=1,i
      if(jatnr(kk).ne.iatnr(nm-i+kk)) goto 145
      ki=ki+1
      if(ki.eq.i) goto 110
      if(ki.eq.1) jst=kk-1
      kkk=katnr(kk-1)+1
      if(ki.eq.1) kkk=katnr(kk-1)+2
      jatnr(kk)=iatnr(kkk)
      katnr(kk)=kkk
  145 continue
  140 continue
  110 continue
c
  410 write(6,3000) sum
      return
c
 1000 format(a1,i2,20i3)
 2000 format(///,' special multi center contribution for the atoms',/,
     1 (20i5))
 3000 format(1x,'multi center shared electron number  =',f15.6,/)
 3500 format(1x,'total contribution =',f15.6)
c
   99 write(6,4000) itest,idmax
 4000 format(1x,'itest = ',i8,'  greater than idmax = ',i8,
     1 ' in subroutine spemu for special multi center contributions',/,
     2 ' enlarge all arrays of dimension nomx')
      return
      end
c
c
      subroutine sen(su,uv,ssu,suv,mos,cint,rna,idim,d,tol,itnr,
     1 ianr,norb,nnc,roc,frocc)
c
c        calculation of shared electron numbers
c
      implicit real*8 (a-h,o-z)
      dimension su(1),uv(1),ssu(1),suv(1),mos(1),cint(1),itnr(1)
      dimension ianr(1),roc(1),frocc(1)
c
c             sort (mao)#*s*c on ssu
      ij=0
      do 10 i=1,norb
      ind0=(i-1)*nnc
      do 20 j=1,idim
      ij=ij+1
      ind=ind0+mos(j)
   20 ssu(ij)=su(ind)
   10 continue
c
c             sort s(mao basis) on suv
      ij=0
      do 30 j=1,idim
      ind0=(mos(j)-1)*nnc
      do 40 i=1,idim
      ind=ind0+mos(i)
      ij=ij+1
   40 suv(ij)=uv(ind)
   30 continue
c
      call osinv1(suv,idim,d,tol,itnr,ianr)
      ij=0
      do 60 i=1,idim
      ind=(i-1)*idim
      do 60 j=1,i
      ij=ij+1
      suv(ij)=suv(ind+j)
   60 continue
      call ukhud(ssu,suv,cint,roc,idim,norb)
c
      rna=0.d0
      do 50 i=1,norb
      roc(i)=roc(i)*frocc(i)
   50 rna=rna+roc(i)
      return
      end
c
c
      subroutine ukhud(a,b,c,d,idim1,idim2)
c
c             u# * triangle matrix * u , diagonalelements only
      implicit real*8(a-h,o-z)
      dimension a(1),b(1),c(1),d(1)
      kl=0
      ij=0
      i10=-idim1
      do 40 j=1,idim2
      i10=i10+idim1
      ij=0
      do 50 k=1,idim1
      cc=0.d0
      ij=ij+1
      i20=(k*k-k)/2
      do 60 l=1,k
   60 cc=cc+a(i10+l)*b(i20+l)
      if(k.eq.idim1) goto 55
      i20=k*(k+3)/2
      do 70 l=k+1,idim1
      cc=cc+a(i10+l)*b(i20)
   70 i20=i20+l
   55 c(ij)=cc
   50 continue
      m1=0
      m2=(j-1)*idim1
      kl=kl+1
      dd=0.0d0
      do 90 n=1,idim1
   90 dd=dd+c(m1+n)*a(m2+n)
      d(kl)=dd
   40 continue
      return
      end
c
      subroutine hypval(jatnr,mos,tol,ipr,ifile,non,
     1 ncmax,ren,jfile,rtot)
c
c         calculation of hypervalence populations
c
      implicit real*8 (a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
      integer centr
      dimension jatnr(1),mos(1)
      common/scfinf/dmh(nomx),roc(nomx)
      common /big/ y(nomx*nomx),su(nomx*nomx),uv(nomx*nomx)
      common/muli/vlist(nocmx,4),ianr(nomx),centr(nocmx),itnr(nomx)
      common /contrc/ chomo,ncont
      iend=0
      rtot=0.d0
      write(6,1000)
      do 100 ia=1,non
      ist=iend+1
      iend=iend+jatnr(ia)
      nnc=ist-1
      if(nnc.eq.0) goto 111
      do 50 i=1,nnc
   50 itnr(i)=i
  111 continue
      do 110 ii=ist,iend
      if(mos(ii).eq.0) goto 110
      nnc=nnc+1
      itnr(nnc)=ii
  110 continue
      if(ia.eq.non) goto 121
      do 120 ii=iend+1,ncont
      nnc=nnc+1
      itnr(nnc)=ii
  120 continue
c
  121 rewind ifile
      read(ifile) (uv(iw),iw=1,ncmax)
      rewind ifile
      rewind jfile
      read(jfile) (su(iw),iw=1,ncmax)
      rewind jfile
      call order(uv,su,itnr,nnc,ncont)
      call osinv1(su,nnc,d,tol,itnr,ianr)
c
c            calculate trace(mao*s*d*s*mao * s**-1)
      rna=0.d0
      do 910 i=1,nnc
      a=0.d0
      i1=i
      i20=(i-1)*nnc
      do 920 j=1,nnc
      a=a+su(i1)*uv(i20+j)
  920 i1=i1+nnc
      rna=rna+a
  910 continue
      rna=ren-rna
      dmh(ia)=rna
      rtot=rtot+rna
      write(6,2000) ia,centr(ia),rna
  100 continue
c
 1000 format(/'       hypervalence contributions'//' unshared atomic',
     1 ' populations')
 2000 format(1x,'atom nr.',i3,1x,a4,5x,f10.6)
      return
      end
c
      subroutine order(a,b,iord,idims,idimb)
      implicit real*8(a-h,o-z)
      dimension a(1),b(1),iord(1)
c           sort matrix a an b according to vector iord
      ij=0
      do 10 j=1,idims
      ind2=(iord(j)-1)*idimb
      do 20 i=1,idims
      ind=ind2+iord(i)
      ij=ij+1
      a(ij)=a(ind)
   20 b(ij)=b(ind)
   10 continue
      return
      end
c
c
      subroutine hm(a,b,c,norb,ncont)
c                   triangle matrix(a) * full matrix(b)
      implicit real*8(a-h,o-z)
      dimension a(1),b(1),c(1)
      ij=0
      do 20 i=1,norb
      ibst=(i-1)*ncont
      do 30 j=1,ncont
      ct=0.d0
      ij=ij+1
      iast=j*(j-1)/2
      do 40 k=1,j
   40 ct=ct+a(iast+k)*b(ibst+k)
      if(j.eq.ncont) goto 29
      iast=j*(j+3)/2
      do 50 k=j+1,ncont
      ct=ct+a(iast)*b(ibst+k)
   50 iast=iast+k
   29 c(ij)=ct
   30 continue
   20 continue
      return
      end
c
      subroutine mkm(a,b,c,idima,idimab,idimb)
      implicit real*8(a-h,o-z)
      dimension a(1),b(1),c(1)
c         (full matrix)# * full matrix
      ij=0
      do 10 i=1,idimb
      iast=0
      ibst=(i-1)*idimab
      do 20 j=1,idima
      ij=ij+1
      v=0.d0
      do 30 k=1,idimab
   30 v=v+a(iast+k)*b(ibst+k)
      c(ij)=v
   20 iast=iast+idimab
   10 continue
      return
      end
c
c
c
c
      subroutine sort(a,b,nr,non,iswi,ncmax)
c         sorting of matrices in atom blocked matrices
      implicit real*8(a-h,o-z)
      dimension a(1),b(1),nr(1)
c                iswi=0: full matrix in atom blocked matrix
      ij=0
      if(iswi.eq.1) goto 30
      jend=0
      do 10 i=1,non
      jst=jend+1
      jend=jend+nr(i)
      do 10 j=jst,jend
      kkst=(j*j-j)/2
      do 10 k=jst,j
      ij=ij+1
      kl=kkst+k
      b(ij)=a(kl)
   10 continue
      return
c                iswi=1: atom blocked matrix in full matrix
   30 iofs=0
      ndim=0
      do 35 i=1,non
   35 ndim=ndim+nr(i)
      do 38 i=1,ncmax
   38 b(i)=0.d0
      iaa=0
      kl=0
      do 40 i=1,non
      ia=nr(i)
      iaa=iaa+ia
      do 50 j=1,ia
      do 60 k=1,ia
      ij=ij+1
      kl=kl+1
   60 b(kl)=a(ij)
      kl=kl+ndim-iaa+iofs
   50 continue
      kl=kl+ia
      iofs=iofs+ia
   40 continue
      return
      end
c
c
c
c
c
c
      subroutine uduk(a,b,c,nr,non,ndim)
c              u# * vector * u
      implicit real*8(a-h,o-z)
      dimension a(1),b(1),c(1),nr(1)
      do 10 i=1,ndim
   10 c(i)=0.d0
      ii=0
      ijst=0
      iist=0
      do 20 l=1,non
      ia=nr(l)
      do 30 i=1,ia
      ii=ii+1
      ij=ijst
      ist=iist+(i-1)*ia
      do 40 j=1,ia
      iu=ist+j
      do 40 k=1,j
      ij=ij+1
      ju=ist+k
   40 c(ij)=c(ij)+a(iu)*b(ii)*a(ju)
   30 continue
      ijst=ij
      iist=iu
   20 continue
      return
      end
c
c
      subroutine ukhu(a,b,c,d,ncont,nc)
c            u# * triangle matrix * u
      implicit real*8(a-h,o-z)
      dimension a(1),b(1),c(1),d(1)
      i10=-ncont
      do 40 j=1,nc
      i10=i10+ncont
      ij=0
      i11=i10+1
      do 50 k=1,ncont
      ij=ij+1
      i20=(k*k-k)/2
      i21=i20+1
      cij = 0.0d0
      eij = 0.0d0
      if (k. eq. 1) go to 110
      do 60 l = 1, k-1, 2
      cij = cij + a(i10+l) * b(i20+l)
   60 eij = eij + a(i11+l) * b(i21+l)
      if (mod(k,2) .eq. 0) go to 120
  110 continue
      cij = cij + a(i10+k) * b(i20+k)
  120 continue
      if(k.eq.ncont) goto 49
      loff = (k+3) * k / 2
      if (k+1 .eq .ncont) go to 130
      do 70 l=k+1,ncont-1,2
      cij=cij+a(i10+l) * b(loff)
      eij=eij+a(i11+l) * b(loff+l)
   70 loff = loff + l + l + 1
      if (mod(ncont-k,2) .eq. 0) go to 49
  130 cij=cij+a(i10+ncont) * b(loff)
   49 c(ij) = cij+eij
   50 continue
      m2=0
      do 80 m=1,nc
      kl=(m-1)*nc+j
      dkl=0.d0
      fkl=0.0d0
      m21 = m2 + 1
      if (ncont .eq. 1) go to 140
      do 90 n=1,ncont-1,2
      dkl=dkl+c(n)*a(m2+n)
   90 fkl=fkl+c(n+1)*a(m21+n)
      if(mod(ncont,2) .eq. 0) go to 150
  140 dkl = dkl + c(ncont) * a(m2+ncont)
  150 continue
      d(kl)=dkl+fkl
      m2=m2+ncont
   80 continue
   40 continue
      return
      end
c
c
      subroutine osinv1(a,n,d,tol,l,m)
c            inversion of a full matrix
      implicit real*8 (a-h,o-z)
      dimension a(1),m(1),l(1)
      if(n.ne.1) goto 5
      a(1)=1.0/a(1)
      return
    5 d=1.d0
      nk=-n
      do 180 k=1,n
      nk=nk+n
      l(k)=k
      m(k)=k
      kk=nk+k
      biga=a(kk)
      do 20 j=k,n
      iz=n*(j-1)
      do 20 i=k,n
      ij=iz+i
      if(abs(biga)-abs(a(ij))) 10,20,20
   10 biga=a(ij)
      l(k)=i
      m(k)=j
   20 continue
      j=l(k)
      if(j-k) 50,50,30
   30 ki=k-n
      do 40 i=1,n
      ki=ki+n
      holo=-a(ki)
      ji=ki-k+j
      a(ki)=a(ji)
   40 a(ji)=holo
   50 i=m(k)
      if(i-k) 80,80,60
   60 jp=n*(i-1)
      do 70 j=1,n
      jk=nk+j
      ji=jp+j
      holo=-a(jk)
      a(jk)=a(ji)
   70 a(ji)=holo
   80 if(abs(biga)-tol) 90,100,100
   90 d=0.d0
      return
  100 continue
      biginv = -1.0d0 / biga
      do 120 i=1,n
      a(nk+i)=a(nk+i)*biginv
  120 continue
      a(nk+k) = a(nk+k) * (-biga)
      joff = 0
      do 155 j=1,n
      if(j. eq. k) go to 154
      ajk = a(joff+k)
      if (k .eq. 1) go to 152
      do 151 i=1,k-1
  151 a(joff+i) = a(joff+i) + a(nk+i) * ajk
  152 continue
      if (k .eq. n) go to 154
      do 153 i=k+1,n
  153 a(joff+i) = a(joff+i) + a(nk+i) * ajk
  154 continue
      joff = joff + n
  155 continue
      kj=k-n
      biginv=-biginv
      do 170 j=1,n
      kj=kj+n
      a(kj)=a(kj)*biginv
  170 continue
      kadd = (k-1)*n + k
      a(kadd) = a(kadd) * biga
      d=d*biga
      a(kk)=biginv
  180 continue
      k=n
  190 k=k-1
      if(k) 260,260,200
  200 i=l(k)
      if(i-k) 230,230,210
  210 jq=n*(k-1)
      jr=n*(i-1)
      do 220 j=1,n
      jk=jq+j
      holo=a(jk)
      ji=jr+j
      a(jk)=-a(ji)
  220 a(ji)=holo
  230 j=m(k)
      if(j-k) 190,190,240
  240 ki=k-n
      do 250 i=1,n
      ki=ki+n
      holo=a(ki)
      ji=ki+j-k
      a(ki)=-a(ji)
  250 a(ji)=holo
      goto 190
  260 return
      end
c
c
      subroutine h1h2h1(a,b,c,d,nr,non,max)
c        matrix * matrix * matrix  (three triangle matrices)
      implicit real*8(a-h,o-z)
      dimension a(1),b(1),c(1),d(1),nr(1)
      do 10 i=1,max
   10 d(i)=0.d0
      iofs=0
      kl=0
      do 20 i=1,non
      ia=nr(i)
      do 30 j=1,ia
      ij=0
      i10=iofs+j*(j-1)/2
      jj=iofs+j
      do 40 k=1,ia
      i2=iofs+k*(k-1)/2
      ij=ij+1
      c(ij)=0.d0
      do 50 l=1,k
      i1=i10+l
      if(l.gt.j) i1=l*(l-1)/2+jj
      i2=i2+1
   50 c(ij)=c(ij)+a(i1)*b(i2)
      if(k.eq.ia) goto 40
      i20=iofs+k
      do 55 l=k+1,ia
      i1=i10+l
      if(l.gt.j) i1=l*(l-1)/2+jj
      i2=i20+l*(l-1)/2
   55 c(ij)=c(ij)+a(i1)*b(i2)
   40 continue
      do 60 m=1,j
      mm=iofs+m
      kl=kl+1
      m20=m*(m-1)/2+iofs
      do 70 n=1,ia
      m2=m20+n
      if(n.gt.m) m2=n*(n-1)/2+mm
   70 d(kl)=d(kl)+c(n)*a(m2)
   60 continue
   30 continue
      iofs=iofs+ia*(ia+1)/2
   20 continue
      return
      end
c
c
      function noteqv(i,j)
      implicit real*8(a-h,o-z)
      noteqv = i-j
      return
      end
c
c
      real*8 function ggntde(i,j,eta,nfmx,nfcol)
      implicit real*8(a-h,o-z)
      dimension eta(nfmx,nfcol)
      dimension aa(3),bb(3),a(10),b(10),e(3)
      data c2/2.0d0/,c3/3.0d0/,const/15.7496099d0/
      do 1 m=1,3
      aa(m) = eta(i,m)
    1 bb(m) = eta(j,m)
      iff1=eta(i,5)+0.0001
      iff2=eta(j,5)+0.0001
      do 2 m=1,10
      a(m) = eta(i,m+5)
    2 b(m) = eta(j,m+5)
      call divpt(aa,eta(i,4),bb,eta(j,4),e,eexp,efact)
      call rhftce(a,aa,e,iff1)
      call rhftce(b,bb,e,iff2)
      fexp=c2*eexp
      ggntde =                          efact*const*(a(1)*b(1)+(a(2)*b(2
     1)+a(3)*b(3)+a(4)*b(4)+a(1)*(b(5)+b(6)+b(7))+b(1)*(a(5)+a(6)+a(7))+
     2 (c3*(a(5)*b(5)+a(6)*b(6)+a(7)*b(7))+a(8)*b(8)+a(9)*b(9)+a(10)*b(1
     30)+a(5)*(b(6)+b(7))+a(6)*(b(5)+b(7))+a(7)*(b(5)+b(6)))/fexp)/fexp)
     4/(fexp*sqrt(fexp))
      return
      end
c
c
      subroutine rhftce(cfs,a,e,iff)
      implicit real*8(a-h,o-z)
      dimension cfs(20) ,a(3),e(3)
      data c2/2.0d0/,c3/3.0d0/
      aex = e(1)-a(1)
      aey = e(2)-a(2)
      aez = e(3)-a(3)
      goto (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20),iff
    1 return
    2 cfs(1)=aex*cfs(2)
      return
    3 cfs(1)=aey*cfs(3)
      return
    4 cfs(1)=aez*cfs(4)
      return
    5 cfs(1)=aex*aex*cfs(5)
      cfs(2)=c2*aex*cfs(5)
      return
    6 cfs(1)=aey*aey*cfs(6)
      cfs(3)=c2*aey*cfs(6)
      return
    7 cfs(1)=aez*aez*cfs(7)
      cfs(4)=c2*aez*cfs(7)
      return
    8 cfs(1)=aex*aey*cfs(8)
      cfs(2)=aey*cfs(8)
      cfs(3)=aex*cfs(8)
      return
    9 cfs(1)=aex*aez*cfs(9)
      cfs(2)=aez*cfs(9)
      cfs(4)=aex*cfs(9)
      return
   10 cfs(1)=aey*aez*cfs(10)
      cfs(3)=aez*cfs(10)
      cfs(4)=aey*cfs(10)
      return
   11 cfs(1)=aex*aex*aex*cfs(11)
      cfs(2)=c3*aex*aex*cfs(11)
      cfs(5)=c3*aex*cfs(11)
      return
   12 cfs(1)=aey*aey*aey*cfs(12)
      cfs(3)=c3*aey*aey*cfs(12)
      cfs(6)=c3*aey*cfs(12)
      return
   13 cfs(1)=aez*aez*aez*cfs(13)
      cfs(4)=c3*aez*aez*cfs(13)
      cfs(7)=c3*aez*cfs(13)
      return
   14 cfs(1)=aex*aex*aey*cfs(14)
      cfs(2)=c2*aex*aey*cfs(14)
      cfs(3)=aex*aex*cfs(14)
      cfs(5)=aey*cfs(14)
      cfs(8)=c2*aex*cfs(14)
      return
   15 cfs(1)=aex*aex*aez*cfs(15)
      cfs(2)=c2*aex*aez*cfs(15)
      cfs(4)=aex*aex*cfs(15)
      cfs(5)=aez*cfs(15)
      cfs(9)=c2*aex*cfs(15)
      return
   16 cfs(1)=aey*aey*aex*cfs(16)
      cfs(2)=aey*aey*cfs(16)
      cfs(3)=c2*aey*aex*cfs(16)
      cfs(6)=aex*cfs(16)
      cfs(8)=c2*aey*cfs(16)
      return
   17 cfs(1)=aey*aey*aez*cfs(17)
      cfs(3)=c2*aey*aez*cfs(17)
      cfs(4)=aey*aey*cfs(17)
      cfs(6)=aez*cfs(17)
      cfs(10)=c2*aey*cfs(17)
      return
   18 cfs(1)=aez*aez*aex*cfs(18)
      cfs(2)=aez*aez*cfs(18)
      cfs(4)=c2*aez*aex*cfs(18)
      cfs(7)=aex*cfs(18)
      cfs(9)=c2*aez*cfs(18)
      return
   19 cfs(1)=aez*aez*aey*cfs(19)
      cfs(3)=aez*aez*cfs(19)
      cfs(4)=c2*aez*aey*cfs(19)
      cfs(7)=aey*cfs(19)
      cfs(10)=c2*aez*cfs(19)
      return
   20 cfs(1)=aex*aey*aez*cfs(20)
      cfs(2)=aez*aey*cfs(20)
      cfs(3)=aex*aez*cfs(20)
      cfs(4)=aex*aey*cfs(20)
      cfs(8)=aez*cfs(20)
      cfs(9)=aey*cfs(20)
      cfs(10)=aex*cfs(20)
      return
      end
c
c
      subroutine dumdum(nuc,noc,vlist,c,dint,nro,nocmx,nt)
      implicit real*8(a-h,o-z)
      dimension vlist(nocmx,4), c(1), dint(nro,1)
      call nuc(noc,vlist,c,dint,nro,nocmx,nt)
      return
      end
c
c
      subroutine nonuc(noc,vlist,c,dint,nro,nocmx,nt)
c           for properties without nuclear part
      implicit real*8(a-h,o-z)
      dimension dint(nro,1)
      dint(1,1)=0.0d0
      nt = 1
      return
      end
c
c
      subroutine fldnuc(noc,vlist,c,dint,nro,nocmx,nt)
c           nuclear part of electric field
      implicit real*8(a-h,o-z)
      dimension dint(nro,1),vlist(nocmx ,4),c(3)  ,fnuc(3)
      do20 i=1,3
   20 fnuc(i)=0.0d0
      do21 i=1,noc
      dst=(vlist(i,1)-c(1))**2+(vlist(i,2)-c(2))**2+(vlist(i,3)-c(3))**2
      if(dst.eq.0.0d0) go to 21
      do 23 j=1,3
      fnuc(j)=fnuc(j)+vlist(i,4)*(vlist(i,j)-c(j))/(sqrt(dst))**3
   23 continue
   21 continue
      do22 l=1,3
   22 dint(1,l) = -fnuc(l)
      nt = 3
      return
      end
c
c
      subroutine potnuc(noc,vlist,c,dint,nro,nocmx,nt)
c                 nuclear part of potential
      implicit real*8(a-h,o-z)
      dimension dint(nro,1),vlist(nocmx ,4),c(3)
      pnuc=0.0d0
      do 12 i=1,noc
      dst=(vlist(i,1)-c(1))**2+(vlist(i,2)-c(2))**2+(vlist(i,3)-c(3))**2
      if(dst.eq.0.0d0)go to 12
      pnuc=pnuc+vlist(i,4)/sqrt(dst)
   12 continue
      dint(1,1)=pnuc
      nt = 1
      return
      end
c
c
      subroutine grdnuc(noc,vlist,c,dint,nro,nocmx,nt)
c           nuclear part of field gradient
      implicit real*8(a-h,o-z)
      dimension dint(nro,*),vlist(nocmx ,4),c(3)  ,gnuc(6)
      do 30 i=1,6
   30 gnuc(i)=0.0d0
      do 31 i=1,noc
      dst=(vlist(i,1)-c(1))**2+(vlist(i,2)-c(2))**2+(vlist(i,3)-c(3))**2
      if(dst.eq.0.0d0) go to 31
      mu=0
      do 32 j=1,3
      do 33 k=1,j
      mu=mu+1
      gnuc(mu)=gnuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))
     1/(sqrt(dst))**5
   33 continue
   32 continue
   31 continue
      dst=gnuc(1)+gnuc(3)+gnuc(6)
      dint(1,1)=3.0d0*gnuc(1)-dst
      dint(1,2)=3.0d0*gnuc(3)-dst
      dint(1,3)=3.0d0*gnuc(6)-dst
      dint(1,4)=gnuc(2)*3.0d0
      dint(1,5)=gnuc(4)*3.0d0
      dint(1,6)=gnuc(5)*3.0d0
      do 40 i=1,6
   40 dint(1,i) = - dint(1,i)
      nt = 6
      return
      end
c
c
      subroutine magnuc(noc,vlist,c,dint,nro,nocmx,nt)
c           nuclear part of second moment
      implicit real*8(a-h,o-z)
      dimension dint(nro,*),vlist(nocmx ,4),c(3)   ,qnuc(6)
      do50 i=1,6
   50 qnuc(i)=0.0d0
      do 51 i=1,noc
      mu=0
      do52 j=1,3
      do53 k=1,j
      mu=mu+1
      qnuc(mu)=qnuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))
   53 continue
   52 continue
   51 continue
      dint(1,1)=qnuc(1)
      dint(1,2)=qnuc(3)
      dint(1,3)=qnuc(6)
      dint(1,4)=qnuc(2)
      dint(1,5)=qnuc(4)
      dint(1,6)=qnuc(5)
      nt = 6
      return
      end
c
c
      subroutine qudnuc(noc,vlist,c,dint,nro,nocmx,nt)
c           nuclear part of quadropole moment
      implicit real*8(a-h,o-z)
      dimension dint(nro,*),vlist(nocmx ,*),c(3)   ,qnuc(6)
      data c3/3.0d0/
      do50 i=1,6
   50 qnuc(i)=0.0d0
      do 51 i=1,noc
      mu=0
      do52 j=1,3
      do53 k=1,j
      mu=mu+1
      qnuc(mu)=qnuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))
   53 continue
   52 continue
   51 continue
      dint(1,1)=qnuc(1)
      dint(1,2)=qnuc(3)
      dint(1,3)=qnuc(6)
      dint(1,4)=qnuc(2)*1.5d0
      dint(1,5)=qnuc(4)*1.5d0
      dint(1,6)=qnuc(5)*1.5d0
      dint(1,7)=dint(1,1)+dint(1,2)+dint(1,3)
      dint(1,1)=(c3*dint(1,1)-dint(1,7))/2.d0
      dint(1,2)=(c3*dint(1,2)-dint(1,7))/2.d0
      dint(1,3) =(c3*dint(1,3)-dint(1,7))/2.d0
      nt=7
      return
      end
c
c
      subroutine dipnuc(noc,vlist,c,dint,nro,nocmx,nt)
c           nuclear part of dipole moment
      implicit real*8(a-h,o-z)
      dimension dint(nro,*),vlist(nocmx ,4),c(3)  ,dnuc(3)
      do40 i=1,3
   40 dnuc(i)=0.0d0
      do41 i=1,noc
      do41 j=1,3
   41 dnuc(j)=dnuc(j)+vlist(i,4)*(vlist(i,j)-c(j))
      do42 k=1,3
   42 dint(1,k)=dnuc(k)
      nt = 3
      return
      end
c
c
      subroutine octnuc(noc,vlist,c,dint,nro,nocmx,nt)
c           nuclear part of third moment
      implicit real*8(a-h,o-z)
      dimension dint(nro,*),vlist(nocmx,4),c(*),onuc(10)
      do60 i=1,10
   60 onuc(i)=0.0d0
      do61 i=1,noc
      mu=0
      do62 j=1,3
      do63 k=1,j
      do64 m=1,k
      mu=mu+1
      onuc(mu)=onuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))*(
     1vlist(i,m)-c(m))
   64 continue
   63 continue
   62 continue
   61 continue
      dint(1,1)=onuc(1)
      dint(1,2)=onuc(4)
      dint(1,3)=onuc(10)
      dint(1,4)=onuc(2)
      dint(1,5)=onuc(5)
      dint(1,6)=onuc(3)
      dint(1,7)=onuc(7)
      dint(1,8)=onuc(8)
      dint(1,9)=onuc(9)
      dint(1,10)=onuc(6)
      nt = 10
      return
      end
c
c
      subroutine sldnuc(noc,vlist,c,dint,nro,nocmx,nt)
c           nuclear part of diamagnetic shielding
      implicit real*8(a-h,o-z)
      dimension dint(nro,*),vlist(nocmx,4),c(3),pnuc(6)
      do 30 i=1,6
   30 pnuc(i)=0.d0
      do 31 i=1,noc
      dst=(vlist(i,1)-c(1))**2+(vlist(i,2)-c(2))**2+(vlist(i,3)-c(3))**2
      if(dst.lt.1.d-10) go to 31
      mu=0
      do 32 j=1,3
      do 33 k=1,j
      mu=mu+1
      pnuc(mu)=pnuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))
     1/(sqrt(dst))**3
  33  continue
  32  continue
  31  continue
      dint(1,1)=pnuc(1)
      dint(1,2)=pnuc(3)
      dint(1,3)=pnuc(6)
      dint(1,4)=pnuc(2)
      dint(1,5)=pnuc(4)
      dint(1,6)=pnuc(5)
      dint(1,7)=dint(1,1)+dint(1,2)+dint(1,3)
      nt=7
      return
      end
c
c
      subroutine propa(name,ij,kl,eta,va,c,nt,nbmx)
      implicit real*8(a-h,o-z)
      common /abfunc/ ra(3),rb(3),ga,gb,ia,ib
c     this common-block is used in subroutine opap4 only
      dimension eta(nbmx,15)
      dimension va(1),c(1)
      dimension a(3),b(3),e(3),d(3),aa(20),bb(20),dd(84),v(10),val(10)
      dimension lin(84),min(84),nin(84)
      data lin/0,1,0,0,2,0,0,1,1,0,3,0,0,2,2,1,0,1,0,1,4,0,0,3,3,1,0,1,0
     1 ,2,2,0,2,1,1,5,0,0,3,3,2,2,0,0,4,4,1,0,0,1,1,3,1,2,2,1,6,0,0,3,
     2 3,0,5,5,1,0,0,1,4,4,2,0,2,0,3,3,1,2,2,1,4,1,1,2/,min/0,0,1,0,0,
     3 2,0,1,0,1,0,3,0,1,0,2,2,0,1,1,0,4,0,1,0,3,3,0,1,2,0,2,1,2,1,0,5,
     4 0,2,0,3,0,3,2,1,0,4,4,1,0,1,1,3,2,1,2,0,6,0,3,0,3,1,0,0,1,5,5,2,0
     5,0,2,4,4,2,1,3,1,3,2,1,4,1,2/,nin/0,0,0,1,0,0,2,0,1,1,0,0,3,0,1,0,
     6 1,2,2,1,0,0,4,0,1,0,1,3,3,0,2,2,1,1,2,0,0,5,0,2,0,3,2,3,0,1,0,1,
     7 4,4,3,1,1,1,2,2,0,0,6,0,3,3,0,1,5,5,1,0,0,2,4,4,0,2,1,2,2,3,1,3,
     8 1,1,4,2/
      do 10 i=1,3
      a(i)=eta(ij,i)
      b(i)=eta(kl,i)
      ra(i)=a(i)
      rb(i)=b(i)
   10 continue
      ga=eta(ij,4)
      gb=eta(kl,4)
      do 20 i=1,20
      aa(i)=eta(ij,i+5)
      bb(i)=eta(kl,i+5)
   20 continue
      iff1=eta(ij,5)+0.0001
      iff2=eta(kl,5)+0.0001
      ia=iff1
      ib=iff2
      call divpt(a,eta(ij,4),b,eta(kl,4),e,gama,efact)
      call rhftce(aa,a,e,iff1)
      call rhftce(bb,b,e,iff2)
      call prod(aa,bb,dd,iff1,iff2)
      do 30 i=1,nt
      val(i)=0.d0
   30 continue
      do 25 i=1,3
      d(i)=e(i)-c(i)
   25 continue
c           name represents an external function
      if(iff1.gt.10.or.iff2.gt.10) goto 110
      if(iff1.gt.4.or.iff2.gt.4) goto 120
      if(iff1.gt.1.or.iff2.gt.1) goto 130
c           s - s
      call name(lin(1),min(1),nin(1),gama,v,d)
      do 50 j=1,nt
      val(j)=dd(1)*v(j)+val(j)
   50 continue
      do 70 i=1,nt
      va(i)=efact*val(i)
   70 continue
      return
  130 continue
      if(iff1.ne.1.and.iff2.ne.1) goto 131
c           s - p
      do 61 i=1,4
      if(abs(dd(i))-1.d-8) 61,61,41
   41 call name(lin(i),min(i),nin(i),gama,v,d)
      do 51 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   51 continue
   61 continue
      do 71 i=1,nt
      va(i)=efact*val(i)
   71 continue
      return
  131 continue
c           p - p
      do 62 i=1,10
      if(abs(dd(i))-1.d-8) 62,62,42
   42 call name(lin(i),min(i),nin(i),gama,v,d)
      do 52 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   52 continue
   62 continue
      do 72 i=1,nt
      va(i)=efact*val(i)
   72 continue
      return
  120 continue
      if(iff1.ne.1.and.iff2.ne.1) goto 121
c           s - d
      do 63 i=1,10
      if(abs(dd(i))-1.d-8) 63,63,43
   43 call name(lin(i),min(i),nin(i),gama,v,d)
      do 53 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   53 continue
   63 continue
      do 73 i=1,nt
      va(i)=efact*val(i)
   73 continue
      return
  121 continue
      if(iff1.gt.4.and.iff2.gt.4) goto 122
c           p - d
      do 64 i=1,20
      if(abs(dd(i))-1.d-8) 64,64,44
   44 call name(lin(i),min(i),nin(i),gama,v,d)
      do 54 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   54 continue
   64 continue
      do 74 i=1,nt
      va(i)=efact*val(i)
   74 continue
      return
  122 continue
c           d - d
      do 65 i=1,35
      if(abs(dd(i))-1.d-8) 65,65,45
   45 call name(lin(i),min(i),nin(i),gama,v,d)
      do 55 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   55 continue
   65 continue
      do 75 i=1,nt
      va(i)=efact*val(i)
   75 continue
      return
  110 continue
      if(iff1.ne.1.and.iff2.ne.1) goto 111
c           s - f
      do 66 i=1,20
      if(abs(dd(i))-1.d-8) 66,66,46
   46 call name(lin(i),min(i),nin(i),gama,v,d)
      do 56 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   56 continue
   66 continue
      do 76 i=1,nt
      va(i)=efact*val(i)
   76 continue
      return
  111 continue
      if(iff1.gt.4.and.iff2.gt.4) goto 112
c           p - f
      do 67 i=1,35
      if(abs(dd(i))-1.d-8) 67,67,47
   47 call name(lin(i),min(i),nin(i),gama,v,d)
      do 57 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   57 continue
   67 continue
      do 77 i=1,nt
      va(i)=efact*val(i)
   77 continue
      return
  112 continue
      if(iff1.gt.10.and.iff2.gt.10) goto 113
c           d - f
      do 68 i=1,56
      if(abs(dd(i))-1.d-8) 68,68,48
   48 call name(lin(i),min(i),nin(i),gama,v,d)
      do 58 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   58 continue
   68 continue
      do 78 i=1,nt
      va(i)=efact*val(i)
   78 continue
      return
  113 continue
c           f - f
      do 69 i=1,84
      if(abs(dd(i))-1.d-8) 69,69,49
   49 call name(lin(i),min(i),nin(i),gama,v,d)
      do 59 j=1,nt
      val(j)=dd(i)*v(j)+val(j)
   59 continue
   69 continue
      do 79 i=1,nt
      va(i)=efact*val(i)
   79 continue
      return
      end
c
c
      real*8 function olap(l,m,n,gama)
      implicit real*8(a-h,o-z)
      dimension dftr(7)
      data pi /3.14159265359d0/
      data dftr /1.d0,1.d0,3.d0,15.d0,105.d0,945.d0,10395.d0/
      lh=l/2
      mh=m/2
      nh=n/2
      if(2*lh-l) 50,1,50
   1  if(2*mh-m) 50,2,50
   2  if(2*nh-n) 50,3,50
   3  hf=.5d0
      olap=(sqrt(pi/gama))**3*(hf/gama)**(lh+mh+nh)*dftr(lh+1)*dftr(mh+
     11)*dftr(nh+1)
      return
  50  olap=0.d0
      return
      end
c
c
      subroutine prod(c,d,s,iff1,iff2)
      implicit real*8(a-h,o-z)
      dimension c(20),d(20),s(84)
      do 10 i=1,84
   10 s(i)=0.d0
      if(iff1.gt.10.or.iff2.gt.10) goto 30
      if(iff1.gt.4.or.iff2.gt.4) goto 20
      s( 1)=c( 1)*d( 1)
c              end of s - s
      if(iff1.eq.1.and.iff2.eq.1) return
      s( 2)=c( 1)*d( 2)+c( 2)*d( 1)
      s( 3)=c( 1)*d( 3)+c( 3)*d( 1)
      s( 4)=c( 1)*d( 4)+c( 4)*d( 1)
c              end of s - p
      if(iff1.eq.1.or.iff2.eq.1) return
      s( 5)=c( 2)*d( 2)
      s( 6)=c( 3)*d( 3)
      s( 7)=c( 4)*d( 4)
      s( 8)=c( 2)*d( 3)+c( 3)*d( 2)
      s( 9)=c( 2)*d( 4)+c( 4)*d( 2)
      s(10)=c( 3)*d( 4)+c( 4)*d( 3)
c              end of p - p
      return
   20 continue
      s( 1)=c( 1)*d( 1)
      s( 2)=c( 1)*d( 2)+c( 2)*d( 1)
      s( 3)=c( 1)*d( 3)+c( 3)*d( 1)
      s( 4)=c( 1)*d( 4)+c( 4)*d( 1)
      s( 5)=c( 1)*d( 5)+c( 5)*d( 1)+c( 2)*d( 2)
      s( 6)=c( 1)*d( 6)+c( 6)*d( 1)+c( 3)*d( 3)
      s( 7)=c( 1)*d( 7)+c( 7)*d( 1)+c( 4)*d( 4)
      s( 8)=c( 1)*d( 8)+c( 8)*d( 1)+c( 2)*d( 3)+c( 3)*d( 2)
      s( 9)=c( 1)*d( 9)+c( 9)*d( 1)+c( 2)*d( 4)+c( 4)*d( 2)
      s(10)=c( 1)*d(10)+c(10)*d( 1)+c( 3)*d( 4)+c( 4)*d( 3)
c              end of s - d
      if(iff1.eq.1.or.iff2.eq.1) return
      s(11)=c( 2)*d( 5)+c( 5)*d( 2)
      s(12)=c( 3)*d( 6)+c( 6)*d( 3)
      s(13)=c( 4)*d( 7)+c( 7)*d( 4)
      s(14)=c( 2)*d( 8)+c( 8)*d( 2)+c( 3)*d( 5)+c( 5)*d( 3)
      s(15)=c( 2)*d( 9)+c( 9)*d( 2)+c( 4)*d( 5)+c( 5)*d( 4)
      s(16)=c( 2)*d( 6)+c( 6)*d( 2)+c( 3)*d( 8)+c( 8)*d( 3)
      s(17)=c( 3)*d(10)+c(10)*d( 3)+c( 4)*d( 6)+c( 6)*d( 4)
      s(18)=c( 2)*d( 7)+c( 7)*d( 2)+c( 4)*d( 9)+c( 9)*d( 4)
      s(19)=c( 3)*d( 7)+c( 7)*d( 3)+c( 4)*d(10)+c(10)*d( 4)
      s(20)=c( 2)*d(10)+c(10)*d( 2)+c( 3)*d( 9)
     1     +c( 9)*d( 3)+c( 4)*d( 8)+c( 8)*d( 4)
c              end of p - d
      if(iff1.lt.5.or.iff2.lt.5) return
      s(21)=c( 5)*d( 5)
      s(22)=c( 6)*d( 6)
      s(23)=c( 7)*d( 7)
      s(24)=c( 5)*d( 8)+c( 8)*d( 5)
      s(25)=c( 5)*d( 9)+c( 9)*d( 5)
      s(26)=c( 6)*d( 8)+c( 8)*d( 6)
      s(27)=c( 6)*d(10)+c(10)*d( 6)
      s(28)=c( 7)*d( 9)+c( 9)*d( 7)
      s(29)=c( 7)*d(10)+c(10)*d( 7)
      s(30)=c( 5)*d( 6)+c( 6)*d( 5)+c( 8)*d( 8)
      s(31)=c( 5)*d( 7)+c( 7)*d( 5)+c( 9)*d( 9)
      s(32)=c( 6)*d( 7)+c( 7)*d( 6)+c(10)*d(10)
      s(33)=c( 5)*d(10)+c(10)*d( 5)+c( 8)*d( 9)+c( 9)*d( 8)
      s(34)=c( 6)*d( 9)+c( 9)*d( 6)+c( 8)*d(10)+c(10)*d( 8)
      s(35)=c( 7)*d( 8)+c( 8)*d( 7)+c( 9)*d(10)+c(10)*d( 9)
c              end of d - d
      return
   30 continue
      s( 1)=c( 1)*d( 1)
      s( 2)=c( 1)*d( 2)+c( 2)*d( 1)
      s( 3)=c( 1)*d( 3)+c( 3)*d( 1)
      s( 4)=c( 1)*d( 4)+c( 4)*d( 1)
      s( 5)=c( 1)*d( 5)+c( 5)*d( 1)+c( 2)*d( 2)
      s( 6)=c( 1)*d( 6)+c( 6)*d( 1)+c( 3)*d( 3)
      s( 7)=c( 1)*d( 7)+c( 7)*d( 1)+c( 4)*d( 4)
      s( 8)=c( 1)*d( 8)+c( 8)*d( 1)+c( 2)*d( 3)+c( 3)*d( 2)
      s( 9)=c( 1)*d( 9)+c( 9)*d( 1)+c( 2)*d( 4)+c( 4)*d( 2)
      s(10)=c( 1)*d(10)+c(10)*d( 1)+c( 3)*d( 4)+c( 4)*d( 3)
      s(11)=c( 2)*d( 5)+c( 5)*d( 2)+c(1)*d(11)+c(11)*d(1)
      s(12)=c( 3)*d( 6)+c( 6)*d( 3)+c(1)*d(12)+c(12)*d(1)
      s(13)=c( 4)*d( 7)+c( 7)*d( 4)+c(1)*d(13)+c(13)*d(1)
      s(14)=c( 2)*d( 8)+c( 8)*d( 2)+c( 3)*d( 5)+c( 5)*d( 3)+c(1)*d(14)+
     1      c(14)*d(1)
      s(15)=c( 2)*d( 9)+c( 9)*d( 2)+c( 4)*d( 5)+c( 5)*d( 4)+c(1)*d(15)+
     1      c(15)*d(1)
      s(16)=c( 2)*d( 6)+c( 6)*d( 2)+c( 3)*d( 8)+c( 8)*d( 3)+c(1)*d(16)+
     1      c(16)*d(1)
      s(17)=c( 3)*d(10)+c(10)*d( 3)+c( 4)*d( 6)+c( 6)*d( 4)+c(1)*d(17)+
     1      c(17)*d(1)
      s(18)=c( 2)*d( 7)+c( 7)*d( 2)+c( 4)*d( 9)+c( 9)*d( 4)+c(1)*d(18)+
     1      c(18)*d(1)
      s(19)=c( 3)*d( 7)+c( 7)*d( 3)+c( 4)*d(10)+c(10)*d( 4)+c(1)*d(19)+
     1      c(19)*d(1)
      s(20)=c( 2)*d(10)+c(10)*d( 2)+c( 3)*d( 9)+c(9)*d(3)+c(4)*d(8)+
     1      c(8)*d(4)+c(1)*d(20)+c(20)*d(1)
c              end of s - f
      if(iff1.eq.1.or.iff2.eq.1) return
      s(21)=c( 5)*d( 5)+c(2)*d(11)+c(11)*d(2)
      s(22)=c( 6)*d( 6)+c(3)*d(12)+c(12)*d(3)
      s(23)=c( 7)*d( 7)+c(4)*d(13)+c(13)*d(4)
      s(24)=c( 5)*d( 8)+c( 8)*d( 5)+c(3)*d(11)+c(11)*d(3)+c(2)*d(14)+
     1      c(14)*d(2)
      s(25)=c( 5)*d( 9)+c( 9)*d( 5)+c(2)*d(15)+c(15)*d(2)+c(4)*d(11)+
     1      c(11)*d(4)
      s(26)=c( 6)*d( 8)+c( 8)*d( 6)+c(2)*d(12)+c(12)*d(2)+c(3)*d(16)+
     1      c(16)*d(3)
      s(27)=c( 6)*d(10)+c(10)*d( 6)+c(3)*d(17)+c(17)*d(3)+c(4)*d(12)+
     1      c(12)*d(4)
      s(28)=c( 7)*d( 9)+c( 9)*d( 7)+c(2)*d(13)+c(13)*d(2)+c(4)*d(18)+
     1      c(18)*d(4)
      s(29)=c( 7)*d(10)+c(10)*d( 7)+c(3)*d(13)+c(13)*d(3)+c(4)*d(19)+
     1      c(19)*d(4)
      s(30)=c( 5)*d( 6)+c( 6)*d( 5)+c( 8)*d( 8)+c(2)*d(16)+c(16)*d(2)+
     1      c(3)*d(14)+c(14)*d(3)
      s(31)=c( 5)*d( 7)+c( 7)*d( 5)+c( 9)*d( 9)+c(2)*d(18)+c(18)*d(2)+
     1      c(4)*d(15)+c(15)*d(4)
      s(32)=c( 6)*d( 7)+c( 7)*d( 6)+c(10)*d(10)+c(3)*d(19)+c(19)*d(3)+
     1      c(4)*d(17)+c(17)*d(4)
      s(33)=c( 5)*d(10)+c(10)*d( 5)+c( 8)*d( 9)+c( 9)*d( 8)+c(3)*d(15)+
     1      c(15)*d(3)+c(4)*d(14)+c(14)*d(4)+c(2)*d(20)+c(20)*d(2)
      s(34)=c( 6)*d( 9)+c( 9)*d( 6)+c( 8)*d(10)+c(10)*d( 8)+c(2)*d(17)+
     1      d(2)*c(17)+c(3)*d(20)+c(20)*d(3)+c(4)*d(16)+c(16)*d(4)
      s(35)=c( 7)*d( 8)+c( 8)*d( 7)+c( 9)*d(10)+c(10)*d( 9)+c(2)*d(19)+
     1      c(19)*d(2)+c(3)*d(18)+c(18)*d(3)+c(4)*d(20)+c(20)*d(4)
c              end of p - f
      if(iff1.eq.2.or.iff2.eq.2) return
      s(36)=c(5)*d(11)+c(11)*d(5)
      s(37)=c(6)*d(12)+c(12)*d(6)
      s(38)=c(7)*d(13)+c(13)*d(7)
      s(39)=c(6)*d(11)+c(11)*d(6)+c(5)*d(16)+c(16)*d(5)+c(8)*d(14)+
     1      c(14)*d(8)
      s(40)=c(7)*d(11)+c(11)*d(7)+c(5)*d(18)+c(18)*d(5)+c(9)*d(15)+
     1      c(15)*d(9)
      s(41)=c(5)*d(12)+c(12)*d(5)+c(6)*d(14)+c(14)*d(6)+c(8)*d(16)+
     1      c(16)*d(8)
      s(42)=c(5)*d(13)+c(13)*d(5)+c(7)*d(15)+c(15)*d(7)+c(9)*d(18)+
     1      c(18)*d(9)
      s(43)=c(7)*d(12)+c(12)*d(7)+c(6)*d(19)+c(19)*d(6)+c(10)*d(17)+
     1      c(17)*d(10)
      s(44)=c(6)*d(13)+c(13)*d(6)+c(7)*d(17)+c(17)*d(7)+c(10)*d(19)+
     1      c(19)*d(10)
      s(45)=c(8)*d(11)+c(11)*d(8)+c(5)*d(14)+c(14)*d(5)
      s(46)=c(9)*d(11)+c(11)*d(9)+c(5)*d(15)+c(15)*d(5)
      s(47)=c(8)*d(12)+c(12)*d(8)+c(6)*d(16)+c(16)*d(6)
      s(48)=c(10)*d(12)+c(12)*d(10)+c(6)*d(17)+c(17)*d(6)
      s(49)=c(10)*d(13)+c(13)*d(10)+c(7)*d(19)+c(19)*d(7)
      s(50)=c(9)*d(13)+c(13)*d(9)+c(7)*d(18)+c(18)*d(7)
      s(51)=c(8)*d(13)+c(13)*d(8)+c(7)*d(20)+c(20)*d(7)+c(9)*d(19)+
     1      c(19)*d(9)+c(10)*d(18)+c(18)*d(10)
      s(52)=c(10)*d(11)+c(11)*d(10)+c(5)*d(20)+c(20)*d(5)+c(9)*d(14)+
     1      c(14)*d(9)+c(8)*d(15)+c(15)*d(8)
      s(53)=c(9)*d(12)+c(12)*d(9)+c(6)*d(20)+c(20)*d(6)+c(10)*d(16)+
     1      c(16)*d(10)+c(8)*d(17)+c(17)*d(8)
      s(54)=c(5)*d(17)+c(17)*d(5)+c(6)*d(15)+c(15)*d(6)+c(14)*d(10)+
     1      d(14)*c(10)+c(9)*d(16)+c(16)*d(9)+c(8)*d(20)+c(20)*d(8)
      s(55)=c(5)*d(19)+c(19)*d(5)+c(7)*d(14)+c(14)*d(7)+c(10)*d(15)+
     1      c(15)*d(10)+c(9)*d(20)+c(20)*d(9)+c(8)*d(18)+c(18)*d(8)
      s(56)=c(6)*d(18)+c(18)*d(6)+c(7)*d(16)+c(16)*d(7)+c(10)*d(20)+
     1      c(20)*d(10)+c(9)*d(17)+c(17)*d(9)+c(8)*d(19)+c(19)*d(8)
c              end of d - f
      if(iff1.eq.3.or.iff2.eq.3) return
      s(57)=c(11)*d(11)
      s(58)=c(12)*d(12)
      s(59)=c(13)*d(13)
      s(60)=c(11)*d(12)+c(12)*d(11)+c(14)*d(16)+c(16)*d(14)
      s(61)=c(11)*d(13)+c(13)*d(11)+c(15)*d(18)+c(18)*d(15)
      s(62)=c(12)*d(13)+c(13)*d(12)+c(17)*d(19)+c(19)*d(17)
      s(63)=c(11)*d(14)+c(14)*d(11)
      s(64)=c(11)*d(15)+c(15)*d(11)
      s(65)=c(13)*d(18)+c(18)*d(13)
      s(66)=c(13)*d(19)+c(19)*d(13)
      s(67)=c(12)*d(17)+d(12)*c(17)
      s(68)=c(12)*d(16)+c(16)*d(12)
      s(69)=c(11)*d(16)+c(16)*d(11)+c(14)*d(14)
      s(70)=c(11)*d(18)+c(18)*d(11)+c(15)*d(15)
      s(71)=c(13)*d(15)+c(15)*d(13)+c(18)*d(18)
      s(72)=c(13)*d(17)+c(17)*d(13)+c(19)*d(19)
      s(73)=c(12)*d(14)+c(14)*d(12)+c(16)*d(16)
      s(74)=c(12)*d(19)+c(19)*d(12)+c(17)*d(17)
      s(75)=c(11)*d(17)+c(17)*d(11)+c(14)*d(20)+c(20)*d(14)+
     1      c(15)*d(16)+c(16)*d(15)
      s(76)=c(11)*d(19)+c(19)*d(11)+c(20)*d(15)+c(15)*d(20)+
     1      c(14)*d(18)+c(18)*d(14)
      s(77)=c(12)*d(18)+c(18)*d(12)+c(16)*d(19)+c(19)*d(16)+
     1      c(17)*d(20)+c(20)*d(17)
      s(78)=c(13)*d(14)+c(14)*d(13)+c(15)*d(19)+c(19)*d(15)+
     1      c(18)*d(20)+c(20)*d(18)
      s(79)=c(12)*d(15)+c(15)*d(12)+c(14)*d(17)+c(17)*d(14)+
     1      c(16)*d(20)+c(20)*d(16)
      s(80)=c(13)*d(16)+c(16)*d(13)+c(17)*d(18)+c(18)*d(17)+
     1      c(19)*d(20)+c(20)*d(19)
      s(81)=c(11)*d(20)+c(20)*d(11)+c(14)*d(15)+c(15)*d(14)
      s(82)=c(12)*d(20)+c(20)*d(12)+c(16)*d(17)+c(17)*d(16)
      s(83)=c(13)*d(20)+c(20)*d(13)+c(18)*d(19)+c(19)*d(18)
      s(84)=c(14)*d(19)+c(19)*d(14)+c(15)*d(17)+c(17)*d(15)+
     1      c(16)*d(18)+c(18)*d(16)+c(20)*d(20)
c              end of f - f
      return
      end
c
c
      subroutine fmc(mvar,xvar,expmx,fmch)
      implicit real*8(a-h,o-z)
       m=mvar
       x=xvar
      expmx=exp(-x)
      if(x -20.d0) 10,10,20
c         x less than or equal to 20.0.
   10  a=m
       a=a+0.5d0
      term=1.d0/a
       ptlsum=term
      do 11 i=2,60
       a=a+1.d0
       term=term*x/a
       ptlsum=ptlsum+term
       if (term/ptlsum- 1.d-8    ) 12, 11, 11
   11  continue
       print 999,m,x
  12  fmch=.5d0*ptlsum*expmx
       return
c         x greater than 20.0.
   20  a=m
      xd=1.d0/x
       b=a+0.5d0
       a=a-0.5d0
      approx=.886226925428d0*(sqrt(xd)*xd**m)
       if (m) 21, 23, 21
   21  do 22 i=1,m
       b=b-1.d0
   22  approx=approx*b
  23  fimult=.5d0*expmx*xd
       fiprop=fimult/approx
       term=1.d0
       ptlsum=term
       notrms=x
       notrms=notrms+m
       do 24 i=2,notrms
       term=term*a*xd
       ptlsum=ptlsum+term
      if (abs(term*fiprop/ptlsum) -  1.d-8    ) 25,25,24
   24  a=a-1.d0
       print 999,m,x
   25  fmch=approx-fimult*ptlsum
       return
  999  format (24h no convergence for fmch, i6, d16.9)
      end
c
c
      subroutine opaa0(l,m,n,ga,v,d)
c           electronic part of potential
      implicit real*8(a-h,o-z)
      dimension v(1),d(3),fnu(7),fn(7),fd(7),dp(4)
      data pi /3.14159265359d0/, fn /1.d0,1.d0,2.d0,6.d0,24.d0,120.d0,
     1 720.d0/
      data fd/1.d0,1.d0,0.5d0,0.1666666666667d0,0.416666666667d-1,
     1 .833333333333d-2,.1388888888889d-2/
      itot=l+m+n
      itoth=itot/2
      pre=(2.d0*pi/ga)*fn(l+1)*fn(m+1)*fn(n+1)
      if(2*itoth-itot) 1,2,1
   1  pre=-pre
   2  del=.25d0/ga
      x=ga*(d(1)**2+d(2)**2+d(3)**2)
      xx=2.d0*x
      call fmc(6,x,expmx,fmch)
      fnu(7)=fmch
      fnu(6)=(expmx+xx*fnu(7))/11.d0
      fnu(5)=(expmx+xx*fnu(6))/9.d0
      fnu(4)=(expmx+xx*fnu(5))/7.d0
      fnu(3)=(expmx+xx*fnu(4))/5.d0
      fnu(2)=(expmx+xx*fnu(3))/3.d0
      fnu(1)=expmx+xx*fnu(2)
      dp(1)=1.d0
      dp(2)=del
      dp(3)=del**2
      v(1)=pre*aainer(0,0,0,l,m,n,d,dp,fnu,fn,fd)*(-1.0d0)
      return
      end
c
c
      subroutine opaa1(l,m,n,ga,v,d)
c           electronic part of electric field
      implicit real*8(a-h,o-z)
      dimension v(3),d(3),fnu(9),fd(9),fn(9),dp(4)
      data pi /3.14159265359d0/, fn /1.d0,1.d0,2.d0,6.d0,24.d0,120.d0,
     1 720.d0,5040.d0,40320.d0/
      data fd /1.d0,1.d0,.5d0,.1666666666667d0,.416666666667d-1,
     1 .833333333333d-2,.1388888888889d-2,.1984126984d-3,0.248015873d-4/
      itot=l+m+n
      itoth=itot/2
      pre=4.d0*pi*fn(l+1)*fn(m+1)*fn(n+1)
       if(2*itoth-itot) 1,2,1
   1  pre=-pre
   2  del=.25d0/ga
      x=ga*(d(1)**2+d(2)**2+d(3)**2)
      xx=2.d0*x
      call fmc(8,x,expmx,fmch)
      fnu(9)=fmch
      fnu(8)=(expmx+xx*fnu(9))/15.d0
      fnu(7)=(expmx+xx*fnu(8))/13.d0
      fnu(6)=(expmx+xx*fnu(7))/11.d0
      fnu(5)=(expmx+xx*fnu(6))/9.d0
      fnu(4)=(expmx+xx*fnu(5))/7.d0
      fnu(3)=(expmx+xx*fnu(4))/5.d0
      fnu(2)=(expmx+xx*fnu(3))/3.d0
      fnu(1)=expmx+xx*fnu(2)
      dp(1)=1.d0
      dp(2)=del
      dp(3)=del**2
      dp(4)=del*dp(3)
      v(1)=pre*aainer(1,0,0,l,m,n,d,dp,fnu,fn,fd)
      v(2)=pre*aainer(0,1,0,l,m,n,d,dp,fnu,fn,fd)
      v(3)=pre*aainer(0,0,1,l,m,n,d,dp,fnu,fn,fd)
      return
      end
c
c
      subroutine opaa2(l,m,n,ga,v,d)
c           electronic part of field gradient
      implicit real*8(a-h,o-z)
      dimension v(6),d(3),fnu(9),fd(9),fn(9),dp(5)
      data pi /3.14159265359d0/, fn /1.d0,1.d0,2.d0,6.d0,24.d0,120.d0,
     1 720.d0,5040.d0,40320.d0/
      data fd /1.d0,1.d0,.5d0,.1666666666667d0,.416666666667d-1,
     1 .833333333333d-2,.1388888888889d-2,.1984126984d-3,0.248015873d-4/
      itot=l+m+n
      itoth=itot/2
      pre=(8.d0*pi*ga)*fn(l+1)*fn(m+1)*fn(n+1)
      if(2*itoth-itot)1,2,1
   1  pre=-pre
   2  del=.25d0/ga
      x=ga*(d(1)**2+d(2)**2+d(3)**2)
      xx=2.d0*x
      call fmc(8,x,expmx,fmch)
      fnu(9)=fmch
      fnu(8)=(expmx+xx*fnu(9))/15.d0
      fnu(7)=(expmx+xx*fnu(8))/13.d0
      fnu(6)=(expmx+xx*fnu(7))/11.d0
      fnu(5)=(expmx+xx*fnu(6))/9.d0
      fnu(4)=(expmx+xx*fnu(5))/7.d0
      fnu(3)=(expmx+xx*fnu(4))/5.d0
      fnu(2)=(expmx+xx*fnu(3))/3.d0
      fnu(1)=expmx+xx*fnu(2)
      dp(1)=1.d0
      dp(2)=del
      dp(3)=del**2
      dp(4)=del*dp(3)
      dp(5)=del*dp(4)
      v(1)=pre*aainer(2,0,0,l,m,n,d,dp,fnu,fn,fd)
      v(2)=pre*aainer(0,2,0,l,m,n,d,dp,fnu,fn,fd)
      v(3)=pre*aainer(0,0,2,l,m,n,d,dp,fnu,fn,fd)
      v(4)=pre*aainer(1,1,0,l,m,n,d,dp,fnu,fn,fd)
      v(5)=pre*aainer(1,0,1,l,m,n,d,dp,fnu,fn,fd)
      v(6)=pre*aainer(0,1,1,l,m,n,d,dp,fnu,fn,fd)
      call opac3(l,m,n,ga,dp(1),d)
      dp(1)=1.333333333333d0*pi*dp(1)
      v(1)=v(1)+dp(1)
      v(2)=v(2)+dp(1)
      v(3)=v(3)+dp(1)
      return
      end
c
c
      subroutine opaa3(l,m,n,ga,v,d)
c           electronic part of diamagnetic shielding
      implicit real*8(a-h,o-z)
      dimension vx(3),d(3),v(7)
      call opaa1(l,m,n,ga,vx,d)
      v(1)=vx(1)*d(1)
      v(2)=vx(2)*d(2)
      v(3)=vx(3)*d(3)
      v(4)=vx(2)*d(1)
      v(5)=vx(3)*d(1)
      v(6)=vx(3)*d(2)
      call opaa1(l+1,m,n,ga,vx,d)
      v(1)=v(1)+vx(1)
      v(4)=v(4)+vx(2)
      v(5)=v(5)+vx(3)
      call opaa1(l,m+1,n,ga,vx,d)
      v(2)=v(2)+vx(2)
      v(6)=v(6)+vx(3)
      call opaa1(l,m,n+1,ga,vx,d)
      v(3)=v(3)+vx(3)
      v(7)=v(1)+v(2)+v(3)
      do 300 i=1,7
      v(i)=-v(i)
 300  continue
      return
      end
c
c
      subroutine opab1(l,m,n,ga,v,d)
c           electronic part of dipole moment
      implicit real*8(a-h,o-z)
      dimension v(3),d(3),g(4)
      g(1)=olap(l,m,n,ga)
      g(2)=olap(l+1,m,n,ga)
      g(3)=olap(l,m+1,n,ga)
      g(4)=olap(l,m,n+1,ga)
      v(1)=g(2)+d(1)*g(1)
      v(2)=g(3)+d(2)*g(1)
      v(3)=g(4)+d(3)*g(1)
      do 300 i=1,3
  300 v(i)=-v(i)
      return
      end
c
c
      subroutine opab2(l,m,n,ga,v,d)
c           electronic part of quadrupole moment
      implicit real*8(a-h,o-z)
      dimension v(7),d(3),g(10)
      g(1)=olap(l,m,n,ga)
      g(2)=olap(l+1,m,n,ga)
      g(3)=olap(l,m+1,n,ga)
      g(4)=olap(l,m,n+1,ga)
      g(5)=olap(l+2,m,n,ga)
      g(6)=olap(l,m+2,n,ga)
      g(7)=olap(l,m,n+2,ga)
      g(8)=olap(l+1,m+1,n,ga)
      g(9)=olap(l+1,m,n+1,ga)
      g(10)=olap(l,m+1,n+1,ga)
      v(1)=g(5)+2.d0*d(1)*g(2)+d(1)**2*g(1)
      v(2)=g(6)+2.d0*d(2)*g(3)+d(2)**2*g(1)
      v(3)=g(7)+2.d0*d(3)*g(4)+d(3)**2*g(1)
      v(4)=(g(8)+d(1)*g(3)+d(2)*g(2)+d(1)*d(2)*g(1))/2.d0
      v(5)=(g(9)+d(1)*g(4)+ d(3)*g(2)+d(1)*d(3)*g(1))/2.d0
      v(6)=(g(10)+d(2)*g(4)+d(3)*g(3)+d(2)*d(3)*g(1))/2.d0
      v(7)=v(1)+v(2)+v(3)
      v(1)=(3.d0*v(1)-v(7))/2.d0
      v(2)=(3.d0*v(2)-v(7))/2.d0
      v(3)=(3.d0*v(3)-v(7))/2.d0
      v(4)=3.d0*v(4)
      v(5)=3.d0*v(5)
      v(6)=3.d0*v(6)
      do 300 i=1,7
 300  v(i)=-v(i)
      return
      end
c
c
      subroutine opab3(l,m,n,ga,v,dd)
c           electronic part of third moment
      implicit real*8(a-h,o-z)
      dimension v(10),dd(3),g(20),d(20)
      g(1)=olap( l,m,n,ga)
      g(2)=olap(l+1,m,n,ga)
      g(3)=olap(l,m+1,n,ga)
      g(4)=olap(l,m,n+1,ga)
      g(5)=olap(l+2,m,n,ga)
      g(6)=olap(l,m+2,n,ga)
      g(7)=olap(l,m,n+2,ga)
      g(8)=olap(l+1,m+1,n,ga)
      g(9)=olap(l+1,m,n+1,ga)
      g(10)=olap(l,m+1,n+1,ga)
      g(11)=olap(l+3,m,n,ga)
      g(12)=olap(l,m+3,n,ga)
      g(13)=olap(l,m,n+3,ga)
      g(14)=olap(l+2,m+1,n,ga)
      g(15)=olap(l+2,m,n+1,ga)
      g(16)=olap(l+1,m+2,n,ga)
      g(17)=olap(l,m+2,n+1,ga)
      g(18)=olap(l+1,m,n+2,ga)
      g(19)=olap(l,m+1,n+2,ga)
      g(20)=olap(l+1,m+1,n+1,ga)
      d(2)=dd(1)
      d(3)=dd(2)
      d(4)=dd(3)
      d(5)=dd(1)**2
      d(6)=dd(2)**2
      d(7)=dd(3)**2
      d(8)=d(2)*d(3)
      d(9)=d(2)*d(4)
      d(10)=d(3)*d(4)
      d(11)=d(5)*d(2)
      d(12)=d(6)*d(3)
      d(13)=d(7)*d(4)
      d(14)=d(5)*d(3)
      d(15)=d(5)*d(4)
      d(16)=d(6)*d(2)
      d(17)=d(6)*d(4)
      d(18)=d(7)*d(2)
      d(19)=d(7)*d(3)
      d(20)=d(8)*d(4)
      v(1)=g(11)+3.d0*d(2)*g(5)+3.d0*d(5)*g(2)+d(11)*g(1)
      v(2)=g(12)+3.d0*d(3)*g(6)+3.d0*d(6)*g(3)+d(12)*g(1)
      v(3)=g(13)+3.d0*d(4)*g(7)+3.d0*d(7)*g(4)+d(13)*g(1)
      v(4)=g(14)+2.d0*d(2)*g(8)+d(3)*g(5)+d(5)*g(3)+2.d0*d(8)*g(2)
     1 +d(14)*g(1)
      v(5)=g(15)+2.d0*d(2)*g(9)+d(4)*g(5)+d(5)*g(4)+2.d0*d(9)*g(2)
     1 +d(15)*g(1)
      v(6)=g(16)+2.d0*d(3)*g(8)+d(2)*g(6)+d(6)*g(2)+2.d0*d(8)*g(3)
     1 +d(16)*g(1)
      v(7)=g(17)+2.d0*d(3)*g(10)+d(4)*g(6)+d(6)*g(4)+2.d0*d(10)*g(3)
     1 +d(17)*g(1)
      v(8)=g(18)+2.d0*d(4)*g(9)+d(2)*g(7)+d(7)*g(2)+2.d0*d(9)*g(4)
     1 +d(18)*g(1)
      v(9)=g(19)+2.d0*d(4)*g(10)+d(3)*g(7)+d(7)*g(3)+2.d0*d(10)*g(4)
     1 +d(19)*g(1)
      v(10)=g(20)+d(2)*g(10)+d(3)*g(9)+d(4)*g(8)+d(10)*g(2)+d(9)*g(3)+
     1d(8)*g(4)+d(20)*g(1)
      do 300 i=1,10
 300  v(i)=-v(i)
      return
      end
c
c
      subroutine opab4(l,m,n,ga,v,d)
c           electronic part of second moment
      implicit real*8(a-h,o-z)
      dimension v(6),d(3),g(10)
      g(1)=olap(l,m,n,ga)
      g(2)=olap(l+1,m,n,ga)
      g(3)=olap(l,m+1,n,ga)
      g(4)=olap(l,m,n+1,ga)
      g(5)=olap(l+2,m,n,ga)
      g(6)=olap(l,m+2,n,ga)
      g(7)=olap(l,m,n+2,ga)
      g(8)=olap(l+1,m+1,n,ga)
      g(9)=olap(l+1,m,n+1,ga)
      g(10)=olap(l,m+1,n+1,ga)
      v(1)=g(5)+2.d0*d(1)*g(2)+d(1)**2*g(1)
      v(2)=g(6)+2.d0*d(2)*g(3)+d(2)**2*g(1)
      v(3)=g(7)+2.d0*d(3)*g(4)+d(3)**2*g(1)
      v(4)=g(8)+d(1)*g(3)+d(2)*g(2)+d(1)*d(2)*g(1)
      v(5)=g(9)+d(1)*g(4)+d(3)*g(2)+d(1)*d(3)*g(1)
      v(6)=g(10)+d(2)*g(4)+d(3)*g(3)+d(2)*d(3)*g(1)
      do 300 i=1,6
 300  v(i)=-v(i)
      return
      end
c
c
      subroutine opac1(l,m,n,ga,v,d)
c           electronic part of planar density
      implicit real*8(a-h,o-z)
      dimension v(3),d(3),dd(3),dftr(6)
      data pi /3.14159265359d0/, dftr /1.d0,1.d0,3.d0,5.d0,15.d0,105.d0/
       do 1 i=1,3
      dd(i)=-d(i)
      v(i)=0.d0
   1  continue
      g=.5d0/ga
      lh=l/2
      mh=m/2
      nh=n/2
      pre=2.d0*g*pi
      if(2*nh-n) 9,3,9
   3  if(2*mh-m) 7,4,7
   4  v(1)=1.d0
      if(l) 39,41,39
  39  v(1)=dd(1)**l
  41  v(1)=pre*v(1)*   exp (-ga*dd(1)**2)*g**(mh+nh)*dftr(mh+1)*dftr(nh
     1+1)
      if(2*lh-l) 50,5,50
   5  v(2)=1.d0
      if(m) 49,51,49
  49  v(2)=dd(2)**m
  51  v(2)=pre*v(2)*   exp (-ga*dd(2)**2)*g**(lh+nh)*dftr(lh+1)*dftr(nh
     1+1)
   6  v(3)=1.d0
      if(n) 59,61,59
  59  v(3)=dd(3)**n
  61  v(3)=pre*v(3)*   exp (-ga*dd(3)**2)*g**(lh+mh)*dftr(lh+1)*dftr(mh
     1+1)
      return
   7  if(2*lh-l) 50,8,50
   8  v(2)=1.d0
      if(m) 79,81,79
  79  v(2)=dd(2)**m
  81  v(2)=pre*v(2)*   exp (-ga*dd(2)**2)*g**(lh+nh)*dftr(lh+1)*dftr(nh
     1+1)
      return
   9  if(2*mh-m) 50,10,50
  10  if(2*lh-l) 50,6,50
  50  return
      end
c
c
      subroutine opac2(l,m,n,ga,v,d)
c           electronic part of linear density
      implicit real*8(a-h,o-z)
      dimension v(3),d(3),dd(3),ds(3),dftr(5)
      data sqrpi /1.772453850906d0/, dftr /1.d0,1.d0,3.d0,15.d0,105.d0/
      do 1 i=1,3
      dd(i)=-d(i)
      ds(i)=d(i)**2
      v(i)=0.d0
   1  continue
      g=.5d0/ga
      pre=sqrpi*sqrt(2.d0*g)
      lh=l/2
      mh=m/2
      nh=n/2
      if(2*nh-n) 4,3,4
   3  v(1)=1.d0
      if(l.le.0.and.dd(1).eq.0.d0) go to 29
  28  v(1)=dd(1)**l
   29 if(m.le.0.and.dd(2).eq.0.d0) go to 31
  30  v(1)=v(1)*dd(2)**m
  31  v(1)=pre*g**nh*v(1)*exp(-ga*(ds(1)+ds(2)))*dftr(nh+1)
   4  if(2*mh-m) 6,5,6
   5  v(2)=1.d0
      if(l.le.0.and.dd(1).eq.0.d0) go to 49
  48  v(2)=dd(1)**l
   49 if(n.le.0.and.dd(1).eq.0.d0) go to 51
  50  v(2)=v(2)*dd(1)**n
  51  v(2)=pre*g**mh*v(2)*exp(-ga*(ds(1)+ds(3)))*dftr(mh+1)
   6  if(2*lh-l) 8,7,8
   7  v(3)=1.d0
      if(m.le.0.and.dd(2).eq.0.d0) go to 69
  68  v(3)=dd(2)**m
   69 if(n.le.0.and.dd(3).eq.0.d0) go to 71
  70  v(3)=v(3)*dd(3)**n
  71  v(3)=pre*g**lh*v(3)*exp(-ga*(ds(2)+ds(3)))*dftr(lh+1 )
   8  return
      end
c
c
      subroutine opac3(l,m,n,ga,v,d)
c           electronic part of charge density
      implicit real*8(a-h,o-z)
      dimension v(1),d(1),dd(3)
      ddsq=0.d0
      do 1 i=1,3
      dd(i)=-d(i)
      ddsq=ddsq+d(i)**2
    1 continue
      v(1) =1.d0
      if(l) 2,3,2
    2 v(1)=dd(1)**l
    3 if(m) 4,5,4
    4 v(1)=v(1)*dd(2)**m
    5 if(n) 6,7,6
    6 v(1)=v(1)*dd(3)**n
    7 v(1)=v(1)*exp(-ga*ddsq)
      return
      end
c
c
      subroutine opad1(l,m,n,ga,v,d)
c           electronic part of overlap
      implicit real*8(a-h,o-z)
      dimension v(1),d(1)
      v(1)=olap(l,m,n,ga)
      return
      end
c
c
      subroutine opap4(l,m,n,gc,v,rc)
c
c     this subroutine calculates the expectation-values of
c     -p**2/2 (kinetic energy) and p**4/(8*cl**2) (relativistic
c     correction of the kinetic energy)
c
c     the program works correct only for s and p and the irreducible
c     linear combinations of the primitive d,f,... functions
c
c     written by s. brode in april,1984
c
c
      implicit real*8 (a-h,o-z)
      logical lodd,modd,nodd,leven,meven,neven
      dimension v(3),rc(3),rca(3),rcb(3),srcab(3),
     &lmni(35),o(25),dftr(7)
      data a0/0.0d0/,a1/1.0d0/,a2/2.0d0/,a4/4.0d0/,a8/8.0d0/,
     &lmni/0,3*1,6*2,10*3,15*4/,cl/137.03604d0/,pi/3.1415926559d0/,
     &dftr/1.0d0,1.0d0,3.0d0,15.0d0,105.0d0,945.0d0,103935.0d0/
      common /abfunc/ ra(3),rb(3),ga,gb,ia,ib
c
c     define the overlap-integral
c
      ola(ll,mm,nn)=gch**(ll/2+mm/2+nn/2)*dftr(ll/2+1)*dftr(mm/2+1)*
     &dftr(nn/2+1)
c
c     some previous calculations
c
      do 10 i=1,25
   10 o(i)=a0
      gch=a1/(a2*gc)
      ropigc=sqrt(pi/gc)*pi/gc
      lodd=mod(l,2).eq.1
      leven=.not.lodd
      modd=mod(m,2).eq.1
      meven=.not.modd
      nodd=mod(n,2).eq.1
      neven=.not.nodd
      ga2=ga*a2
      gb2=gb*a2
      gab4=ga2*gb2
      prcaa=a0
      prcbb=a0
      do 20 i=1,3
      rca(i)=rc(i)-ra(i)
      rcb(i)=rc(i)-rb(i)
      srcab(i)=rca(i)+rcb(i)
      prcaa=prcaa+rca(i)*rca(i)
      prcbb=prcbb+rcb(i)*rcb(i)
   20 continue
      zwa=ga2*prcaa-dfloat(2*lmni(ia)+3)
      zwb=gb2*prcbb-dfloat(2*lmni(ib)+3)
c
c     calculation of the overlap-integrals
c
      if(lodd.or.modd.or.nodd) go to 110
      o( 1)=ola(l   ,m  ,n  )
      o( 5)=ola(l+2,m   ,n  )
      o( 6)=ola(l   ,m+2,n  )
      o( 7)=ola(l   ,m  ,n+2)
      o(11)=ola(l+4,m   ,n  )
      o(12)=ola(l   ,m+4,n  )
      o(13)=ola(l   ,m  ,n+4)
      o(20)=ola(l+2,m+2,n   )
      o(21)=ola(l+2,m   ,n+2)
      o(22)=ola(l   ,m+2,n+2)
  110 continue
c
      if(leven.or.modd.or.nodd) go to 120
      o( 2)=ola(l+1,m   ,n  )
      o( 8)=ola(l+3,m   ,n  )
      o(14)=ola(l+1,m+2,n   )
      o(15)=ola(l+1,m   ,n+2)
  120 continue
c
      if(lodd.or.meven.or.nodd) go to 130
      o( 3)=ola(l   ,m+1,n  )
      o( 9)=ola(l   ,m+3,n  )
      o(16)=ola(l+2,m+1,n   )
      o(17)=ola(l   ,m+1,n+2)
  130 continue
c
      if(lodd.or.modd.or.neven) go to 140
      o( 4)=ola(l   ,m  ,n+1)
      o(10)=ola(l   ,m  ,n+3)
      o(18)=ola(l+2,m   ,n+1)
      o(19)=ola(l   ,m+2,n+1)
  140 continue
c
      if(leven.or.meven.or.nodd) go to 150
      o(23)=ola(l+1,m+1,n   )
  150 continue
      if(leven.or.modd.or.neven) go to 160
      o(24)=ola(l+1,m   ,n+1)
  160 continue
      if(lodd.or.meven.or.neven) go to 170
      o(25)=ola(l   ,m+1,n+1)
  170 continue
c
c     calculation of kinetic energy
c
      p2=-ga*(zwa*o(1)
     &       +ga2*(a2*(rca(1)*o(2)+rca(2)*o(3)+rca(3)*o(4))
     &       +o(5)+o(6)+o(7)))
c
c     calculation of relativistic correction
c
      p4=(gab4*(o(11)+o(12)+o(13)+a2*(o(20)+o(21)+o(22))
     &         +a2*(srcab(1)*(o( 8)+o(16)+o(18))
     &             +srcab(2)*(o(14)+o( 9)+o(19))
     &             +srcab(3)*(o(15)+o(17)+o(10)))
     &         +a4*(rca(1)*rcb(1)*o(5)
     &             +rca(2)*rcb(2)*o(6)
     &             +rca(3)*rcb(3)*o(7)
     &             +a2*((rca(1)*rcb(2)+rca(2)+rcb(1))*o(23)
     &                 +(rca(1)*rcb(3)+rca(3)+rcb(1))*o(24)
     &                 +(rca(2)*rcb(3)+rca(3)+rcb(2))*o(25))))
     &   +(gb2*zwa+ga2*zwb)*(o(5)+o(6)+o(7))
     &   +a2*(gb2*zwa*(rcb(1)*o(2)+rcb(2)*o(3)+rcb(3)*o(4))
     &       +ga2*zwb*(rca(1)*o(2)+rca(2)*o(3)+rca(3)*o(4)))
     &   +zwa*zwb*o(1))
     &   *(-gab4)/(a8*cl**2)
c
c     store p2 and p4 to the correct storage-locations
c
      v(1)=p2*ropigc
      v(2)=p4*ropigc
      v(3)=v(1)+v(2)
      return
      end
c
c
      real*8 function aainer(r,s,t,l,m,n,d,dp,fnu,fn,fd)
      implicit real*8(a-h,o-z)
      dimension  d(3),dp(1),fnu(1),fn(1),fd(1)
      integer r,s,t,u,v,w,uvwt,rstt,uvwth
      rstt=r+s+t
      lmnt=l+m+n
      lmnrst=rstt+lmnt
      lh=l/2
      mh=m/2
      nh=n/2
      aainer=0.0d0
      last1 = lh+1
      do 6 ii=1,last1
      i = ii-1
      il=l-2*i+1
      ilr=r+il
      ilrh=(ilr-1)/2
      fi=fn(ilr)*fd(il)*fd(i+1)
      last2 = mh+1
      do 5 jj=1,last2
      j = jj-1
      jm=m-2*j+1
      jms=s+jm
      jmsh=(jms-1)/2
      fij=fn(jms)*fd(jm)*fd(j+1)*fi
      last3 = nh+1
      do 4 kk=1,last3
      k = kk-1
      kn=n-2*k+1
      knt=t+kn
      knth=(knt-1)/2
      ijkt=i+j+k
      fijk=fn(knt)*fd(kn)*fd(k+1)*dp(ijkt+1)*fij
      lmrsij=lmnrst-2*ijkt
      last4 = ilrh+1
      do 3 iu=1,last4
      u = iu-1
      ilru=ilr-2*u
      fu=fd(u+1)*fd(ilru)
      if(abs(d(1))-0.0000001d0)10,10,11
   10 if(ilru-1) 12,12, 3
   11 fu=fu*d(1)**(ilru-1)
   12 last5 = jmsh+1
      do 2 iv=1,last5
      v = iv-1
      jmsv=jms-2*v
      fuv=fu*fd(v+1)*fd(jmsv)
      if(abs(d(2))-0.0000001d0)20,20,21
   20 if(jmsv-1) 22,22,2
   21 fuv=fuv*d(2)**(jmsv-1)
   22 last6 = knth+1
      do 1 iw=1,last6
      w = iw-1
      kntw=knt-2*w
      fuvw=fuv*fd(w+1) *fd(kntw)
      if(abs(d(3))-0.0000001d0)30,30,31
   30 if(kntw-1) 32,32,1
   31 fuvw=fuvw*d(3)**(kntw-1)
   32 uvwt=u+v+w
      uvwth=uvwt/2
      if(2*uvwth-uvwt) 33,34,33
   33 fuvw=-fuvw
   34 nuindx=lmrsij-uvwt
      fuvw=fuvw*fnu(nuindx  +1)*dp(uvwt+1)
      aainer=fijk*fuvw+aainer
    1 continue
    2 continue
    3 continue
    4 continue
    5 continue
    6 continue
      return
      end
c
c
      subroutine divpt(a,alpha,b,beta,ecoord,eexp,efact)
      implicit real*8(a-h,o-z)
c         this computes the centre, exponent, and multiplying factor of
c         a single gaussian which can replace the product of two gaussia
c         centres a and b, and exponents alpha and beta.
      dimension a(3), b(3), ecoord(3)
      do 1 i=1,3
    1 ecoord(i)=(alpha*a(i)+beta*b(i))/(alpha+beta)
      eexp=alpha+beta
      absqd=0.0d0
      do 2 i=1,3
    2 absqd=absqd+(b(i)-a(i))*(b(i)-a(i))
      efact=exp(-alpha*beta*absqd/(alpha+beta))
      return
      end
c
c
      subroutine elcmom(nt,frocc,dint,nro,norb,tot,cint,nomx)
c           total property in mo basis * occupation numbers
      implicit real*8(a-h,o-z)
      dimension frocc(*),dint(nro,*),tot(*),cint(nomx,*)
      do 11 nu=1,nt
      tot(nu)=0.0d0
      do 12 n=1,norb
      tot(nu)=tot(nu)  +  frocc(n)*cint(n,nu)
      dint(2,nu)=tot(nu)
   12 continue
   11 continue
      return
      end
c
c
      subroutine total(nt,dint,nro)
c           nuclear + electronic part of property
      implicit real*8(a-h,o-z)
      dimension dint(nro,*)
      do 20 k=1,nt
   20 dint(3,k)=dint(1,k)+dint(2,k)
      return
      end
c
c
      subroutine output(dint,norb,n,nt,nro,frocc,tot,c,ilbl,a,
     &                  title,cint,buff)
      implicit real*8(a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
      real*8    c(1), ilbl(1), a(1), title(1)
      parameter (mloc=250)
      dimension dint(nro,*),frocc(*),tot(nt),cint(nomx,*),buff(lbuff)
c
c     dimensions for localization have to be:
c     help(ncont,norb,nt), xlv(ncont,norb), op(nnorb,nt), d(norb,norb)
c     eig(norb), mon(norb)
c
      dimension mon(mloc),xlv(2),repv(2)
      logical ok
      common /big/ y(nomx*nomx), help(nomx*mloc*3),
     &  d(mloc*(mloc+1)*3/2),op(mloc*mloc)
      common /eig/    eig(mloc)
      common /prptyp/ m
      common /contrc/ chomo,ncont
      equivalence     (help(1),xlv(1),repv(1))
c
      if(m.eq.15) go to 17
      call cntrac(cint,norb,n,nt,nro,buff)
      call elcmom(nt,frocc,dint,nro,norb,tot,cint,nomx)
      call total(nt,dint,nro)
      call prnt(dint,nro,nt,norb,c,ilbl,a,title,cint,nomx)
      return
c
   17 write(6,6000)
      do 10 i=1,norb
   10 mon(i)=i
c
      read(5,1000) nmo
      if(nmo) 300,200,100
  100 read(5,1000) (mon(i),i=1,nmo)
c
      ok=.true.
      ialt=mon(1)
      do 220 ip=2,nmo
      if(mon(ip).gt.ialt) go to 210
      ok=.false.
      write(6,2222) ip,mon(ip),ialt
  210 ialt=mon(ip)
  220 continue
      if (ok) go to 250
      stop
c
  200 read(5,2000) thr
      write(6,6002) thr
      ii=0
      do 20 i=1,norb
      if(eig(i).lt.thr) go to 20
      ii=ii+1
      mon(ii)=i
   20 continue
      nmo=ii
c
  250 ij=0
      do 30 i=1,nmo
      eig(i)=eig(mon(i))
      ist=(mon(i)-1)*ncont+1
      ie=ist+ncont-1
      do 40 k=ist,ie
      ij=ij+1
   40 y(ij)=y(k)
   30 continue
      norb=nmo
c
  300 write(6,6001) norb,(mon(i),i=1,norb)
      nnorb=norb*(norb+1)/2
      call local(nt,norb,ncont,op,nnorb,buff,lbuff,y,help,d,xlv,repv)
      return
c
 1000 format(5x,5i5)
 2000 format(5x,f10.4)
 2222 format('  ====> danger : ordering of mos illegal ',i5,5x,2i5)
 6000 format(//' boys localization '//)
 6001 format(1x,i5,' occupied mo s selected :',5i5/(31x,5i5))
 6002 format(6x,' threshold for selection of mo s :',f10.4/)
      end
c
c
      subroutine local(nt,norb,ncont,op,nnorb,buff,lbuff,y,help,d,xlv,
     1 repv)
      implicit real*8(a-h,o-z)
      integer opn(3)
      logical ok
      parameter (mloc=250)
      dimension y(ncont,norb),help(ncont,norb,nt),buff(lbuff),
     &          op(nnorb,nt),d(norb,norb),xlv(ncont,norb),repv(2)
      common/holl/blabel(19),alabel(19),fmt(4),jl,lomao,iopen,nsym
      common/supp/nsupp
      common/eig/eig(mloc)
c
c     help(i,j), xlv(i,j) and repv(i) are equivalenced in subroutine
c     output
c
      rewind 2
      do 10 m=1,nt
      do 10 i=1,norb
      do 10 k=1,ncont
   10 help(k,i,m)=0.0d0
      do 20 j=1,nt
   20 opn(j)=j
c
      ix=lbuff
      do 30 k=1,ncont
      do 30 l=1,k
      fac=1.0d0
      if(k.eq.l) fac=0.5d0
      do 50 m=1,nt
      ix=ix+1
      if(ix.le.lbuff)go to 40
      ix=1
      read(2) buff
   40 continue
      temp=-buff(ix)*fac
      do 60 i=1,norb
      help(k,i,m)=help(k,i,m)+temp*y(l,i)
   60 help(l,i,m)=help(l,i,m)+temp*y(k,i)
   50 continue
   30 continue
c
      do 70 m=1,nt
      ij=0
      do 80 i=1,norb
      do 80 j=1,i
      ij=ij+1
      top=0.0d0
      do 90 k=1,ncont
   90 top=top+y(k,i)*help(k,j,m)
      op(ij,m)=top
   80 continue
   70 continue
c
c
      read(5,5000) nop
      if(nop.eq.0) go to 100
      read(5,5000) (opn(nn),nn=1,nop)
c
      if(nop.eq.1) go to 130
      ok=.true.
      ialt=opn(1)
      do 110 ip=2,nop
      if(opn(ip).gt.ialt) go to 120
      ok=.false.
      write(6,6000) ip,opn(ip),ialt
  120 ialt=opn(ip)
  110 continue
      if (ok) go to 130
      stop
c
  130 if(nsupp.lt.2) go to 140
      write(6,6010)
      do 150 i=1,nnorb
  150 write(6,6020) (op(i,kk),kk=1,nt)
c
  140 ii=0
      do 160 n=1,nop
      nopn=opn(n)
      do 160 l=1,nnorb
  160 op(l,n)=op(l,nopn)
      nt=nop
c
  100 write(6,6030) nt,(opn(nn),nn=1,nt)
      thi=1.0d-13
      call boys(norb,nt,nnorb,op,thi,d)
c
      ij=0
      do 300 i=1,norb
      do 300 j=1,i
      ij=ij+1
      repv(ij)=0.d0
      do 300 k=1,norb
  300 repv(ij)=repv(ij)+d(k,i)*eig(k)*d(k,j)
      write(6,6140)
      call outpak(repv,norb,1,2)
c
      do 170 i=1,norb
      do 180 k=1,ncont
      xlv(k,i)=y(k,1)*d(1,i)
      if(norb.eq.1) go to 180
      do 190 j=2,norb
  190 xlv(k,i)=xlv(k,i)+y(k,j)*d(j,i)
  180 continue
  170 continue
      nxlv=norb*ncont
c
      if(nsupp.lt.2) go to 200
      write(6,6040)
      do 210 i=1,nnorb
  210 write(6,6020) (op(i,kk),kk=1,nt)
c
  200 write(6,6050)
      write(6,6060) (opn(i),i=1,nt)
      ii=0
      do 220 i=1,norb
      ii=ii+i
  220 write(6,6070) i,(op(ii,kk),kk=1,nt)
c
      write(6,6080)
      do 230 i=1,norb
  230 write(6,6090) i,(d(j,i),j=1,norb)
c
      if(nsupp.lt.2) go to 240
      write(6,6100)
      do 250 i=1,norb
  250 write(6,6110) i,(xlv(ll,i),ll=1,ncont)
c
  240 write(lomao,6120) (alabel(i1),i1=1,19)
      write(lomao,6120) (blabel(i1),i1=1,19)
      write(lomao,6130) norb,ncont,nxlv
      write(lomao,6120) (fmt(i1),i1=1,4)
      write(lomao,fmt) ((xlv(ll,i),ll=1,ncont),i=1,norb)
c
      return
 5000 format(5x,5i5)
 6000 format(' ====> danger : ordering of operators illegal ',i5,5x,2i5)
 6010 format(/' operators in mo basis '/)
 6020 format(4f20.14)
 6030 format(/1x,i5,' operators selected     :',10i5)
 6040 format(/' operators in lmo basis '/)
 6050 format(//' diagonalelements of the operators in lmo basis :'/)
 6060 format(8x,'lmo',12x,3(i2,12x)/)
 6070 format(5x,i5,5x,3f14.8)
 6080 format(/' lmo s as linearcombination of mo s ')
 6090 format(/' lmo nb.',i5/(1x,5f14.8))
 6100 format(/' lmo vectors '/)
 6110 format(' lmo nb.',i5/(1x,4f14.8))
 6120 format(19a4)
 6130 format(20i4)
 6140 format(/' fockmatrix in lmo basis',/)
      end
c
c
      subroutine boys(n,no,nn1,op,thi,d)
c
c     subroutine for boys localization  r.ahlrichs 2/85
c     to determine vectors |i>, i=1,n
c     which maximize
c     f = sum(i,j=1,n)sum(k=1,no) (<i|op(k)|i>-<j|op(k)|j>)**2
c     equivalent to maximization of
c     g = sum(i=1,n)sum(k=1,no) <i|op(k)|i>**2 + <j|op(k)|j>**2
c     equivalent to minimization of
c     h = sum(i<j)sum(k) <i|op(k)|j>**2
c     stationarity condition
c     sum(k) op(ij,k)*(op(ii,k)-op(jj,k)) = 0 all i.ne.j
c
      implicit real*8(a-h,o-z)
      dimension d(n,n),op(nn1,no)
c
c     input
c           n = dimension of lin. space
c           no = # of operators in f or g
c           op(ij,k) matrix representation of k'th operator
c                symmetric and stored triangular,nn1=n(n+1)/2
c           thi = rotation threshold (1.-6 to 1.-12 sngl or dbl prc)
c     result
c            d contains the solution vectors as columns
c            op the transformed operators
c
      a0=0.d0
      a1=1.d0
      a2=2.d0
      a10=10.d0
      aqu=0.25d0
      accr=1.d-15
      th=dmax1(thi,accr)
c
c     set d to unit matrix
c
      do 20 i=1,n
      do 10 j=1,n
  10  d(j,i)=a0
  20  d(i,i)=a1
      if(n.eq.1) return
c
c     opn = operator norm needed for cut off threshold
c
      opn=a0
      do 22 k=1,no
      do 21 i=1,nn1
  21  opn=opn+op(i,k)**2
  22  continue
      if(opn.eq.a0) return
      opn=sqrt(opn/dfloat(nn1*no))*accr
      ths=aqu
      it=0
      do 200 isw=1,9999
c
c     loop over sweeps
c
      thsw=a0
      ii=1
      ij=2
      do 110 i=2,n
      ii=ii+i
      jj=0
      im1=i-1
      ip1=i+1
      li=ii-i
      do 100 j=1,im1
c
c     loop over i,j rotations
c     maximize sum(k) <i|op(k)|i>**2 + <j|op(k)|j>**2
c     like in jacobi
c
      jj=jj+j
c
c     get ax , b determining the tan of the rotation angle
c
      ax=a0
      b=a0
      do 30 k=1,no
      aux1=op(jj,k)-op(ii,k)
      aux2=op(ij,k)
      ax=ax+aux1*aux2
  30  b=b+aqu*aux1**2-aux2**2
c
c     rotation angle x0
c
      if (abs(ax)+abs(b).lt.opn) go to 100
      x0=aqu*datan2(ax,b)
      thsw=dmax1(thsw,abs(x0))
      if (abs(x0).lt.ths) go to 100
c
c     rotation
c
      it=it+1
c     c=dcos(x0)
      s=dsin(x0)
      s2=s*s
      c2=a1-s2
      c=sqrt(c2)
      do 40 l=1,n
      aux1=d(l,j)
      d(l,j)=c*d(l,j)+s*d(l,i)
  40  d(l,i)=c*d(l,i)-s*aux1
      jm1=j-1
c     c2=c*c
c     s2=s*s
      cs=c*s
      c2ms2=c2-s2
      lj=jj-j
      jp1=j+1
      do 95 k=1,no
      aux1=op(ii,k)
      aux2=op(jj,k)
      op(ii,k)=aux1*c2+aux2*s2-a2*op(ij,k)*cs
      op(jj,k)=aux1*s2+aux2*c2+a2*op(ij,k)*cs
      op(ij,k)=cs*(aux1-aux2)+op(ij,k)*c2ms2
      if(jm1.le.0) go to 60
      do 50 l=1,jm1
      aux1=c*op(li+l,k)-s*op(lj+l,k)
      op(lj+l,k)=s*op(li+l,k)+c*op(lj+l,k)
  50  op(li+l,k)=aux1
  60  if (im1.eq.j) go to 80
      jl=jj+j
      do 70 l=jp1,im1
      aux1=c*op(li+l,k)-s*op(jl,k)
      op(jl,k)=c*op(jl,k)+s*op(li+l,k)
      op(li+l,k)=aux1
  70  jl=jl+l
  80  if (i.eq.n) go to 95
      jl=ii+j
      il=ii+i
      do 90 l=ip1,n
      aux1=c*op(jl,k)+s*op(il,k)
      op(il,k)=c*op(il,k)-s*op(jl,k)
      op(jl,k)=aux1
      il=il+l
  90  jl=jl+l
  95  continue
 100  ij=ij+1
 110  ij=ij+1
c
c     check convergence
c
      if (thsw.le.th) go to 210
      ths=dmin1(thsw**2,ths*aqu)
      ths=dmax1(th,ths)
c
c     write (6,199) isw,it,thsw,ths
c199  format (' isw,it,thsw,,ths',2i5,2d16.8)
c
 200  continue
      write (6,300) isw,thsw
      return
 210  write (6,310) it,isw,thsw
      return
 300  format (///' boys not converged in',i6,' sweeps thsw =',
     *d16.8/)
 310  format (/' boys converged after',i6,' pivots, in ',i6
     *,' sweeps, thsw = ',d16.8/)
      end
c
c
      subroutine mmult(a,b,bs,n)
      implicit real*8(a-h,o-z)
      dimension a(n,1), b(n,1), bs(n,1)
      do 6 i=1,n
      do 6 l=1,n
      bs(i,l) = 0.0d0
      do 3 j=1,n
      do 3 k=1,n
      bs(i,l) = bs(i,l) + a(j,i)*b(j,k)*a(k,l)
    3 continue
    6 continue
      do 9 i=1,n
      do 9 j=1,n
    9 b(i,j) = bs(i,j)
      return
      end
c
c
      subroutine prnt(dint,nro,nt,norb,c,ilbl,a,title,cint,nomx)
c
c     this subroutine transforms into the principal axis system and
c     produces the final output
c
      implicit real*8(a-h,o-z)
      dimension dint(nro,*),c(3),a(13,3),title(15,10),cint(nomx,*)
      real*8 ilbl(12)
      dimension pr1(3,3), pr2(3,3), u(3,3)
      dimension name(3)
      data name/2hx*,2hy*,2hz*/
      data audeb/2.541770d0/
      common/ioind/isl1, isl2
      common/prptyp/mm
      common/ncnt/ilnmcm
      common/supp/nsupp
c       file shout short output, nshout=13
      nshout=13
      mmm=0
      if(mm.eq.13) mmm=2
c     the titels for the relat. correction are stored at title(15,j)
c     not at title(13,j) !
      write(6,108) (ilbl(i), i=1,12)
      if(nsupp.lt.1) goto 50
      write(6,104) (a(mm,i), i=1,3)
      write(6,201)
      write(6,300) (title(mm+mmm,i), i=1,nt)
      write(6,301)
      do 1 i=1,norb
      write(6,106)i,i,(cint(i,k),k=1,nt)
    1 continue
   50 write(6,108) (ilbl(i), i=1,12)
      write(6,200) (a(mm,i), i=1,3), ilnmcm, (c(j), j=1,3)
      write(nshout,1200) (a(mm,i), i=1,3), ilnmcm, (c(j), j=1,3)
      if(mm.eq.13) write(6,109)
      write(6,201)
      write(nshout,1201)
      write(6,300) (title(mm+mmm,i), i=1,nt)
      write(nshout,1300) (title(mm+mmm,i), i=1,nt)
      write(6,301)
      if(mm.ne.13) go to 6000
      write(6,6101) (dint(1,l), l=1,nt)
      write(6,6102) (dint(2,m), m=1,nt)
      write(6,6103) (dint(3,n), n=1,nt)
      go to 6001
6000  write(6,101) (dint(1,l), l=1,nt)
      write(6,102) (dint(2,m), m=1,nt)
      write(6,103) (dint(3,n), n=1,nt)
      write(nshout,1101) (dint(1,l), l=1,nt)
      write(nshout,1102) (dint(2,m), m=1,nt)
      write(nshout,1103) (dint(3,n), n=1,nt)
c     output of electrostatic potential on file 21
      if(mm.eq.1) write(21,20000) c(1),c(2),c(3),dint(3,1)
20000 format(4f10.6)
6001  continue
      if(mm.ne.4) goto 61
      sum=0.d0
      do 60 i=1,3
   60 sum=sum+dint(3,i)**2
      sum=sqrt(sum)*audeb
      write(6,403) sum
      write(nshout,1403) sum
c     only the following properties are transformed into the principal
c     axis system: field gradient, quadrupole moment, second moment,
c     diamagnetic shielding
   61 if(mm.eq.3) go to 12
      if(mm.eq.5) go to 12
      if(mm.eq.11) go to 12
      if(mm.eq.12)go to 12
      return
   12 flag =abs(dint(3,4)) +abs(dint(3,5)) +abs(dint(3,6))
      if(flag.lt.3.0d-05) return
      if(mm.eq.3) ll=12
      if(mm.eq.5) ll=13
      if(mm.eq.11) ll=14
      if(mm.eq.12) ll=12
      if(nsupp.lt.1) goto 51
      write(6,108) (ilbl(i), i=1,12)
      write(6,104) (a(mm,i), i=1,3)
      write(6,400)
      write(6,300) (title(ll,i), i=1,nt)
      write(6,1400)
      write(6,1300) (title(ll,i), i=1,nt)
      write(6,301)
   51 pr1(1,1) = dint(3,1)
      pr1(1,2) = dint(3,4)
      pr1(1,3) = dint(3,5)
      pr1(2,2) = dint(3,2)
      pr1(2,3) = dint(3,6)
      pr1(3,3) = dint(3,3)
      call hdiag(pr1,3,0,u,nr)
c     transform electronic, total, and nuclear parts.
      do 14 i=1,3
      pr2(1,1)=dint(i,1)
      pr2(1,2)=dint(i,4)
      pr2(1,3)=dint(i,5)
      pr2(2,2)=dint(i,2)
      pr2(2,3)=dint(i,6)
      pr2(3,3)=dint(i,3)
      do 18 ii=1,3
      do 18 ij=ii,3
   18 pr2(ij,ii)=pr2(ii,ij)
      call mmult(u,pr2,pr1,3)
      dint(i,1)=pr2(1,1)
      dint(i,4)=pr2(1,2)
      dint(i,5)=pr2(1,3)
      dint(i,2)=pr2(2,2)
      dint(i,6)=pr2(2,3)
      dint(i,3)=pr2(3,3)
   14 continue
c     transform mo integrals.
      do 22 i=1,norb
      pr2(1,1)=cint(i,1)
      pr2(1,2)=cint(i,4)
      pr2(1,3)=cint(i,5)
      pr2(2,2)=cint(i,2)
      pr2(2,3)=cint(i,6)
      pr2(3,3)=cint(i,3)
      do 15 ii=1,3
      do 15 ij=ii,3
   15 pr2(ij,ii) = pr2(ii,ij)
      call mmult(u,pr2,pr1,3)
      cint(i,1)=pr2(1,1)
      cint(i,4)=pr2(1,2)
      cint(i,5)=pr2(1,3)
      cint(i,2)=pr2(2,2)
      cint(i,6)=pr2(2,3)
      cint(i,3)=pr2(3,3)
      if(nsupp.lt.1) goto 22
      write(6,106)i,i,(cint(i,k),k=1,nt)
   22 continue
      write(6,108) (ilbl(i), i=1,12)
      write(6,200) (a(mm,i), i=1,3), ilnmcm, (c(k), k=1,3)
      write(6,400)
      write(6,304) (title(ll,i), i=1,nt)
      write(6,101) (dint(1,l), l=1,6)
      write(6,102) (dint(2,l), l=1,6)
      write(6,103) (dint(3,l), l=1,6)
   37 write(6,303)
      write(6,401)
      do 24 i=1,3
   24 write(6,202) name(i), (u(j,i), j=1,3)
      return
  200 format(   //,1x, 3a6,/,1x,38hcalculated with reference to the poin
     1t,5x,a6,/,1x,11hcoordinates,/, 6x,3hx =,f12.6,5x,3hy =,f12.6,5x,
     23hz =,f12.6)
 1200 format(1x, 3a6,/,1x,38hcalculated with reference to the point,
     15x,a6,/,1x,11hcoordinates,/, 6x,3hx =,f12.6,5x,3hy =,f12.6,5x,
     23hz =,f12.6)
  201 format(/,1x,42hin atomic units and for global coordinates)
 1201 format(1x,42hin atomic units and for global coordinates)
  202 format(1x,a2,1x,3d15.6)
  300 format(//,10x,7(9x,a6),/,10x,3(9x,a6))
 1300 format(10x,7(9x,a6),/,10x,3(9x,a6))
  301 format(/)
  303 format(///)
  304 format(//,12x,7(9x,a6),/,12x,3(9x,a))
  400 format(/,1x,49hin atomic units and for the principle axis system)
 1400 format(1x,49hin atomic units and for the principle axis system)
  401 format(1x,67hthe new coordinate system, (x*,y*,z*), is related to
     1the old one by,//,16x,1hx,14x,1hy,14x,1hz,//)
  403 format(1x,'norm ',8x,d15.7,'  debeye')
 1403 format(1x,'norm ',8x,d15.7,'  debeye')
  101 format(/,1x,7hnuclear,5x,7(1x,d15.7),/,13x,3(1x,d15.7))
  102 format(1x,10helectronic,2x,7(1x,d15.7),/,13x,3(1x,d15.7))
  103 format(/,1x,5htotal,7x,7(1x,d15.7),/,13x,3(1x,d15.7))
 1101 format(1x,7hnuclear,5x,7(1x,d15.7),/,13x,3(1x,d15.7))
 1102 format(1x,10helectronic,2x,7(1x,d15.7),/,13x,3(1x,d15.7))
 1103 format(1x,5htotal,7x,7(1x,d15.7),/,13x,3(1x,d15.7))
  104 format(   //,1x,23hintegrals over mo s for,5x,3a6)
  106 format(1x,2i5,7(1x,d15.7),/,11x,3(1x,d15.7))
  108 format(1h0,12a6)
  109 format(/,1x,'the sum of the nuclear part means the darwin-term.',
     &' the relativisticly corrected '/' kinetic energy is printed as ',
     &'the total sum.')
 6101 format(/,1x,7hdarwin ,5x,3(1x,f25.6))
 6102 format(1x,10helectronic,2x,3(1x,f25.6))
 6103 format(/,1x,5htotal,7x,3(1x,f25.6))
      end
c
c
      subroutine hdiag(h,n,iegen,u,nr)
c     diagonalization of a real symmetric matrix by the jacobi method
c            where h is the array to be diagonalized.
c     n is the order of the matrix, h.
c     returns eigval in x
c     iegen must be set unequal to zero if only eigenvalues are
c            to be computed.
c     iegen must be set equal to zero if eigenvalues and eigenvectors
c            are to be computed.
c     u is the unitary matrix used for formation of the eigenvectors.
c     nr is the number of rotations.
c     the subroutine operates only on the elements of h that are to the
c            right of the main diagonal.
      implicit real*8(a-h,o-z)
      dimension h(3,3), u(3,3), x(3), iq(3)
    9 if (iegen) 15,10,15
 10   do 14 i=1,n
      do 14 j=i,n
      if(i-j)12,11,12
 11   u(i,j)=1.0d0
      go to 14
 12   u(i,j)=0.d0
      u(j,i) = 0.d0
 14   continue
 15   nr = 0
      if (n-1) 2001,2001,17
c     scan for largest off diagonal element in each row
c     x(i) contains largest element in ith row
c     iq(i) holds second subscript defining position of element
 17   nmi1=n-1
      hdtest = 0.0d0
      do 30 i=1,nmi1
      x(i) = 0.d0
      ipl1=i+1
      do 30 j=ipl1,n
      hdtest = hdtest + h(i,j)**2
      if (x(i).gt.abs(h(i,j))) go to 30
   20 x(i) =abs(h(i,j))
      iq(i)=j
 30   continue
      hdtest = 2.0d0*hdtest
c     set indicator for shut-off.rap=2**-27,nr=no. of rotations
      rap=1.d-12
      do 35 i=1,n
   35 hdtest = hdtest + h(i,i)**2
      hdtest =rap*(sqrt(hdtest/dfloat(n)) + 1.0d0)
c     rap*(rms eigval + 1)
c     find maximum of x(i) s for pivot element and
c     test for end of problem
 40   do  70  i=1,nmi1
      if (i-1) 60,60,45
 45   if ( xmax- x(i)) 60,70,70
 60   xmax=x(i)
      ipiv=i
      jpiv=iq(i)
 70   continue
c     return if max.h(i,j)less than(2**-27)absf(h(k,k)-min)
      if (hdtest- xmax) 148,1000,1000
 148  nr = nr+1
c     compute tangent, sine and cosine,h(i,i),h(j,j)
  150 hhdifo=0.5d0*(h(ipiv,ipiv)-h(jpiv,jpiv))
      hhdifn=dsign(1.0d0,hhdifo)*sqrt(hhdifo**2+h(ipiv,jpiv)**2)
      cosine=sqrt(0.5d0*(hhdifo+hhdifn)/hhdifn)
      sine=h(ipiv,jpiv)/hhdifn*0.5d0/cosine
      hhsum=0.5d0*(h(ipiv,ipiv)+h(jpiv,jpiv))
      h(ipiv,ipiv)=hhsum+hhdifn
      h(jpiv,jpiv)=hhsum-hhdifn
      h(ipiv,jpiv)=0.0d0
c     inspect the iqs between i+1 and n-1 to determine
c     whether a new maximum value should be computed since
c     the present maximum is in the i or j row.
      do 350 i=1,nmi1
      if(i-ipiv)210,350,200
 200  if(i-jpiv)210,350,210
 210  if(iq(i)-ipiv)230,240,230
 230  if(iq(i)-jpiv)350,240,350
  240 htempi = h(i,ipiv)
      htempj = h(i,jpiv)
      h(i,ipiv) = 0.0d0
      h(i,jpiv) = 0.0d0
      ipl1=i+1
      x(i) =0.d0
c     search in depleted row for new maximum
      do 320 j=ipl1,n
      if (x(i).ge.abs(h(i,j))) go to 320
  300 x(i) =abs(h(i,j))
      iq(i)=j
 320  continue
      h(i,ipiv)=htempi
      h(i,jpiv)=htempj
 350  continue
      x(jpiv) =0.d0
      x(ipiv) =0.d0
c     change the other elements of h
      do 530 i=1,n
      if(i-ipiv)370,530,420
 370  htemp = h(i,ipiv)
      h(i,ipiv) = cosine*htemp + sine*h(i,jpiv)
      if (x(i).ge.abs(h(i,ipiv))) go to 390
      x(i) =abs(h(i,ipiv))
      iq(i) = ipiv
 390  h(i,jpiv) = -sine*htemp + cosine*h(i,jpiv)
      if (x(i).ge.abs(h(i,jpiv))) go to 530
  400 x(i) =abs(h(i,jpiv))
      iq(i) = jpiv
      go to 530
 420  if(i-jpiv)430,530,480
 430  htemp = h(ipiv,i)
      h(ipiv,i) = cosine*htemp + sine*h(i,jpiv)
      if (x(ipiv).ge.abs(h(ipiv,i))) go to 450
  440 x(ipiv) =abs(h(ipiv,i))
      iq(ipiv) = i
 450  h(i,jpiv) = -sine*htemp + cosine*h(i,jpiv)
      if (x(i)-abs(h(i,jpiv))) 400,530,530
 480  htemp = h(ipiv,i)
      h(ipiv,i) = cosine*htemp + sine*h(jpiv,i)
      if (x(ipiv).ge.abs(h(ipiv,i))) go to 500
      x(ipiv) =abs(h(ipiv,i))
      iq(ipiv) = i
 500  h(jpiv,i) = -sine*htemp + cosine*h(jpiv,i)
      if (x(jpiv).ge.abs(h(jpiv,i))) go to 530
      x(jpiv) =abs(h(jpiv,i))
      iq(jpiv) = i
 530  continue
c     test for computation of eigenvectors
      if(iegen)40,540,40
 540  do 550 i=1,n
      htemp=u(i,ipiv)
      u(i,ipiv)=cosine*htemp+sine*u(i,jpiv)
 550  u(i,jpiv)=-sine*htemp+cosine*u(i,jpiv)
      go to 40
 1000 continue
      call schort(u,n)
 2000 return
 2001                  x(1) = h(1,1)
      go to 2000
      end
c
c
      subroutine schort(u,n)
      implicit real*8(a-h,o-z)
      dimension u(3,3)
      do 605 i=1,n
      if (i.eq.1) go to 601
      im1 = i-1
      do 603 j=1,im1
      term = 0.0d0
      do 602 k=1,n
  602 term = term + u(k,i)*u(k,j)
      do 603 k=1,n
  603 u(k,i) = u(k,i) - term*u(k,j)
  601 term = 0.0d0
      do 604 k=1,n
  604 term = term + u(k,i)**2
      term = 1.0d0/sqrt(term)
      do 605 k=1,n
  605 u(k,i) = u(k,i)*term
      return
      end
      subroutine inpsym(mcu,msu,mru,mctu,mcsu,icomx,mgcsu)
c
c     this subroutine reads the data from the karlsruhe integral input
c     program (file ft24)
c
      implicit real*8(a-h,o-z)
      parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
     &          lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
     &          lpru=18000,lsu=50,lstu=8,lconu=12,
     &          lcu=8,lccu=182)
      common/supp/nsupp
      common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
     1mcxu,ngcs,mdum(45),mccu,mblu,ns,ng,nsf,nst,ncons
      common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
     1  ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
     2    eta(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
     3  z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
     4  la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
     5  mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu),iic(lctu,lcsu)
      dimension blabel(19)
      dimension ipq(256),lfmt(2)
      data a0,a1/0.0d0,1.0d0/
      itape=24
c        file aoinp input for ao basis
c        itape=24
      open(unit=itape,file='argosin',status='old',
     1access='sequential',form='formatted')
      ntape=0
      inp=itape
      inp2=31
c        file geomfl
c        inp2=31
**      open(unit=inp2,file='geomfl',status='unknown',
c**     1access='sequential',form='formatted')
c**      rewind inp2
      inpn=25
c        file inpn is a copy of aoinp where the general
c        contractions are resolved into segmented ones
c**      open(unit=inpn,file='aoinpn',status='unknown',
c**     1access='sequential',form='formatted')
      call inpnew(mcu,msu,mru,mctu,mcsu,icomx,mgcsu,inp,inp2,inpn)
      itape=inpn
      inp=itape
      rewind itape
      ipq(1)=0
      do 1000 i=1,255
 1000 ipq(i+1)=ipq(i)+i
      read(itape,890)blabel
c**      read(itape,8900) (lfmt(i),i=1,2)
 8900 format(19a4)
 7900 format(1x,19a4)
      read(itape,*)ngen,ns,naords,ncons,ngcs,itol,icut,idummy,idummy,
c    & icp, nkoord
     1 itgeom
c**      if(itgeom.gt.0) read(inp2,999) niter
c**      if(itgeom.gt.0) inp=inp2
      if(ns.le.msu) goto 2
      write(6,970) ns
      stop
    2 if(naords.le.kaords) goto 4
      write(6,971) naords
      stop
    4 if(ncons.le.mconsu) goto 6
      write(6,972) ncons
      stop
    6 if(ngcs.le.mgcsu) goto 8
      write(6,973) ngcs
      stop
    8 if(nsupp.eq.2) write (6,790) blabel
c     if (nkoord .eq. 1) write(6,8770)
c8770 format(/'atomic coordinates were read from file 28'/)
c**      if(nsupp.eq.2) write (6,7900) (lfmt(i),i=1,2)
      if(nsupp.eq.2) write (6,800) ngen, ns, naords, ncons, ngcs, itol,
     1 icut, idummy, nsupp, itgeom
      read(itape,901)nst,(nd(ist),ityp(ist),ist=1,nst)
      if(nsupp.eq.2) write (6,801) nst, (nd(ist), ityp(ist), ist=1,nst)
      if(nst.le.mstu) goto 9
      write(6,974) nst
      stop
    9 read(itape,*) idmx
      if(idmx.eq.0) goto 11
      do 1346 ixj=1,idmx
 1346 read(itape,900)
   11 do 10 ist=1,nst
   10 nso(ist)=0
      do 20 iaords=1,naords
      read(itape,*)iru,(la(ir,iaords),ir=1,iru)
      if(iru.le.mru) goto 20
      write(6,975) iru
      stop
   20 nir(iaords)=iru
      do 30 igcs=1,ngcs
      read(itape,*)icsu,ictu,maords(igcs)
      if(nsupp.eq.2) write (6,800) icsu, ictu, maords(igcs)
      if(icsu.le.mcsu) goto 26
      write(6,976) icsu
      stop
   26 if(ictu.le.mctu) goto 28
      write(6,977) ictu
      stop
   28 do 30 ics=1,icsu
      read(itape,*)(iic(ict,ics),ict=1,ictu)
      if(nsupp.eq.2) write (6,900) (iic(ict,ics), ict=1,ictu)
      do  30 ict=1,ictu
      c1=dfloat(iic(ict,ics))
      c(ict,ics,igcs)=c1
      ccc=abs(c1)
   30 if(ccc.ne.a0.and.ccc.ne.a1)c(ict,ics,igcs)=dsign(sqrt(ccc),c1)
c        read in exponents and contraction coefficients
      do 40 icons=1,ncons
      read (itape,'(i4,i3)') iconu,lmnp1(icons)
      read(itape,*)(zet(icon,icons), eta(icon,icons), icon=1,iconu)
      if(iconu.le.icomx) goto 55
      write(6,2000) icons,icomx
 2000 format(1x,'too many primitive functions in contraction nr. ',i5,
     1 '    icomx = ',i5)
       stop
   55 continue
      do 39 icon=1,iconu
   39 if(eta(icon,icons).eq.0.d0.and.iconu.eq.1)
     1eta(icon,icons)=1.d0
      if(nsupp.eq.2) write (6,803) iconu,lmnp1(icons),(zet(icon,icons),
     1   eta(icon,icons), icon=1,iconu)
   40 ncon(icons)=iconu
      ngenp1=ngen+1
      ng=ngenp1
c        read in atomic labels, charges, coordinates, and basis sets
      isf=0
      isfr=0
      do 75 is=1,ns
      read(inp,904)mtype(is),nf(is),nc(is),chg(is)
      if(nsupp.eq.2) write (6,804) mtype(is), nf(is), nc(is), chg(is)
      icu=nc(is)
      if(icu.le.mcu) goto 46
      write(6,979) icu
      stop
   46 read(inp,905)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
c     if (nkoord .eq. 1)
c    &               read(28,*)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
      if(nsupp.eq.2)write(6,805)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
      if(icu.eq.1) go to 60
      do 50 ig=2,ng
      read(itape,*)(ica(ic,is,ig),ic=1,icu)
   50 if(nsupp.eq.2) write (6,800) (ica(ic,is,ig), ic=1,icu)
   60 ifu=nf(is)
      if(ifu.eq.0) go to 75
      do 70 if=1,ifu
      isf=isf+1
      read(itape,*)mcons(isf),igcs
      if(nsupp.eq.2) write (6,800) mcons(isf), igcs
      iaords=maords(igcs)
      iru= nir(iaords)
      mgcs(isf)=igcs
      ilmnp1=lmnp1(mcons(isf))
      nnt(isf)=ipq(ilmnp1+1)
      ntu(isf)=(nnt(isf)*(ilmnp1+2))/3
      ntl(isf)=ntu(isf)-nnt(isf)+1
      nct(isf)=icu*nnt(isf)
      do 70 ir=1,iru
      isfr=isfr+1
      lai=la(ir,iaords)
      nso(lai)=nso(lai)+1
   70 continue
   75 continue
  790 format(1x,19a4)
  970 format(14h msu too small,i3)
  971 format(17h kaords too small,i3)
  972 format(17h mconsu too small,i4)
  973 format(16h mgcsu too small,i4)
  974 format(15h mstu too small,i3)
  975 format(14h mru too small,i3)
  976 format(15h mcsu too small,i3)
  977 format(15h mctu too small,i3)
  978 format(16h mconu too small,i3)
  979 format(14h mcu too small,i3)
  800 format(1x,26i3)
  801 format(i4,12(i3,a3))
  802 format(1x,26f3.0)
  803 format(i4,i3/(1x,2g14.8))
  804 format(1x,a5,2i3,f3.0)
  805 format(1x,3f14.8)
  890 format(19a4)
  900 format(26i3)
  901 format(i3,12(i3,a3))
  902 format(26f3.0)
  903 format(2i3/(2g14.8))
  904 format(a3,2i3,f3.0)
  905 format(3f14.8)
  949 format(6(/),                                  /////68h0primitive a
     1o integrals neglected if exponential factor below 10**(-,i2,1h)/63
     2h contracted ao and so integrals neglected if value below 10**(-,i
     32,1h)/43h symmetry orbital integrals written on unit,i3)
  951 format(11h degeneracy,16x,12(i4,4x))
  952 format(6h label,21x,12(2x,a4,2x))
  953 format(1h1,12x,a4,6h atoms//15h0nuclear charge,f7.2//7h0center,12x
     1,1hx,15x,1hy,15x,1hz/(i4,7x,3f16.8))
  954 format(/9h0operator,6x,19hcenter interchanges/2x,6h1(id.),6x,24i3)
  955 format(i3,6h(gen.),5x,24i3)
  956 format(i3,11x,24i3)
  958 format(//10x,a2,9h orbitals/18h0orbital exponents,7x,7(1pg15.7))
  959 format(25h contraction coefficients,7(1pg15.7))
  960 format(10h0sym. orb.,12x,15hcenter, ao type/14x,(18(i3,1h,,a2)))
  961 format(i3,a4,i1,6x,(18f6.2))
  999 format(20i4)
      return
      end
      subroutine inpnew(mcu, msu, mru, mctu, mcsu, icomx, mgcsu,
     1  inp,inp2,inpn)
c
c     this subroutine reads the data from the karlsruhe integral input
c     program (file ft24) and makes a segmented contraction copy on
c     file inpn=25
c
c     irun=0   first scan through ao-input-file
c     irun=1   create new ao-input-file (inpn)
c
      implicit real*8(a-h,o-z)
      parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
     &          lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
     &          lpru=18000,lsu=50,lstu=8,lconu=12,
     &          lcu=8,lccu=182,mrcru=8,ione=1,izero=0)
      common/supp/nsupp
      common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
     1mcxu,ngcs,mdum(45),mccu,mblu,ns,ng,nsf,nst,ncons
      common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
     1  ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
     2  etax(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
     3  z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
     4  la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
     5  mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu),iic(lctu,lcsu)
      dimension blabel(19), nrcr(lconsu), icnoff(lconsu)
      dimension ipq(256),lfmt(2),eta(mrcru,lsu,lconsu)
      dimension nfn(lsu)
      data a0,a1/0.0d0,1.0d0/
      irun=0
      ipq(1)=0
      do 1000 i=1,255
 1000 ipq(i+1)=ipq(i)+i
      itape=inp
 1100 continue
      read(itape,890)blabel
c**      read(itape,8900) (lfmt(i),i=1,2)
 8900 format(19a4)
 7900 format(1x,19a4)
      read(itape,*)ngen,ns,naords,ncons,ngcs,itol,icut,idummy,idummy,
c    & icp, nkoord
     1 itgeom
c**      if(itgeom.gt.0) read(inp2,999) niter
c**      if(itgeom.gt.0) inp=inp2
      if(ns.le.msu) goto 2
      write(6,970) ns
      stop
    2 if(naords.le.kaords) goto 4
      write(6,971) naords
      stop
    4 if(ncons.le.mconsu) goto 6
      write(6,972) ncons
      stop
    6 if(ngcs.le.mgcsu) goto 8
      write(6,973) ngcs
      stop
    8 if(irun.gt.0)write (inpn,790) blabel
c     if (nkoord .eq. 1) write(inpn,8770)
c8770 format(/'atomic coordinates were read from file 28'/)
      read(itape,901)nst,(nd(ist),ityp(ist),ist=1,nst)
      if(irun.gt.0)then
c**      write (inpn,7900) (lfmt(i),i=1,2)
      write (inpn,900) ngen, ns, naords, nconsn, ngcs, itol,
     1 icut, idummy, idummy, izero
      write (inpn,901) nst, (nd(ist), ityp(ist), ist=1,nst)
      endif
      if(nst.le.mstu) goto 9
      write(6,974) nst
      stop
    9 read(itape,*) idmx
      if(irun.gt.0)write(inpn,900) idmx
      if(idmx.eq.0) goto 11
      do 1346 ixj=1,idmx
      if(irun.gt.0) write(inpn,900)
 1346 read(itape,900)
   11 do 10 ist=1,nst
   10 nso(ist)=0
      do 20 iaords=1,naords
      read(itape,*)iru,(la(ir,iaords),ir=1,iru)
      if(irun.gt.0) write(inpn,900)iru,(la(ir,iaords),ir=1,iru)
      if(iru.le.mru) goto 20
      write(6,975) iru
      stop
   20 nir(iaords)=iru
      do 30 igcs=1,ngcs
      read(itape,*)icsu,ictu,maords(igcs)
      if(irun.gt.0) write (inpn,900) icsu, ictu, maords(igcs)
      if(icsu.le.mcsu) goto 26
      write(6,976) icsu
      stop
   26 if(ictu.le.mctu) goto 28
      write(6,977) ictu
      stop
   28 do 30 ics=1,icsu
      read(itape,*)(iic(ict,ics),ict=1,ictu)
      if(irun.gt.0) write (inpn,900) (iic(ict,ics), ict=1,ictu)
      do  30 ict=1,ictu
      c1=dfloat(iic(ict,ics))
      c(ict,ics,igcs)=c1
      ccc=abs(c1)
   30 if(ccc.ne.a0.and.ccc.ne.a1)c(ict,ics,igcs)=dsign(sqrt(ccc),c1)
c        read in exponents and contraction coefficients
      nconsn=0
      do 41 icons=1,ncons
c     read (itape,*) iconu, lmnp1(icons), ircru
      read (itape,*) iconu, lmnp1(icons), ircru
  903 format(3i3)
      if(ircru.eq.0) ircru=1
c     write (6,900) iconu, lmnp1(icons), ircru
      if(ircru.gt.mrcru)then
          write (6,983) ircru
          stop 'program error'
      endif
      if(iconu.gt.icomx)then
          write (6,982) iconu
          stop 'program error'
      endif
      nrcr(icons)=ircru
      ncon(icons)=iconu
c
      do 40 icon=1,iconu
      read(itape,*)zet(icon,icons),(eta(ircr,icon,icons),ircr=1,ircru)
c     write (6,84) zet(icon,icons), (eta(ircr,icon,icons),
c    1   ircr=1,ircru)
40    continue
      icnoff(icons)=nconsn
      do 100 ircr=1,ircru
      nconsn=nconsn+1
      if(irun.gt.0) write(inpn,900)iconu,lmnp1(icons),ione
      do 39 icon=1,iconu
      if(eta(ircr,icon,icons).eq.0.d0.and.ircru.eq.1.and.iconu.eq.1)
     1eta(ircr,icon,icons)=1.d0
c     write(6,2000)ircr,icon,icons,eta(ircr,icon,icons)
      if(irun.gt.0) write(inpn,806)zet(icon,icons),eta(ircr,icon,icons)
39    continue
100   continue
      nrcr(icons)=ircru
41    continue
      ngenp1=ngen+1
      ng=ngenp1
c        read in atomic labels, charges, coordinates, and basis sets
      isf=0
      isfr=0
      do 75 is=1,ns
      if(irun.eq.0) nfn(is)=0
      read(inp,904)mtype(is),nf(is),nc(is),chg(is)
      if(irun.gt.0) write (inpn,904) mtype(is), nfn(is), nc(is), chg(is)
      icu=nc(is)
      if(icu.le.mcu) goto 46
      write(6,979) icu
      stop
   46 read(inp,905)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
c     if (nkoord .eq. 1)
c    &               read(28,905)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
      if(irun.gt.0)write(inpn,905)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
      if(icu.eq.1) go to 60
      do 50 ig=2,ng
      read(itape,*)(ica(ic,is,ig),ic=1,icu)
      if(irun.gt.0) write (inpn,900) (ica(ic,is,ig), ic=1,icu)
   50 continue
   60 ifu=nf(is)
      if(ifu.eq.0) go to 75
      do 70 if=1,ifu
      isf=isf+1
      read(itape,*)mcons(isf),igcs
c     write(6,900)mcons(isf),igcs
      ircru=nrcr(mcons(isf))
      irst=icnoff(mcons(isf))
      do 200 ircr=1,ircru
      if(irun.gt.0) write(inpn,900)irst+ircr,igcs
      if(irun.eq.0) then
      nfn(is)=nfn(is)+1
      endif
200   continue
   70 continue
   75 continue
      if(irun.eq.0) then
      irun=1
      rewind inp
      rewind itape
      goto 1100
      endif
   84 format(1x,5g14.8)
  790 format(1x,19a4)
  970 format(14h msu too small,i3)
  971 format(17h kaords too small,i3)
  972 format(17h mconsu too small,i4)
  973 format(16h mgcsu too small,i4)
  974 format(15h mstu too small,i3)
  975 format(14h mru too small,i3)
  976 format(15h mcsu too small,i3)
  977 format(15h mctu too small,i3)
  979 format(14h mcu too small,i3)
  982 format(24h symnew: icomx too small,i3)
  983 format(24h symnew: mrcru too small,i3)
  806 format(2g14.8)
  890 format(19a4)
  900 format(26i3)
  901 format(i3,12(i3,a3))
  904 format(a3,2i3,f3.0)
  905 format(3f14.8)
  949 format(6(/),                                  /////68h0primitive a
     1o integrals neglected if exponential factor below 10**(-,i2,1h)/63
     2h contracted ao and so integrals neglected if value below 10**(-,i
     32,1h)/43h symmetry orbital integrals written on unit,i3)
  951 format(11h degeneracy,16x,12(i4,4x))
  952 format(6h label,21x,12(2x,a4,2x))
  953 format(1h1,12x,a4,6h atoms//15h0nuclear charge,f7.2//7h0center,12x
     1,1hx,15x,1hy,15x,1hz/(i4,7x,3f16.8))
  954 format(/9h0operator,6x,19hcenter interchanges/2x,6h1(id.),6x,24i3)
  955 format(i3,6h(gen.),5x,24i3)
  958 format(//10x,a2,9h orbitals/18h0orbital exponents,7x,7(1pg15.7))
  959 format(25h contraction coefficients,7(1pg15.7))
  960 format(10h0sym. orb.,12x,15hcenter, ao type/14x,(18(i3,1h,,a2)))
  999 format(20i4)
c2000 format(3i4,f10.4)
      return
      end
c
c
      subroutine getran(nbfao,frocc,norb,itrnao,idxxt)
c
c     this subroutine calculates the symmetry transformtion matrix p
c     for the transformation from sao in ao basis
c
      implicit real*8(a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
      parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
     &          lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
     &          lpru=18000,lsu=50,lstu=8,lconu=12,
     &          lcu=8,lccu=182)
      common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
     1   mcxu, ngcs, mdum(45), mccu, mblu, ns, ng,nsf,nst
      common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
     1  ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
     2  ccoef(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
     3  z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
     4  la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
     5  mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu)
      common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
      dimension frocc(1)
      dimension itrnao(1),nsot(14),idxxt(1)
      isfxo=0
      isfp=0
      ibfao=0
      root3=sqrt(3.0d0)
      root15=sqrt(15.d0)
      do 160 is=1,ns
      icu=nc(is)
      ifu=nf(is)
      do 140 ic=1,icu
      isfx=isfxo
      isf=isfp
      if(ifu.eq.0)go to 140
      do 120 if=1,ifu
      isf=isf+1
      icons=mcons(isf)
      jtype=lmnp1(icons)
      nbu=(jtype*(jtype+1))/2
      if(nbu*icu.ne.nct(isf))go to 920
      isft=isfx+(ic-1)*nbu
      isfx=isfx+nbu*icu
      do 100 nb=1,nbu
      ibfao=ibfao+1
      isft=isft+1
      itrnao(isft)=ibfao
      idxxt(ibfao)=0
      if(jtype.le.2) goto 100
      if(jtype.eq.3.and.nb.le.3) idxxt(ibfao)=1
      if(jtype.eq.4.and.nb.gt.3.and.nb.le.9) idxxt(ibfao)=1
      if(jtype.eq.4.and.nb.le.3) idxxt(ibfao)=2
  100 continue
  120 continue
  140 continue
       isfxo=isfx
      isfp=isfp+ifu
  160 continue
      nsotx=0
      do 180 ist=1,nst
      nsot(ist)=nsotx
      nsotx=nsotx+nso(ist)
  180 continue
      nbfso=nsotx
      if(nsotx.ne.nbfso)go to 920
      do 200 i=1,nbfao
      do 200 j=1,nbfso
      p(i,j)=0.0d0
  200 continue
      isf=0
      isfr=0
      isftx=0
      do 260 is=1,ns
      ifu=nf(is)
      if(ifu.eq.0)go to 260
      do 255 if=1,ifu
      isf=isf+1
      igcs=mgcs(isf)
      iaords=maords(igcs)
      iru=nir(iaords)
      ictu=nct(isf)
      ics=0
      do 250 ir=1,iru
      isfr=isfr+1
      isft=isftx
      lai=la(ir,iaords)
      nsot(lai)=nsot(lai)+1
      ibfso=nsot(lai)
      iau=nd(lai)
      if(iau.ne.1)go to 920
      do 250 ia=1,iau
      ics=ics+1
      do 248 ict=1,ictu
      isft=isft+1
      c1=c(ict,ics,igcs)
      ibfao=itrnao(isft)
c
c     integral routine (r.pitzer) and property package use different
c     normalizations of d(xx),d(yy),d(zz),f(xxy),f(xxz),f(yyx),f(yyz),
c     f(zzx),f(zzy) -factor=sqrt(3)- , and for f(xxx),f(yyy),f(zzz)
c     -factor=sqrt(15)-
c
      if(idxxt(ibfao).eq.1) c1=c1*root3
      if(idxxt(ibfao).eq.2) c1=c1*root15
      p(ibfao,ibfso)=c1
  248 continue
  250 continue
      isftx=isftx+ictu
  255 continue
  260 continue
      call scfinp(nbfao,frocc,norb)
      return
  920 print 925
  925 format(' nct(isf) is not properly set')
      stop
      end
c
c
      subroutine scfinp(nbfao,frocc,nnorb)
c
c     this subroutine reads the data as written out by routine outgrd
c     of the scf-program (file ft08)
c
      implicit real*8(a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
      parameter (mloc=250,lstu=8,lsfru=250)
      common/eig/eig(mloc)
      common/supp/nsupp
      common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
      common/scfinf/cm(2,mloc),norbs(lstu),nc(lstu),nlamda(lstu),
     1 no(lstu),ityp(lstu)
      common /contrc/ chomo,ncont
      dimension frocc(1)
      common/holl/blabel(19),alabel(19),fmt(4),jl,lomao,iopen,nsym
      data ir/ 8/,ncol/10/, ika/2hc(/, a0 /0.0d0/
c  number of ireps
      mstu=lstu
c  lsfru as in integral program
c  total number of functions , sum(norb) , lsfru=nomx
      msfru=lsfru
c  max number of functions/irep , max(norb)
      mxbfpi=lsfru
      mxcoef=nomx*nomx
c  iopen = -1   : ci or mcscf case
c  iopen =  0   : closed shell case
c  iopen =  1   : open shell case
      ipmax=0
      icmax=0
      idmax=0
      if(nsym.le.mstu)go to 41
      print 3300,mstu
      stop
   41 continue
      do 106 l=1,nsym
      read (ir,1001) ityp(l)
      read (ir,6410) lll,norbs(l),nc(l),no(l)
      nlamda(l)=lll
      norb=norbs(l)
      if(norb.eq.0) goto 107
      ipmin = ipmax + 1
      ipmax = ipmax + norb
      if(ipmax.le.msfru)go to 58
      write(6,2800) msfru
      stop
   58 read (ir,6430) (cm(1,kl), kl=ipmin,ipmax)
      if(iopen.le.0) go to 52
      read (ir,6430) (cm(2,kl), kl=ipmin,ipmax)
      go to 56
   52 do 54 kl=ipmin,ipmax
   54 cm(2,kl)=a0
   56 icl=norbs(l)
      jcl=nlamda(l)
      icmin=icmax+1
      icmax = icmax + icl*jcl
      if(icmax.le.mxcoef) goto 63
      print 2700, mxcoef
      stop
   63 continue
  237 icx = icmin - 1
      do 224 ic=1,icl
      icxi = icx + 1
      icx = icx + jcl
  224 read (ir,fmt) (su(jc), jc=icxi,icx)
   66 continue
      if(iopen.lt.0) goto 67
      read(ir,fmt) (eig(ie),ie=ipmin,ipmax)
   67 if(nsupp.eq.2) print 1100
      if(nsupp.eq.2) print 2000, ityp(l)
      do 105 mk = 1, icl, ncol
      kin = mk + ncol - 1
      if (kin .gt. icl) kin = icl
      if(nsupp.eq.2) print 5200,(ika,i, i=mk,kin)
      nlow = icmin + jcl*mk - jcl - 1
      nhigh = icmin + kin *jcl - 1
      do 510 i1 = 1, jcl
      nlow = nlow + 1
      if(nsupp.eq.2) print 3500, (su(ic), ic=nlow,nhigh,lll)
  510 continue
      klmin=ipmin+mk-1
      klmax=ipmin+kin-1
      if(nsupp.eq.2.and.iopen.gt.-1)print 2500,(cm(1,kl),kl=klmin,klmax)
      if(nsupp.eq.2.and.iopen.eq.-1)print 2600,(cm(1,kl),kl=klmin,klmax)
      if(nsupp.eq.2.and.iopen.gt.0) print 5000,(cm(2,kl),kl=klmin,klmax)
 105  continue
      goto 106
  107 read(ir,6430) dummy
      if(iopen.gt.0) read(ir,6430) dummy
      read(ir,fmt) (dummy,ie=1,lll)
      read(ir,fmt) dummy
  106 continue
  273 continue
      do 109 i=1,ipmax
  109 chomo=chomo-cm(1,i)-cm(2,i)
      write(6,6700) chomo
      print 6500, alabel, blabel
      call comtrn(nsym,p,nbfao,nnorb,frocc)
      if(ir.ne.5)rewind ir
      return
c                        format statements
 1000 format(19a4)
 1001 format(a4)
 1100 format(47h1symmetry, basis, and configuration information)
 2000 format(/1h0,6x,6hirrep ,a4,19x,20horbital coefficients)
 2500 format(30h0number of electrons occupying/
     1   24h closed shell orbitals =,f6.1,9f11.1)
 2600 format(30h0number of electrons occupying/
     1   24h natural orbitals      =,f8.5,9f11.5)
 2700 format(38h0vector components number greater than,i5)
 2800 format(1x,' more than ',i5,' mos (lmos), change dimensions in',
     1 ' common eig and scfinf',/, 'for dimensions for localization',
     2 ' see subroutine output')
 3000 format(44h0number of orbitals in a given block exceeds, i4)
 3300 format(40h0maximum value of lambda is greater than,i3)
 3500 format(22x,10f11.6)
 5000 format(30h0number of electrons occupying/
     1   22h open shell orbitals =,f8.1,9f11.1)
 5200 format(21x,10(6x,a2,i2,1h)))
 6410 format(20i4)
 6430 format(6f12.9)
 6500 format(//,1x,19a4,/,19a4,//)
 6700 format(/,' charge of molecule',f6.2)
      end
c
c
      subroutine comtrn(nsym,p,nbfao,norb,frocc)
c
c     this subroutine transforms the eigenvectors from sao in ao basis
c     and calculates occupation numbers
c     only occupied mos are transformed
c     dimension of y  = occupied orbitals * total orbitals
c
      implicit real*8(a-h,o-z)
      parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
      parameter (mloc=250,lstu=8)
      common/scfinf/cm(2,mloc),norbs(lstu),nc(lstu),nlamda(lstu),
     1 no(lstu),ityp(lstu)
      common /big/ y(nomx*nomx),su(nomx*nomx),uv(nomx*nomx)
      common/supp/nsupp
      dimension p(nomx,1),frocc(1),idummy(100)
      data idummy/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
     1            21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,
     2            38,39,40,41,42,43,44,45,46,47,48,49,50,
     3            0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     4            0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     5            0,0,0,0,0,0,0,0,0,0/
      ipr=0
      do 17 i=1,nsym
      nn=norbs(i)
   17 ipr=ipr+nn*nbfao
      maxuv=nomx*nomx
      if(ipr.le.maxuv) goto 200
      write(6,1000) ipr,maxuv
 1000 format(1x,' dimension in comtrn too small, ipr =',i7,'  maxuv =',
     1 i7)
      stop
  200 continue
      ibfsox=0
      max=0
      iorb=0
      ic=0
      delt=1.0d-6
      nst=1
      if(nsupp.eq.3) write(6,2000)
      do 100 lam=1,nsym
      nnn=nlamda(lam)
      nn=norbs(lam)
      if(nn.eq.0) goto 100
      do 90 ii=1,nn
      iorb=iorb+1
      min=max+1
      max=max+nnn
      do 80 ibfao=1,nbfao
      ic=ic+1
      ibfso=ibfsox
      x=0.0d0
      do 70 iq=min,max
      ibfso=ibfso+1
      x=x+p(ibfao,ibfso)*su(iq)
   70 continue
      y(ic)=x
   80 continue
   90 continue
      if(nsupp.ne.3) goto 100
      write(6,3000) ityp(lam)
      call prtblk(y(nst),nbfao,nn,' ao ',' mo ',idummy)
      nst=nst+nbfao*nn
  100 ibfsox=ibfsox+nnn
c           frocc = sum of closed and open shell occupation numbers
      norb=iorb
      do 140 i=1,norb
  140 frocc(i)=cm(1,i)+cm(2,i)
c
 2000 format(/,' eigenvectors in ao basis',/)
 3000 format(/,' irrep ',a4)
      return
      end
c
c
      subroutine erduw(a,b,na,epslon)
c taken mainly from version 4, aug., 1971, of jacscf, u. of wa.
c   matrix diagonalization by the jacobi method.
c     a = real symmetric matrix to be diagonalized. it is stored
c       by columns with all sub-diagonal elements omitted, so a(i,j) is
c       stored as a((j*(j-1))/2+i).
c     b = matrix to be multiplied by the matrix of eigenvectors.
c     na = dimension of the matrices.
c      epslon is the convergence criterion for off-diagonal elements.
c     modified nov 82 r.a.
c     programm keeps track of non-diagonal norm and adjusts thresh
c     dynamically. the eventually quadratic convergence of jacobi is
c     exploited to reduce thresh faster at the end.
c
      implicit real*8(a-h,o-z)
      dimension a(1), b(1)
      data a0, a1s2 /0.0d0, 0.5d0/
      if(na.eq.1) return
      loopc=64
      sumnd=a0
      sum=a0
      ij=1
      do 24 i=1,na
      do 16 j=1,i
      term=a(ij)*a(ij)
      sum=sum+term
      if(i.eq.j) go to 16
      sumnd=sumnd+term
   16 ij=ij+1
   24 continue
      thrshg=sqrt((sum+sumnd)/na)*epslon
      small=sumnd*epslon*na*na
   32 if (sumnd.lt.small) go to 35
      thresh=sqrt(sumnd+sumnd)/na
      go to 37
   35 thresh=thresh*0.03d0
   37 thresh=dmax1(thresh,thrshg)
      n=0
      ij=2
      jj=1
      do 112 j=2,na
      jj=jj+j
      jm1=j-1
      ii=0
      do 104 i=1,jm1
      ii=ii+i
      if(abs(a(ij)).lt.thresh) go to 104
      n=n+1
      sumnd=sumnd-a(ij)*a(ij)
      sum=a1s2*(a(jj)+a(ii))
      term=a1s2*(a(jj)-a(ii))
      amax=dsign(sqrt(term*term+a(ij)*a(ij)),term)
      c=sqrt((amax+term)/(amax+amax))
      s=a(ij)/(c*(amax+amax))
      a(ii)=sum-amax
      a(jj)=sum+amax
      a(ij)=a0
      im1=i-1
      if(im1) 40,56,40
   40 ki=ii-i
      kj=jj-j
      do 48 k=1,im1
      ki=ki+1
      kj=kj+1
      term=c*a(ki)-s*a(kj)
      a(kj)=s*a(ki)+c*a(kj)
   48 a(ki)=term
   56 if(jm1.eq.i) go to 72
      ip1=i+1
      ik=ii+i
      kj=ij
      do 64 k=ip1,jm1
      kj=kj+1
      term=c*a(ik)-s*a(kj)
      a(kj)=s*a(ik)+c*a(kj)
      a(ik)=term
   64 ik=ik+k
   72 if(j.eq.na) go to 88
      jp1=j+1
      ik=jj+i
      jk=jj+j
      do 80 k=jp1,na
      term=c*a(ik)-s*a(jk)
      a(jk)=s*a(ik)+c*a(jk)
      a(ik)=term
      ik=ik+k
   80 jk=jk+k
   88 ki=im1*na
      kj=jm1*na
      do 96 k=1,na
      ki=ki+1
      kj=kj+1
      term=c*b(ki)-s*b(kj)
      b(kj)=s*b(ki)+c*b(kj)
   96 b(ki)=term
  104 ij=ij+1
  112 ij=ij+1
      loopc=loopc-1
      if(loopc) 120,1024,120
  120 if(thresh.le.thrshg.and.n.eq.0) return
      go to 32
 1024 write (6,2048)
      stop
 2048 format(12h error erduw)
      end
c
c
      subroutine outpak (matrix,nrow,nctl,isw)
c...........version = 09/05/73/03
c.......................................................................
c
c outpak prints a real*8 symmetric matrix stored in row-packed lower
c
c triangular form (see diagram below) in formatted form with numbered
c
c rows and columns.  the input is as follows:
c
c        matrix(*)...........packed matrix
c
c        nrow................number of rows to be output
c
c        nctl................carriage control flag: 1 for single space,
c                                                   2 for double space,
c                                                   3 for triple space.
c
c the matrix elements are arranged in storage as follows:
c
c
c        2
c        4    5
c        7    8    9
c       11   12   13   14
c       16   17   18   19   20
c       22   23   24   25   26   27
c
c and so on.
c
c outpak is set up to handle 8 columns/page with a 8f15.8 format
c
c for the columns.  if a different number of columns is required, change
c
c formats 1000 and 2000, and initialize kcol with the new number of
c
c columns.
c
c author:  nelson h.f. beebe, quantum theory project, university of
c          florida, gainesville
c
c.......................................................................
      real*8 matrix,column,columm
      integer begin,asa,blank,ctl
      dimension matrix(1),asa(3)
      data kcol/7/, column/8h atom nr/, asa/4h    , 4h0   , 4h-   /,
     &     blank/4h    /, zero/0.d+00/, columm/8h  lmo nr/
      ctl = blank
      if ((nctl.le.3).and.(nctl.gt.0)) ctl = asa(nctl)
c
c last is the last column number in the row currently being printed
c
      last = min0(nrow,kcol)
c
c begin is the first column number in the row currently being printed.
c
      ncol=1
c.....begin non standard do loop.
      begin=1
 1050 ncol = 1
      if(isw.eq.1) write (6,1000) (column,i,i = begin,last)
      if(isw.eq.2) write (6,1000) (columm,i,i = begin,last)
      do 40 k = begin,nrow
      ktotal = (k*(k-1))/2 + begin - 1
      do 10 i = 1,ncol
      if (matrix(ktotal+i) .ne. zero) go to 20
   10 continue
      go to 30
   20 iplk=ktotal+ncol
      kstart=ktotal+1
      if(isw.eq.1) write (6,2000) ctl,k,(matrix(ii),ii=kstart,iplk)
      if(isw.eq.2) write (6,3000) ctl,k,(matrix(ii),ii=kstart,iplk)
   30 if (k .lt. (begin+kcol-1)) ncol = ncol + 1
   40 continue
      last = min0(last+kcol,nrow)
      begin=begin+ncol
      if (begin.le.nrow) go to 1050
 1000 format (/14x,6(3x,a8,i4),4x,a6,i4)
 2000 format (a1,'atom nr',i4,2x,8f15.8)
 3000 format (a1,' lmo nr',i4,2x,8f15.8)
      return
      end
c
c
      subroutine prtblk(c,nbfn,nomx,lrow,lcol,mos)
c....................................................................
c  subroutine prtblk prints orbital coefficients by symmetry blocks.
c  it accepts as argumemts the coefficient matrix c, the number of
c  functions in the symmetry block nbfn, and the number of orbitals
c  that have been printed previoulsy nskip.
c.....................................................................
      implicit real*8(a-h,o-z)
      dimension c(nbfn,nomx),mos(1)
      data ncol/8/
      jlast=min0(nomx,ncol)
      do 400 jstrt=1,nomx,ncol
      write(6,100)(lcol,mos(j),j=jstrt,jlast)
100   format(/8x,8(6x,a4,i4,1x))
      do 300 i=1,nbfn
c
c  check for non-zero element before printing row
c
      do 120 j=jstrt,jlast
      if(c(i,j).ne.0.d0)go to 150
120   continue
c
c  exit from loop means all c's are zero.  skip the write statement
c
      go to 300
150   continue
      write(6,200)lrow,i,(c(i,j),j=jstrt,jlast)
200   format(1x,a4,i4,8f15.6)
300   continue
      jlast=min0(jlast+ncol,nomx)
400   continue
      return
      end
      subroutine orbval(eta,nbmx,norbmx,nomx,nstart,nstop,contrk,
     1ftype,frocc,c,d,dinx,cin,v,ddsq,kpmax,norbpl,mdens,nocc,y,iplot)
      implicit real*8 (a-h,o-z)
      character*4 ftype
      common/contrc/chomo,ncont
      dimension y(2),c(kpmax,3),d(kpmax,3),dinx(kpmax)
      dimension cin(kpmax,norbmx),v(kpmax),ddsq(kpmax),a(3),ftype(1)
      dimension eta(nbmx,25),frocc(1),
     1nstart(1),nstop(1),contrk(nomx,6),nocc(2)
      norb1=norbpl+1
      do 15 k=1,kpmax
      do 15 i=1,norb1
   15 cin(k,i)=0.0
      do 50 i=1,ncont
      istart=nstart(i)
      istop=nstop(i)
      iprim=0
      do 51 k=1,kpmax
   51 dinx(k)=0.0
      do 40 ii=istart,istop
      iprim=iprim+1
      do 11 l=1,3
   11 a(l)=eta(ii,l)
      do 100 k=1,kpmax
      ddsq(k)=0.0
      do 25 j=1,3
      d(k,j)=c(k,j)-a(j)
   25 ddsq(k)=ddsq(k)+d(k,j)**2
      ddsq(k)=-ddsq(k)*eta(ii,4)
      if(ddsq(k).lt.-290.) go to 40
      ddsq(k)=exp(ddsq(k))*contrk(i,iprim)
  100 v(k)=1.0
      if(ftype(ii).eq.'s   ') go to 68
      if(ftype(ii).eq.'x   ') go to 2
      if(ftype(ii).eq.'y   ') go to 3
      if(ftype(ii).eq.'z   ') go to 4
      if(ftype(ii).eq.'xx  ') go to 5
      if(ftype(ii).eq.'yy  ') go to 6
      if(ftype(ii).eq.'zz  ') go to 7
      if(ftype(ii).eq.'xy  ') go to 8
      if(ftype(ii).eq.'xz  ') go to 9
      if(ftype(ii).eq.'yz  ') go to 10
    2 do 102 k=1,kpmax
  102 v(k)=d(k,1)
      go to 68
    3 do 103 k=1,kpmax
  103 v(k)=d(k,2)
      go to 68
    4 do 104 k=1,kpmax
  104 v(k)=d(k,3)
      go to 68
    5 do 105 k=1,kpmax
  105 v(k)=d(k,1)**2
      go to 68
    6 do 106 k=1,kpmax
  106 v(k)=d(k,2)**2
      go to 68
    7 do 107 k=1,kpmax
  107 v(k)=d(k,3)**2
      go to 68
    8 do 108 k=1,kpmax
  108 v(k)=d(k,1)*d(k,2)
      go to 68
    9 do 109 k=1,kpmax
  109 v(k)=d(k,1)*d(k,3)
      go to 68
   10 do 110 k=1,kpmax
  110 v(k)=d(k,2)*d(k,3)
   68 do 69 k=1,kpmax
   69 dinx(k)=v(k)*ddsq(k)+dinx(k)
   40 continue
      do 70 nn=1,norbpl
      n=nocc(nn)
      kk=ncont*(n-1)+i
      do 80 k=1,kpmax
   80 cin(k,nn)=dinx(k)*y(kk)+cin(k,nn)
   70 continue
   50 continue
      if(mdens.gt.0) go to 150
      do 850 n=1,norbpl
c     write(iplot) nocc(n)
  850 write(iplot,99)(cin(k,n),k=1,kpmax)
 99   format(9f8.5)
      if(mdens.eq.0) return
  150 do 90 n=1,norbpl
      do 90 k=1,kpmax
      cin(k,n)=cin(k,n)*cin(k,n)
   90 cin(k,norb1)=cin(k,norb1)+cin(k,n)*frocc(n)
      do 870 n=1,norb1
c     write(iplot) nocc(n)
  870 write(iplot,99)(cin(k,n),k=1,kpmax)
c     do 860 k=1,kpmax
c 860 write(6,861) cin(k,norbpl),(c(k,i),i=1,3)
  861 format(2x,4f16.8)
      return
      end
      subroutine plotin(nmaxx,nmaxy,kpmax,plane,norbpl,mdens,
     1nocc,norbmx,iplot,iat)
      implicit real*8 (a-h,o-z)
      dimension plane(2),nocc(2)
      common /vec/ vecx(3),vecy(3),center(3)
c
c     plane   cartesian coordinates for each grid point
c     norbpl  number of orbitals to be plotted
c     mdens   =0  orbital values are calculated
c             =1  orbital density and sum of densities is calculated
c     nocc    indices of orbitals to be plotted
c
c     vecx,vecy vectors defining plane for which density is
c             calculated
c     center   center of the plane
c     dxmax    boundary of drawing in x direction
c     dymax    boundary of drawing in y direction
c     nx       number of points in x direction
c     ny       number of points in y direction
c
      read(5,100)vecx
  100 format(8f10.0)
      read(5,100)vecy
      read(5,100)center
      read(5,105)dxmax,dymax,nx,ny
  105 format(2f10.0,2i4)
      write(6,200)
  200 format('0plot input'//' vectors defining plane of plot'//
     1' original'/)
      write(6,210)vecx,vecy
  210 format(/' x direction: ',3f10.3,' y direction:',3f10.3)
      write(6,220)
  220 format(' normalized')
      sumx=0.
      sumy=0.
      do 1000 i=1,3
      sumx=sumx+vecx(i)**2
      sumy=sumy+vecy(i)**2
 1000 continue
      sumx=sqrt(sumx)
      sumy=sqrt(sumy)
      do 1010 i=1,3
      vecx(i)=vecx(i)/sumx
      vecy(i)=vecy(i)/sumy
 1010 continue
      write(6,210)vecx,vecy
      write(6,225)nx,ny
  225 format(' number of points: x direction',i4,
     1' y direction',i4)
      write(6,230)center
  230 format(' origin: ',3f10.5)
      write(6,240)dxmax,dymax
  240 format(/' boundaries: +/- xmax ',f10.5,' +/-ymax ',f10.5/)
      if(nx.gt.nmaxx)nx=nmaxx
      if(ny.gt.nmaxy)ny=nmaxy
      if(mod(nx,2).eq.0)nx=nx-1
      if(mod(ny,2).eq.0)ny=ny-1
      write(6,250)nx,ny
  250 format(' number of grid points: nx ',i4,' ny ',i4)
      write(iplot,99) nx,ny
      write(iplot,98) dxmax,dymax
   99 format(2i4)
   98 format(2f10.5)
      kpmax=nx*ny
      distx=2.*dxmax/(nx-1)
      disty=2.*dymax/(ny-1)
      call gridxy(plane,nx,ny,distx,disty,vecx,vecy,center,dxmax,dymax)
      read(5,260)norbpl,mdens
  260 format(20i4)
      read(5,260)(nocc(i),i=1,norbpl)
      if(mdens.eq.0)write(6,280)norbpl
  280 format(1x,i4,' orbitalcontours will be plotted'/
     1' orbital numbers:')
      if(mdens.eq.1)write(6,285)norbpl
  285 format(1x,i4,' orbitaldensities will be plotted'/
     1' orbital numbers:')
      if(norbpl.le.norbmx)go to 1020
      norbpl=norbmx
      write(6,270)norbpl
  270 format(' number of orbitals to be plotted reduced to',i4)
 1020 continue
      write(iplot,99) norbpl
      write(6,290)(nocc(i),i=1,norbpl)
  290 format(1x,30i4)
      read(5,260) iat
      write(iplot,99) iat
      if(iat.eq.0) write(6,300)
  300 format(/,' the atoms will be not ploted')
      if(iat.eq.1) write(6,301)
  301 format(/,' the atoms will be ploted together with the orbitals')
      if(iat.eq.2) write(6,302)
  302 format(/,' the atoms will be ploted on separat filed')
      return
      end
      subroutine gridxy(plane,nx,ny,distx,disty,vecx,vecy,center)
      implicit real*8 (a-h,o-z)
      dimension plane(nx,ny,3),vecx(3),vecy(3),center(3),start(3)
      nx1=nx/2
      ny1=ny/2
      do 10 i=1,3
      start(i)=center(i)-nx1*distx*vecx(i)-ny1*disty*vecy(i)
   10 continue
      do 40 k=1,3
      do 30 j=1,ny
      do 20 i=1,nx
      plane(i,j,k)=start(k)+(j-1)*distx*vecx(k)+(i-1)*disty*vecy(k)
   20 continue
   30 continue
   40 continue
      return
      end
      subroutine molcor(iplot)
      implicit real*8(a-h,o-z)
      parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
     &          lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
     &          lpru=18000,lsu=50,lstu=8,lconu=12,
     &          lcu=8,lccu=182)
      common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
     1mcxu,ngcs,mdum(45),mccu,mblu,ns,ng,nsf,nst,ncons
      common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
     1  ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
     2    eta(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
     3  z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
     4  la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
     5  mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu),iic(lctu,lcsu)
      common/vec/ vecx(3),vecy(3),center(3)
      parameter(maxatm=100)
      dimension xatom(maxatm),yatom(maxatm)
      natom=0
      do 75 is=1,ns
      icu=nc(is)
      do 76 ic=1,icu
      natom=natom+1
      if(natom.gt.maxatm) stop 'maxatm is not eneught'
      xatom(natom)=(x(ic,is)-center(1))*vecx(1)
     +            +(yy(ic,is)-center(2))*vecx(2)
     +            +(z(ic,is)-center(3))*vecx(3)
      yatom(natom)=(x(ic,is)-center(1))*vecy(1)
     +            +(yy(ic,is)-center(2))*vecy(2)
     +            +(z(ic,is)-center(3))*vecy(3)
   76 continue
   75 continue
      write(iplot,99) natom,(xatom(i),yatom(i),i=1,natom)
   99 format(i4,/,(9f8.5))
      return
      end