      subroutine deft(nreal,ninteger,nbuffer,nmaxatoms,nmaxcontractions,
     &                ndim,nauxfunctions,nmaxpoints,nesppts,ncontacts,
     &                iwkvec,nptsatom,natomtype,izmat1,izmat2,izmat3,
     &                neutral,isys,nequivvec,nfuncatom,ilfunc,icfunc,
     &                mtloca,ngaussians,nshels,nshelp,nsheld,icdlfunc,
     &                icdcfunc,ixclfunc,ixccfunc,ibfreeze,iafreeze,
     &                idfreeze,iifreeze,wkvec,coord,coordsave,charge,
     &                coordptcharge,ptcharge,grad,coeffa,coeffb,overlap,
     &                dmta,dmtb,core,focka,fockb,focke,alpha,coeff,
     &                alphacd,coefscd,coefpcd,coefdcd,alphaxc,coefsxc,
     &                coefpxc,coefdxc,cdfitc,xcfita,xcfitb,xcfite,
     &                xcfitsave,tvector,tveca,tvecb,tvece,cutoff,buffer,
     &                inputstring,atomsymbol)

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

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

      include "mpif.h"

      character inputstring(nmaxatoms,*)*30,atomsymbol(*)*4

      dimension ibmode(2),iamode(3),idmode(4),iimode(4)

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

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

      data zero,one,two,three,four,five,seven/0.0d0,1.0d0,2.0d0,3.0d0,
     &                                        4.0d0,5.0d0,7.0d0/
      data smallnumber,bignumber/0.000001d0,1.0d+25/
      data bohr,evolt/0.529177d0,27.2116d0/

      call startup

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      call mpi_buffer_attach(buffer,nbuffer*8,ierror)

      if(myid.gt.0) goto 1014

      imode=0
      imesp=0
      irestart=0
      ioptimize=0
      ncenters=0
      nptcharges=0
      iptsflag=1
      icharge=0
      imultiplicity=1
      ngridtype=5454
      nfunctional=0
      idiis=1
      idiffuse=0
      ixcfit=0
      iptsweight=1
      nmaxiterations=50
      nmaxgeometries=30
      ivibrations=0
      ilocalize=0
      nfdpbt=0
      nequivalences=0
      nbneutrals=0
      nbsystems=1
      level=1
      ibfrozen=0
      iafrozen=0
      idfrozen=0
      iifrozen=0

      do 1001 j=1,nmaxatoms
      do 1001 i=1,nmaxatoms
      isys(i+1,1)=i
      nequivvec(i,1)=0
      nequivvec(i,2)=0
 1001 neutral(i,j)=0

      dmixing=0.100
      dmaxstep=0.300
      econvergence=0.00000001
      gconvergence=0.00180
      elevelshift=zero
      epsilon=one
      cutoff(1)=five
      cutoff(2)=seven
      cutoff(3)=bignumber

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

      if(nbsystems.gt.1) nbsystems=nbsystems-1

      if(nbsystems.eq.1) isys(1,1)=ncenters

      write(6,1002)
 1002 format(//)

      dnuclearcharge=zero

      do 1003 i=1,ncenters
 1003 dnuclearcharge=dnuclearcharge+charge(i)

      nelectrons=nint(dnuclearcharge)-icharge

      if((mod(imultiplicity,2)).eq.1) then

                  nalpha=(nelectrons/2)+((imultiplicity-1)/2)
                  nbeta=(nelectrons/2)-((imultiplicity-1)/2)

                                      else

                  nalpha=((nelectrons+1)/2)+((imultiplicity-2)/2)
                  nbeta=((nelectrons-1)/2)-((imultiplicity-2)/2)

                                      endif

      if((nalpha+nbeta).ne.nelectrons) then

             write(6,1004)
 1004        format(' the charge and multiplicity do not make sense')
             stop

                                       endif

      savedmixing=dmixing

      call inputbases(nmaxatoms,ncenters,nconts,ncontp,ncontd,
     &                ncontractions,ncds,ncdsets,ncdfunctions,nxcs,
     &                nxcsets,nxcfunctions,idiffuse,ixcfit,wkvec,
     &                wkvec(1*nauxfunctions+1),wkvec(2*nauxfunctions+1),
     &                wkvec(3*nauxfunctions+1),icfunc,ilfunc,ngaussians,
     &                mtloca,nshels,nshelp,nsheld,icdcfunc,icdlfunc,
     &                ixccfunc,ixclfunc,natomtype,alpha,coeff,alphacd,
     &                coefscd,coefpcd,coefdcd,alphaxc,coefsxc,coefpxc,
     &                coefdxc,inputstring)

      do 1005 i=1,ncdfunctions
 1005 cdfitc(i)=zero

      icds=0
      icdspd=ncds

      do 1006 i=1,ncenters
 1006 if(natomtype(i).ne.0) call start(natomtype(i),icds,icdspd,cdfitc)

      do 1007 i=1,nxcfunctions
      xcfita(i)=-(one/three)*cdfitc(i)
      xcfitb(i)=-(one/three)*cdfitc(i)
 1007 xcfite(i)=-(one/three)*(three/four)*cdfitc(i)

      ja=1
      jb=ja+9*ncenters**2
      jc=jb+9*ncenters**2
      jd=jc+3*ncenters
      je=jd+3*ncenters
      jf=je+3*ncenters

      if(ioptimize.ne.0) then
                               do 1008 j=1,ncenters
                               do 1008 i=1,3
 1008                          coordsave(i,j)=coord(i,j)
                         endif

      if(irestart.ne.0) call readt(ncenters,ncdfunctions,nxcfunctions,
     &                             nmaxcontractions,ncontractions,coord,
     &                             cdfitc,xcfita,xcfitb,xcfite,coeffa,
     &                             coeffb)

      if(intcoordflag.eq.1) call zmatrix(ncenters,coord,izmat1,izmat2,
     &                                   izmat3,atomsymbol)

      write(6,1009)
 1009 format(' cartesian coordinates (with basis sets):',/)

      if(idiffuse.eq.1) write(6,1010)
 1010 format('   N.B. on heavy atoms, diffuse functions added to',
     &         ' the following bases ',/)
      if(idiffuse.eq.2) write(6,1011)
 1011 format('   N.B. on all atoms, diffuse functions added to',
     &         ' the following bases ',/)

      do 1012 i=1,ncenters
 1012 write(6,1013) atomsymbol(i),coord(1,i)*bohr,coord(2,i)*bohr,
     &              coord(3,i)*bohr,inputstring(i,1),inputstring(i,2)
 1013 format(1x,a4,2x,3f12.7,6x,a30,/,49x,a30,/)

      if(irestart.eq.3) goto 1028

      call rijmat(ncenters,atomsymbol,coord)

 1014 continue

      call mpi_bcast(ioptimize,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ivibrations,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(imesp,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(irestart,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(intcoordflag,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(iptsflag,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ixcfit,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(iptsweight,1,mpi_integer,0,mpi_comm_world,ierr)

      if(ivibrations.eq.1) then

          if(myid.eq.0) call savet(ncenters,ncdfunctions,nxcfunctions,
     &                             nmaxcontractions,ncontractions,coord,
     &                             cdfitc,xcfita,xcfitb,xcfite,coeffa,
     &                             coeffb)

          call vibrations(nreal,ninteger,nmaxatoms,nmaxcontractions,
     &                    nauxfunctions,nmaxpoints,ncenters,
     &                    nptcharges,nbsystems,level,nalpha,nbeta,
     &                    nconts,ncontp,ncontd,ncds,ncdsets-ncds,nxcs,
     &                    nxcsets-nxcs,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,coordsave,
     &                    charge,grad,coordptcharge,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,tvector,
     &                    tveca,tvecb,tvece,atomsymbol)

            return

                      endif

      if(ioptimize.eq.0) then

          call scfcontrol(nfdpbt,nreal,ninteger,nmaxatoms,
     &                    nmaxcontractions,nauxfunctions,nmaxpoints,
     &                    ncontacts,nesppts,nmaxiterations,nfunctional,
     &                    idiis,ngridtype,iptsflag,ncenters,nptcharges,
     &                    nbsystems,level,nalpha,nbeta,nconts,ncontp,
     &                    ncontd,ncds,ncdsets-ncds,nxcs,nxcsets-nxcs,
     &                    ixcfit,natomtype,nptsatom,nshels,nshelp,
     &                    nsheld,ilfunc,icfunc,ngaussians,mtloca,
     &                    nfuncatom,icdlfunc,icdcfunc,ixclfunc,ixccfunc,
     &                    isys,iwkvec,epsilon,econvergence,dmixing,
     &                    elevelshift,cutoff,coord,charge,coordptcharge,
     &                    ptcharge,core,focka,fockb,focke,dmta,dmtb,
     &                    coeffa,coeffb,overlap,alpha,coeff,alphacd,
     &                    coefscd,coefpcd,coefdcd,alphaxc,coefsxc,
     &                    coefpxc,coefdxc,cdfitc,xcfita,xcfitb,xcfite,
     &                    xcfitsave,tvector,tveca,tvecb,tvece,wkvec)

          goto 1028

                         endif

      gradfactor=bignumber

      call mpi_bcast(nmaxgeometries,1,mpi_integer,0,mpi_comm_world,ierr)

      do 1015 igeom=1,nmaxgeometries

      call mpi_bcast(gradfactor,1,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(gconvergence,1,mpi_real8,0,mpi_comm_world,ierr)

      if(gradfactor.lt.gconvergence) goto 1015

      if(myid.eq.0) then
                          write(6,1016) igeom
 1016                     format(////,' geometry # ',i4,/)

                          if(intcoordflag.eq.1) call zmatrix(ncenters,
     &                                                       coord,
     &                                                       izmat1,
     &                                                       izmat2,
     &                                                       izmat3,
     &                                                       atomsymbol)

                          write(6,1017)
 1017                     format(' cartesian coordinates :',/)

                          do 1018 i=1,ncenters
 1018                     write(6,1019) atomsymbol(i),coord(1,i)*bohr,
     &                                                coord(2,i)*bohr,
     &                                                coord(3,i)*bohr
 1019                     format(5x,a4,2x,3f15.7)

                          write(6,1020)
 1020                     format()
                    endif

      call scfcontrol(nfdpbt,nreal,ninteger,nmaxatoms,nmaxcontractions,
     &                nauxfunctions,nmaxpoints,ncontacts,nesppts,
     &                nmaxiterations,nfunctional,idiis,ngridtype,
     &                iptsflag,ncenters,nptcharges,nbsystems,level,
     &                nalpha,nbeta,nconts,ncontp,ncontd,ncds,
     &                ncdsets-ncds,nxcs,nxcsets-nxcs,ixcfit,natomtype,
     &                nptsatom,nshels,nshelp,nsheld,ilfunc,icfunc,
     &                ngaussians,mtloca,nfuncatom,icdlfunc,icdcfunc,
     &                ixclfunc,ixccfunc,isys,iwkvec,epsilon,
     &                econvergence,dmixing,elevelshift,cutoff,coord,
     &                charge,coordptcharge,ptcharge,core,focka,fockb,
     &                focke,dmta,dmtb,coeffa,coeffb,overlap,alpha,coeff,
     &                alphacd,coefscd,coefpcd,coefdcd,alphaxc,coefsxc,
     &                coefpxc,coefdxc,cdfitc,xcfita,xcfitb,xcfite,
     &                xcfitsave,tvector,tveca,tvecb,tvece,wkvec)

      if(ixcfit.eq.0) then
                            nxcs=0
                            nxcsets=0
                      endif

      call mpi_bcast(ncontractions,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(savedmixing,1,mpi_real8,0,mpi_comm_world,ierr)

      dmixing=savedmixing

      nsqdim=ncontractions*(ncontractions+1)/2
      ja=3*nsqdim

      call mpi_bcast(ncdsets,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nfuncatom,ncontractions,
     &               mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(izmat1,ncenters,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(izmat2,ncenters,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(izmat3,ncenters,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(dmaxstep,1,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(coordsave,3*ncenters,
     &               mpi_real8,0,mpi_comm_world,ierr)

      call move(0,igeom,nreal-ja,ninteger,nmaxpoints,nmaxcontractions,
     &          nsqdim,ncenters,nalpha,nbeta,ncontractions,nconts,
     &          ncontp,ncontd,ncds,ncdsets-ncds,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,gradfactor,dmaxstep,
     &          gconvergence,coordsave,coord,charge,dmta,dmtb,coeffa,
     &          coeffb,alpha,coeff,cdfitc,alphacd,coefscd,coefpcd,
     &          coefdcd,grad,wkvec(1),wkvec(ja+1),atomsymbol)

 1015 continue

      if(gradfactor.lt.gconvergence) then
                                           if(myid.eq.0) write(6,1021)
                                     else
                                           if(myid.eq.0) write(6,1022)
                                     endif

 1021 format(/,' geometry optimized!!! :-) :-) :-) ',//)
 1022 format(/,' geometry not optimized??? :-( :-( :-( ',//)

      if(myid.eq.0) write(6,1023)
 1023 format(/,' final geometry:',/)

      if((intcoordflag.eq.1).and.(myid.eq.0)) call zmatrix(ncenters,
     &                                                     coord,
     &                                                     izmat1,
     &                                                     izmat2,
     &                                                     izmat3,
     &                                                     atomsymbol)

      if(myid.eq.0) write(6,1024)
 1024 format(' cartesian coordinates :',/)

      do 1025 i=1,ncenters
 1025 if(myid.eq.0) write(6,1026) atomsymbol(i),coord(1,i)*bohr,
     &                            coord(2,i)*bohr,coord(3,i)*bohr
 1026 format(5x,a4,2x,3f15.7)

      if(myid.eq.0) write(6,1027)
 1027 format()

 1028 continue

      if(myid.eq.0) call rijmat(ncenters,atomsymbol,coord)

      if(myid.eq.0) call writet(ncenters,ncdfunctions,nxcfunctions,
     &                          nmaxcontractions,ncontractions,coord,
     &                          cdfitc,xcfita,xcfitb,xcfite,coeffa,
     &                          coeffb)

      if(myid.ne.0) return

      ncontractions=nconts+3*ncontp+6*ncontd

      ja=1
      jb=ja+3*nesppts
      jc=jb+4*nmaxatoms*nmaxatoms
      jd=jc+2*nmaxatoms
      je=jd+nmaxatoms
      jf=je+3*nmaxatoms
      jg=jf+3*nmaxatoms
      jh=jg+3*ncontacts

      if(imesp.eq.1) then

               do 2001 i=1,ncontractions
               do 2001 j=1,i
               imax=max(i,j)
               imin=min(i,j)
               ij=((imax*(imax-1))/2)+imin
               temp_value=zero
               do 2002 k=1,nalpha
 2002          temp_value=temp_value+coeffa(i,k)*coeffa(j,k)
               if(i.ne.j) temp_value=two*temp_value
 2001          dmta(ij)=temp_value

               do 3001 i=1,ncontractions
               do 3001 j=1,i
               imax=max(i,j)
               imin=min(i,j)
               ij=((imax*(imax-1))/2)+imin
               temp_value=zero
               do 3002 k=1,nalpha
 3002          temp_value=temp_value+coeffb(i,k)*coeffb(j,k)
               if(i.ne.j) temp_value=two*temp_value
 3001          dmtb(ij)=temp_value

               call mespgen(nreal-jh,ninteger-2*nesppts,iptsflag,
     &                      nmaxatoms,nesppts,ncontacts,nauxfunctions,
     &                      level,ncenters,nbsystems,nequivalences,
     &                      nbneutrals,nalpha,nbeta,nconts,ncontp,
     &                      ncontd,ncds,ncdsets-ncds,nshels,nshelp,
     &                      nsheld,ilfunc,ngaussians,icfunc,mtloca,
     &                      icdlfunc,icdcfunc,natomtype,nequivvec,
     &                      neutral,isys,iwkvec,iwkvec(nesppts+1),
     &                      iwkvec(2*nesppts+1),coord,charge,dmta,dmtb,
     &                      alpha,coeff,alphacd,coefscd,coefpcd,coefdcd,
     &                      cdfitc,wkvec(ja),wkvec(jb),wkvec(jc),
     &                      wkvec(jd),wkvec(je),wkvec(jf),wkvec(jg),
     &                      wkvec(jh))

                     endif

      if(irestart.eq.3) return

      ncontractions=nconts+3*ncontp+6*ncontd

      if(myid.eq.0) then
                          ncntsq=(ncontractions*(ncontractions+1))/2

                          rewind 43

                          call fastrd(43,overlap,ncntsq)
                    endif

      if(5*(ncontractions**2)+1.gt.nreal) stop 'mulliken'

      call mulliken(ncontractions,ncenters,nconts,ncontp,ncontd,nshels,
     &              nshelp,nsheld,icfunc,nfuncatom,natomtype,overlap,
     &              dmta,dmtb,wkvec(0*(ncontractions**2)+1),
     &              wkvec(1*(ncontractions**2)+1),
     &              wkvec(2*(ncontractions**2)+1),
     &              wkvec(3*(ncontractions**2)+1),
     &              wkvec(4*(ncontractions**2)+1))

      ja=1
      jb=ja+ncontractions*ncontractions
      jc=jb+ncontractions*ncontractions
      jd=jc+ncontractions*ncontractions
      je=jd+ncontractions*ncontractions
      jf=je+ncontractions
      jg=jf+ncontractions

      if(ilocalize.ge.2) then

                   call boys_mm(nreal-jf,ninteger,ncontractions,nconts,
     &                          ncontp,ncontd,nshels,nshelp,nsheld,
     &                          mtloca,ilfunc,ngaussians,icfunc,iwkvec,
     &                          coord,alpha,coeff,wkvec(ja),wkvec(jb),
     &                          wkvec(jc),wkvec(jd))

                   call boys_loc(ncontractions,nmaxcontractions,nalpha,
     &                           coeffa,wkvec(ja),wkvec(jb),wkvec(jc),
     &                           wkvec(jd),wkvec(je),wkvec(jf),
     &                           wkvec(jg))

                   call orbitals(ncontractions,nmaxcontractions,nalpha,
     &                           nconts,ncontp,ncontd,ilfunc,icfunc,
     &                           ngaussians,nshels,nshelp,nsheld,
     &                           nfuncatom,coord,alpha,coeff,overlap,
     &                           coeffa,wkvec(ja),wkvec(jg),wkvec(jb),
     &                           wkvec(jc))

                   write(6,1029)
 1029              format(///,' boys localized orbitals generated ',
     &                        '(spin alpha)',/)

                   do 1030 iii=1,nalpha

                   write(6,1031) iii
 1031              format('   orbital # ',i5)

                   call boys_print(ncontractions,ncenters,nconts,ncontp,
     &                             ncontd,nshels,nshelp,nsheld,icfunc,
     &                             nfuncatom,natomtype,overlap,
     &                             coeffa(1,iii),
     &                             wkvec(0*(ncontractions**2)+1),
     &                             wkvec(1*(ncontractions**2)+1),
     &                             wkvec(2*(ncontractions**2)+1))

 1030              continue

                   if(nbeta.eq.nalpha) goto 1035

                   call boys_mm(nreal-jf,ninteger,ncontractions,nconts,
     &                          ncontp,ncontd,nshels,nshelp,nsheld,
     &                          mtloca,ilfunc,ngaussians,icfunc,iwkvec,
     &                          coord,alpha,coeff,wkvec(ja),wkvec(jb),
     &                          wkvec(jc),wkvec(jd))

                   call boys_loc(ncontractions,nmaxcontractions,nbeta,
     &                           coeffb,wkvec(ja),wkvec(jb),wkvec(jc),
     &                           wkvec(jd),wkvec(je),wkvec(jf),
     &                           wkvec(jg))

                   call orbitals(ncontractions,nmaxcontractions,nbeta,
     &                           nconts,ncontp,ncontd,ilfunc,icfunc,
     &                           ngaussians,nshels,nshelp,nsheld,
     &                           nfuncatom,coord,alpha,coeff,overlap,
     &                           coeffb,wkvec(ja),wkvec(jg),wkvec(jb),
     &                           wkvec(jc))

                   write(6,1032)
 1032              format(///,' boys localized orbitals generated ',
     &                        '(spin beta)',/)

                   do 1033 iii=1,nbeta

                   write(6,1034) iii
 1034              format('   orbital # ',i5)

                   call boys_print(ncontractions,ncenters,nconts,ncontp,
     &                             ncontd,nshels,nshelp,nsheld,icfunc,
     &                             nfuncatom,natomtype,overlap,
     &                             coeffb(1,iii),
     &                             wkvec(0*(ncontractions**2)+1),
     &                             wkvec(1*(ncontractions**2)+1),
     &                             wkvec(2*(ncontractions**2)+1))

 1033              continue

 1035              continue

                         endif

      if((ilocalize.gt.0).and.(ilocalize.le.2)) then

                   write(6,1036)
 1036              format(///,' canonical orbitals generated ',
     &                        '(spin alpha)')

                   ka=1
                   kb=ka+ncontractions
                   kc=kb+ncontractions*ncontractions
                   kd=kc+ncontractions*ncontractions
                   ke=kd+ncontractions*ncontractions
                   kf=ke+ncontractions

                   call canonical(nmaxcontractions,ncontractions,nalpha,
     &                            nconts,ncontp,ncontd,nfuncatom,iwkvec,
     &                            focka,coeffa,wkvec(ka),wkvec(kb),
     &                            wkvec(kc),wkvec(kd),wkvec(ke),
     &                            wkvec(kf))

                   if(nbeta.eq.nalpha) goto 1038

                   write(6,1037)
 1037              format(///,' canonical orbitals generated ',
     &                        '(spin beta)')

                   ka=1
                   kb=ka+ncontractions
                   kc=kb+ncontractions*ncontractions
                   kd=kc+ncontractions*ncontractions
                   ke=kd+ncontractions*ncontractions
                   kf=ke+ncontractions

                   call canonical(nmaxcontractions,ncontractions,nbeta,
     &                            nconts,ncontp,ncontd,nfuncatom,iwkvec,
     &                            fockb,coeffb,wkvec(ka),wkvec(kb),
     &                            wkvec(kc),wkvec(kd),wkvec(ke),
     &                            wkvec(kf))

 1038              continue

                                                endif

      na=ncontractions*(ncontractions+1)/2

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

      call mpi_buffer_detach(buffer,nbuffer*8,ierror)

      return
      end
