c---------------------------------------------------------------------
      subroutine max_accel
c
c Find the maximum difference in the accelerations of the the sites
c at the beginning and the end of the time step.
c
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CONTROL'
      include 'COMMON.VAR'
      include 'COMMON.MD'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.IOUNITS'
      double precision aux(3),accel(3),accel_old(3),dacc
      do j=1,3
c        aux(j)=d_a(j,0)-d_a_old(j,0)
         accel_old(j)=d_a_old(j,0)
         accel(j)=d_a(j,0)
      enddo 
      amax=0.0d0
      do i=nnt,nct
c Backbone
        if (i.lt.nct) then
c 7/3/08 changed to asymmetric difference
          do j=1,3
c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
            accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
            accel(j)=accel(j)+0.5d0*d_a(j,i)
c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
            if (dabs(accel(j)).gt.dabs(accel_old(j))) then
              dacc=dabs(accel(j)-accel_old(j))
              if (dacc.gt.amax) amax=dacc
            endif
          enddo
        endif
      enddo
c Side chains
      do j=1,3
c        accel(j)=aux(j)
        accel_old(j)=d_a_old(j,0)
        accel(j)=d_a(j,0)
      enddo
      if (nnt.eq.2) then
        do j=1,3
          accel_old(j)=accel_old(j)+d_a_old(j,1)
          accel(j)=accel(j)+d_a(j,1)
        enddo
      endif
      do i=nnt,nct
        if (itype(i).ne.10) then
          do j=1,3 
c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
            accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
            accel(j)=accel(j)+d_a(j,i+nres)
          enddo
        endif
        do j=1,3
c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
          if (dabs(accel(j)).gt.dabs(accel_old(j))) then
            dacc=dabs(accel(j)-accel_old(j))
            if (dacc.gt.amax) amax=dacc
          endif
        enddo
        do j=1,3
          accel_old(j)=accel_old(j)+d_a_old(j,i)
          accel(j)=accel(j)+d_a(j,i)
c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
        enddo
      enddo
      return
      end