      subroutine vibrations(nreal,ninteger,nmaxatoms,nmaxorbs,nmaxaux,
     &                      nmaxpoints,ncenters,nptcharges,nbsystems,
     &                      level,nalpha,nbeta,nconts,ncontp,ncontd,
     &                      ncds,ncdspd,nxcs,nxcspd,idiis,
     &                      nmaxiterations,nfunctional,ngridtype,
     &                      ibfrozen,iafrozen,idfrozen,iifrozen,ixcfit,
     &                      iptsweight,imode,iwkvec,izmat1,izmat2,
     &                      izmat3,natomtype,nfuncatom,nptsatom,icfunc,
     &                      ilfunc,ngaussians,mtloca,nshels,nshelp,
     &                      nsheld,icdcfunc,icdlfunc,ixccfunc,ixclfunc,
     &                      isys,ibfreeze,iafreeze,idfreeze,iifreeze,
     &                      ibmode,iamode,idmode,iimode,econvergence,
     &                      dmixing,elevelshift,cutoff,wkvec,coord,
     &                      coordfix,charge,grad,coordpc,ptcharge,
     &                      overlap,core,focka,fockb,focke,dmta,dmtb,
     &                      coeffa,coeffb,alpha,coeff,alphacd,coefscd,
     &                      coefpcd,coefdcd,alphaxc,coefsxc,coefpxc,
     &                      coefdxc,cdfitc,xcfita,xcfitb,xcfite,
     &                      xcfitsave,tvec,tveca,tvecb,tvece,labelatoms)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

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

      character labelatoms(*)*4

      dimension iwkvec(*),izmat1(*),izmat2(*),izmat3(*),natomtype(*),
     &          nfuncatom(*),nptsatom(*),icfunc(*),ilfunc(*),
     &          ngaussians(*),mtloca(*),nshels(*),nshelp(*),nsheld(*),
     &          icdcfunc(*),icdlfunc(*),ixccfunc(*),ixclfunc(*),
     &          isys(nmaxatoms,*),ibfreeze(2,*),iafreeze(3,*),
     &          idfreeze(4,*),iifreeze(4,*),ibmode(*),iamode(*),
     &          idmode(*),iimode(*)

      dimension cutoff(*),wkvec(*),coord(3,*),coordfix(3,*),charge(*),
     &          grad(3,*),coordpc(3,*),ptcharge(*),overlap(nmaxorbs,*),
     &          core(*),focka(*),fockb(*),focke(*),dmta(*),dmtb(*),
     &          coeffa(nmaxorbs,*),coeffb(nmaxorbs,*),alpha(*),coeff(*),
     &          alphacd(*),coefscd(*),coefpcd(*),coefdcd(*),alphaxc(*),
     &          coefsxc(*),coefpxc(*),coefdxc(*),cdfitc(*),xcfita(*),
     &          xcfitb(*),xcfite(*),xcfitsave(*),tvec(*),tveca(*),
     &          tvecb(*),tvece(*)

      data small,step,bohr/0.00001d0,0.0025d0,0.529177d0/

      rewind 73
      rewind 74

      ncontractions=nconts+3*ncontp+6*ncontd

      displacement=step/bohr

      nstep=0

      do 1001 i=1,ncenters
      do 1001 j=1,3
 1001 coordfix(j,i)=coord(j,i)

      do 1002 i=1,2
      do 1002 j=1,ncenters
      do 1002 k=1,3

      if(i.eq.1) coord(k,j)=coord(k,j)+displacement
      if(i.eq.2) coord(k,j)=coord(k,j)-displacement

      nstep=nstep+1

      write(6,1003) nstep
 1003 format(///////////,' geometry number',i4,
     &                   ' of the vibrational analysis',/)

      do 1004 l=1,ncenters
 1004 write(6,1005) l,coord(1,l)*bohr,coord(2,l)*bohr,coord(3,l)*bohr
 1005 format(' atom # ',i3,5x,3f12.7)

      write(6,1006)
 1006 format(/)

      call scfrun(nreal,ninteger,nmaxatoms,nmaxorbs,nmaxaux,nmaxpoints,
     &            nfunctional,nmaxiterations,idiis,ngridtype,ncenters,
     &            nptcharges,nbsystems,level,nalpha,nbeta,nconts,ncontp,
     &            ncontd,ncds,ncdspd,nxcs,nxcspd,ixcfit,natomtype,
     &            nptsatom,nshels,nshelp,nsheld,ilfunc,icfunc,mtloca,
     &            ngaussians,nfuncatom,icdlfunc,icdcfunc,ixclfunc,
     &            ixccfunc,isys,iwkvec,econvergence,dmixing,elevelshift,
     &            diiserror,cutoff,coord,charge,coordpc,ptcharge,core,
     &            focka,fockb,focke,dmta,dmtb,coeffa,coeffb,alpha,coeff,
     &            alphacd,coefscd,coefpcd,coefdcd,alphaxc,coefsxc,
     &            coefpxc,coefdxc,cdfitc,xcfita,xcfitb,xcfite,xcfitsave,
     &            tvec,tveca,tvecb,tvece,overlap,wkvec)

      n=ncontractions*(ncontractions+1)/2

      call dipolemm(nreal-n-5,ninteger,ncenters,nconts,ncontp,ncontd,
     &              nshels,nshelp,nsheld,mtloca,ilfunc,ngaussians,
     &              icfunc,iwkvec,wkvec(1),coord,charge,alpha,coeff,
     &              dmta,dmtb,wkvec(5),wkvec(n+5))

      call fastwr(73,wkvec,3)

      ndim=(ncontractions*(ncontractions+1))/2

      ja=3*ndim

      call move(1,nstep,nreal-ja,ninteger,nmaxpoints,nmaxorbs,ndim,
     &          ncenters,nalpha,nbeta,ncontractions,nconts,ncontp,
     &          ncontd,ncds,ncdspd,ibfrozen,iafrozen,idfrozen,iifrozen,
     &          imode,ngridtype,iptsweight,natomtype,nptsatom,nfuncatom,
     &          izmat1,izmat2,izmat3,nshels,nshelp,nsheld,ilfunc,icfunc,
     &          ngaussians,mtloca,icdlfunc,icdcfunc,ibfreeze,iafreeze,
     &          idfreeze,iifreeze,ibmode,iamode,idmode,iimode,iwkvec,
     &          dmaxfactor,step,small,coordfix,coord,charge,dmta,dmtb,
     &          coeffa,coeffb,alpha,coeff,cdfitc,alphacd,coefscd,
     &          coefpcd,coefdcd,grad,wkvec(1),wkvec(ja+1),labelatoms)

      write(6,1007)
 1007 format(/,' gradient :',/)

      do 1008 l=1,ncenters
 1008 write(6,1009) l,grad(1,l),grad(2,l),grad(3,l)
 1009 format(' atom # ',i3,5x,3f17.12)

      do 1010 l=1,ncenters
 1010 write(74) grad(1,l),grad(2,l),grad(3,l)

      do 1011 l=1,ncenters
      do 1011 m=1,3
 1011 coord(m,l)=coordfix(m,l)

 1002 continue

      ja=1
      jb=ja+ncenters*ncenters*9
      jc=jb+ncenters*ncenters*9
      jd=jc+ncenters*9
      je=jd+ncenters*9
      jf=je+ncenters*3
      jg=jf+ncenters*9
      jh=jg+ncenters*3
      ji=jh+ncenters*3

      call irspectrum(ncenters,natomtype,displacement,coord,wkvec(ja),
     &                wkvec(jb),wkvec(jc),wkvec(jd),wkvec(je),wkvec(jf),
     &                wkvec(jg),wkvec(jh),wkvec(ji))

      return
      end
