      subroutine move(iswitch,istep,nreal,ninteger,nmaxpoints,ndim1,
     &                ndim2,natoms,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,stepsize,dmaxmove,coordsave,
     &                coord,charge,dmata,dmatb,coeffa,coeffb,alpha,
     &                coeffs,cdfitc,alphacd,coefscd,coefpcd,coefdcd,
     &                grad,overlap,wkvec,lsym)

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"

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

      dimension coordsave(3,*),coord(3,*),charge(*),dmata(*),dmatb(*),
     &          coeffa(ndim1,*),coeffb(ndim1,*),alpha(*),coeffs(*),
     &          cdfitc(*),alphacd(*),coefscd(*),coefpcd(*),coefdcd(*),
     &          grad(3,*),overlap(ndim2,*),wkvec(*)

      character*4 lsym(*)

      data zero,big/0.0d0,999999.9d0/

      save ifollow

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

      igeometry=istep

      call gradmm(nreal,ninteger,ndim2,natoms,nconts,ncontp,ncontd,ncds,
     &            ncdspd,nshels,nshelp,nsheld,mtloca,ilfunc,ngaussians,
     &            icfunc,icdlfunc,icdcfunc,iwkvec,coord,charge,alpha,
     &            coeffs,cdfitc,alphacd,coefscd,coefpcd,coefdcd,grad,
     &            overlap,dmata,dmatb,wkvec(1),wkvec(ndim2+1))

      call cdgrmm(natoms,ncds,ncdspd,icdlfunc,icdcfunc,cdfitc,alphacd,
     &            coefscd,coefpcd,coefdcd,coord,grad,wkvec)

      call gradxc(nreal-6*natoms,nmaxpoints,ngridtype,iptsweight,
     &            ncontractions,natoms,nconts,ncontp,ncontd,natomtype,
     &            nptsatom,nfuncatom,ilfunc,icfunc,ngaussians,nshels,
     &            nshelp,nsheld,iwkvec,wkvec,wkvec(6*natoms+1),dmata,
     &            dmatb,coord,grad,alpha,coeffs)

      ia=1
      ib=ia+ndim2
      ic=ib+nalpha*nalpha
      id=ic+ndim2
      ie=id+ndim2

      if(ie.gt.nreal) stop 'not enough room for pulay forces'

      call pulayforce(ndim1,ndim2,natoms,nalpha,nbeta,nconts,ncontp,
     &                ncontd,icfunc,mtloca,nshels,nshelp,nsheld,grad,
     &                wkvec(ia),overlap,coeffa,coeffb,wkvec(ib),
     &                wkvec(ic),wkvec(id))

      if((iswitch.eq.1).or.(myid.gt.0)) return

 1001 continue

      if(igeometry.eq.1) then

            ia=1
            ib=ia+natoms*natoms
            ic=ib+natoms*natoms*natoms
            id=ic+natoms*natoms*natoms
            ie=id+natoms*natoms*natoms

            call find_ints(natoms,nbonds,nangles,ndihedrals,nimpropers,
     &                     natomtype,iwkvec(ia),iwkvec(ib),iwkvec(ic),
     &                     iwkvec(id),iwkvec(ie),coord)

            rewind 31

            do 1002 i=1,ibfrozen
            iwkvec(ia+2*(nbonds+i-1)+0)=ibfreeze(1,i)
 1002       iwkvec(ia+2*(nbonds+i-1)+1)=ibfreeze(2,i)

            do 1003 i=1,iafrozen
            iwkvec(ib+3*(nangles+i-1)+0)=iafreeze(1,i)
            iwkvec(ib+3*(nangles+i-1)+1)=iafreeze(2,i)
 1003       iwkvec(ib+3*(nangles+i-1)+2)=iafreeze(3,i)

            do 1004 i=1,idfrozen
            iwkvec(ic+4*(ndihedrals+i-1)+0)=idfreeze(1,i)
            iwkvec(ic+4*(ndihedrals+i-1)+1)=idfreeze(2,i)
            iwkvec(ic+4*(ndihedrals+i-1)+2)=idfreeze(3,i)
 1004       iwkvec(ic+4*(ndihedrals+i-1)+3)=idfreeze(4,i)

            do 1005 i=1,iifrozen
            iwkvec(id+4*(nimpropers+i-1)+0)=iifreeze(1,i)
            iwkvec(id+4*(nimpropers+i-1)+1)=iifreeze(2,i)
            iwkvec(id+4*(nimpropers+i-1)+2)=iifreeze(3,i)
 1005       iwkvec(id+4*(nimpropers+i-1)+3)=iifreeze(4,i)

            nbonds=nbonds+ibfrozen
            nangles=nangles+iafrozen
            ndihedrals=ndihedrals+idfrozen
            nimpropers=nimpropers+iifrozen

            if(imode.eq.1) then
                                 iwkvec(ia+2*nbonds+0)=ibmode(1)
                                 iwkvec(ia+2*nbonds+1)=ibmode(2)

                                 nbonds=nbonds+1
                           endif

            if(imode.eq.2) then
                                 iwkvec(ib+3*nangles+0)=iamode(1)
                                 iwkvec(ib+3*nangles+1)=iamode(2)
                                 iwkvec(ib+3*nangles+2)=iamode(3)

                                 nangles=nangles+1
                           endif

            if(imode.eq.3) then
                                 iwkvec(ic+4*ndihedrals+0)=idmode(1)
                                 iwkvec(ic+4*ndihedrals+1)=idmode(2)
                                 iwkvec(ic+4*ndihedrals+2)=idmode(3)
                                 iwkvec(ic+4*ndihedrals+3)=idmode(4)

                                 ndihedrals=ndihedrals+1
                           endif

            if(imode.eq.4) then
                                 iwkvec(id+4*nimpropers+0)=iimode(1)
                                 iwkvec(id+4*nimpropers+1)=iimode(2)
                                 iwkvec(id+4*nimpropers+2)=iimode(3)
                                 iwkvec(id+4*nimpropers+3)=iimode(4)

                                 nimpropers=nimpropers+1
                           endif

            if(imode.eq.1) ifollow=nbonds
            if(imode.eq.2) ifollow=nbonds+nangles
            if(imode.eq.3) ifollow=nbonds+nangles+ndihedrals
            if(imode.eq.4) ifollow=nbonds+nangles+ndihedrals+nimpropers

            write(31) nbonds,nangles,ndihedrals,nimpropers

            write(31) (iwkvec(ia+i),i=0,2*nbonds-1)
            write(31) (iwkvec(ib+i),i=0,3*nangles-1)
            write(31) (iwkvec(ic+i),i=0,4*ndihedrals-1)
            write(31) (iwkvec(id+i),i=0,4*nimpropers-1)

                         endif

      rewind 31

      read(31) nbonds,nangles,ndihedrals,nimpropers

      nmodes=nbonds+nangles+ndihedrals+nimpropers

      ia=1
      ib=ia+1+nbonds*2
      ic=ib+1+nangles*3
      id=ic+1+ndihedrals*4
      ie=id+1+nimpropers*4

      read(31) (iwkvec(ia+i),i=0,2*nbonds-1)
      read(31) (iwkvec(ib+i),i=0,3*nangles-1)
      read(31) (iwkvec(ic+i),i=0,4*ndihedrals-1)
      read(31) (iwkvec(id+i),i=0,4*nimpropers-1)

      ja=1
      jb=ja+nmodes
      jc=jb+nmodes
      jd=jc+nmodes
      je=jd+natoms*3
      jf=je+nmodes*natoms*3
      jg=jf+nmodes*nmodes
      jh=jg+nmodes*nmodes
      ji=jh+natoms*3
      jj=ji+nmodes
      jk=jj+nmodes
      jl=jk+nmodes*nmodes
      jm=jl+nmodes*natoms*3
      jn=jm+nmodes*nmodes
      jo=jn+nmodes*nmodes
      jp=jo+nmodes*nmodes
      jq=jp+natoms*3

      ifrozen=ibfrozen+iafrozen+idfrozen+iifrozen

      ii=ie-1

      kk=nbonds-ibfrozen

      do 1006 i=nbonds-ibfrozen+1,nbonds
      ii=ii+1
      kk=kk+1
 1006 iwkvec(ii)=kk

      kk=nbonds+nangles-iafrozen

      do 1007 i=nangles-iafrozen+1,nangles
      ii=ii+1
      kk=kk+1
 1007 iwkvec(ii)=kk

      kk=nbonds+nangles+ndihedrals-idfrozen

      do 1008 i=ndihedrals-idfrozen+1,ndihedrals
      ii=ii+1
      kk=kk+1
 1008 iwkvec(ii)=kk

      kk=nbonds+nangles+ndihedrals+nimpropers-iifrozen

      do 1009 i=nimpropers-iifrozen+1,nimpropers
      ii=ii+1
      kk=kk+1
 1009 iwkvec(ii)=kk

      if(imode.eq.0) then

               call redundant(igeometry,natoms,nmodes,nbonds,nangles,
     &                        ndihedrals,nimpropers,ifrozen,natomtype,
     &                        iwkvec(ia),iwkvec(ib),iwkvec(ic),
     &                        iwkvec(id),iwkvec(ie),stepsize,grms,gmax,
     &                        crms,cmax,coord,wkvec(ja),wkvec(jb),
     &                        wkvec(jc),wkvec(jd),wkvec(je),wkvec(jf),
     &                        wkvec(jg),grad,wkvec(ji),wkvec(jj),
     &                        wkvec(jk),wkvec(jl),wkvec(jm),wkvec(jn),
     &                        wkvec(jo),wkvec(jp),wkvec(jq))

                     else

               call tstate(igeometry,ifollow,natoms,nmodes,nbonds,
     &                     nangles,ndihedrals,nimpropers,ifrozen,
     &                     natomtype,iwkvec(ia),iwkvec(ib),iwkvec(ic),
     &                     iwkvec(id),iwkvec(ie),stepsize,grms,gmax,
     &                     crms,cmax,coord,wkvec(ja),wkvec(jb),
     &                     wkvec(jc),wkvec(jd),wkvec(je),wkvec(jf),
     &                     wkvec(jg),grad,wkvec(ji),wkvec(jj),wkvec(jk),
     &                     wkvec(jl),wkvec(jm),wkvec(jn),wkvec(jo),
     &                     wkvec(jp),wkvec(jq))

                     endif

      if(igeometry.eq.0) then

             igeometry=1

             write(6,1010)
 1010        format(//,' redefining the redundant internal',
     &                 ' coordinates, angle approaching 180!!',//)

             goto 1001

                         endif

      write(6,1011) grms,gmax,crms,cmax
 1011 format(/,' rms component of gradient =      ',f10.6,/,
     &         ' maximum component of gradient =  ',f10.6,/,
     &         ' rms change in coordinates =      ',f10.6,/,
     &         ' maximum change in coordinates =  ',f10.6,/////)

      if((grms.lt.0.00030).and.
     &   (gmax.lt.0.00045).and.
     &   (crms.lt.0.00120).and.
     &   (cmax.lt.dmaxmove))    then
                                      dmaxfactor=zero

                                      call stringcopy(3*natoms,
     &                                                coord,wkvec(jp))
                                else
                                      dmaxfactor=big
                                endif

      return
      end
