c-------------------------------------------------------------------------------
      subroutine RESPA_step(itime)
c-------------------------------------------------------------------------------
c  Perform a single RESPA step.
c-------------------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
      integer IERROR,ERRCODE
#endif
      include 'COMMON.SETUP'
      include 'COMMON.CONTROL'
      include 'COMMON.VAR'
      include 'COMMON.MD'
#ifndef LANG0
      include 'COMMON.LANGEVIN'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.IOUNITS'
      include 'COMMON.NAMES'
      include 'COMMON.TIME1'
      double precision energia_short(0:n_ene),
     & energia_long(0:n_ene)
      double precision cm(3),L(3),vcm(3),incr(3)
      double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
     & d_a_old0(3,0:maxres2)
      integer ilen,count,rstcount
      external ilen
      character*50 tytul
      integer maxcount_scale /10/
      common /gucio/ cm
      double precision stochforcvec(MAXRES6)
      common /stochcalc/ stochforcvec
      integer itime
      logical scale
      common /cipiszcze/ itt
      itt=itime
      if (ntwe.ne.0) then
      if (large.and. mod(itime,ntwe).eq.0) then
        write (iout,*) "***************** RESPA itime",itime
        write (iout,*) "Cartesian and internal coordinates: step 0"
c        call cartprint
        call pdbout(0.0d0,"cipiszcze",iout)
        call intout
        write (iout,*) "Accelerations from long-range forces"
        do i=0,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
     &      (d_a(j,i+nres),j=1,3)
        enddo
        write (iout,*) "Velocities, step 0"
        do i=0,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
     &      (d_t(j,i+nres),j=1,3)
        enddo
      endif
      endif
c
c Perform the initial RESPA step (increment velocities)
c      write (iout,*) "*********************** RESPA ini"
      call RESPA_vel
      if (ntwe.ne.0) then
      if (mod(itime,ntwe).eq.0 .and. large) then
        write (iout,*) "Velocities, end"
        do i=0,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
     &      (d_t(j,i+nres),j=1,3)
        enddo
      endif
      endif
c Compute the short-range forces
#ifdef MPI
      tt0 =MPI_Wtime()
#else
      tt0 = tcpu()
#endif
C 7/2/2009 commented out
c      call zerograd
c      call etotal_short(energia_short)
c      call cartgrad
c      call lagrangian
C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
        do i=0,2*nres
          do j=1,3
            d_a(j,i)=d_a_short(j,i)
          enddo
        enddo
      if (ntwe.ne.0) then
      if (large.and. mod(itime,ntwe).eq.0) then
        write (iout,*) "energia_short",energia_short(0)
        write (iout,*) "Accelerations from short-range forces"
        do i=0,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
     &      (d_a(j,i+nres),j=1,3)
        enddo
      endif
      endif
#ifdef MPI
        t_enegrad=t_enegrad+MPI_Wtime()-tt0
#else
        t_enegrad=t_enegrad+tcpu()-tt0
#endif
      do i=0,2*nres
        do j=1,3
          dc_old(j,i)=dc(j,i)
          d_t_old(j,i)=d_t(j,i)
          d_a_old(j,i)=d_a(j,i)
        enddo
      enddo 
c 6/30/08 A-MTS: attempt at increasing the split number
      do i=0,2*nres
        do j=1,3
          dc_old0(j,i)=dc_old(j,i)
          d_t_old0(j,i)=d_t_old(j,i)
          d_a_old0(j,i)=d_a_old(j,i)
        enddo
      enddo 
      if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
      if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
c
      scale=.true.
      d_time0=d_time
      do while (scale)

      scale=.false.
c      write (iout,*) "itime",itime," ntime_split",ntime_split
c Split the time step
      d_time=d_time0/ntime_split 
c Perform the short-range RESPA steps (velocity Verlet increments of
c positions and velocities using short-range forces)
c      write (iout,*) "*********************** RESPA split"
      do itsplit=1,ntime_split
        if (lang.eq.1) then
          call sddir_precalc
        else if (lang.eq.2 .or. lang.eq.3) then
#ifndef LANG0
          call stochastic_force(stochforcvec)
#else
          write (iout,*) 
     &      "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
#ifdef MPI
          call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
#endif
          stop
#endif
        endif
c First step of the velocity Verlet algorithm
        if (lang.eq.2) then
#ifndef LANG0
          call sd_verlet1
#endif
        else if (lang.eq.3) then
#ifndef LANG0
          call sd_verlet1_ciccotti
#endif
        else if (lang.eq.1) then
          call sddir_verlet1
        else
          call verlet1
        endif
c Build the chain from the newly calculated coordinates	
        call chainbuild_cart
        if (rattle) call rattle1
        if (ntwe.ne.0) then
        if (large.and. mod(itime,ntwe).eq.0) then
          write (iout,*) "***** ITSPLIT",itsplit
          write (iout,*) "Cartesian and internal coordinates: step 1"
          call pdbout(0.0d0,"cipiszcze",iout)
c          call cartprint
          call intout
          write (iout,*) "Velocities, step 1"
          do i=0,nres
            write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
     &        (d_t(j,i+nres),j=1,3)
          enddo
        endif
        endif
#ifdef MPI
        tt0 = MPI_Wtime()
#else
        tt0 = tcpu()
#endif
c Calculate energy and forces
        call zerograd
        call etotal_short(energia_short)
#ifdef TIMING_ENE
#ifdef MPI
        t_eshort=t_eshort+MPI_Wtime()-tt0
#else
        t_eshort=t_eshort+tcpu()-tt0
#endif
#endif
        call cartgrad
c Get the new accelerations
        call lagrangian
C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
        do i=0,2*nres
          do j=1,3
            d_a_short(j,i)=d_a(j,i)
          enddo
        enddo
        if (ntwe.ne.0) then
        if (large.and. mod(itime,ntwe).eq.0) then
          write (iout,*)"energia_short",energia_short(0)
          write (iout,*) "Accelerations from short-range forces"
          do i=0,nres
            write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
     &        (d_a(j,i+nres),j=1,3)
          enddo
        endif
        endif
c 6/30/08 A-MTS
c Determine maximum acceleration and scale down the timestep if needed
        call max_accel
        amax=amax/ntime_split**2
        call predict_edrift(epdrift)
        if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) 
     &   write (iout,*) "amax",amax," damax",damax,
     &   " epdrift",epdrift," epdriftmax",epdriftmax
c Exit loop and try with increased split number if the change of
c acceleration is too big
        if (amax.gt.damax .or. epdrift.gt.edriftmax) then
          if (ntime_split.lt.maxtime_split) then
            scale=.true.
            ntime_split=ntime_split*2
            do i=0,2*nres
              do j=1,3
                dc_old(j,i)=dc_old0(j,i)
                d_t_old(j,i)=d_t_old0(j,i)
                d_a_old(j,i)=d_a_old0(j,i)
              enddo
            enddo 
            write (iout,*) "acceleration/energy drift too large",amax,
     &      epdrift," split increased to ",ntime_split," itime",itime,
     &       " itsplit",itsplit
            exit
          else
            write (iout,*) 
     &      "Uh-hu. Bumpy landscape. Maximum splitting number",
     &       maxtime_split,
     &      " already reached!!! Trying to carry on!"
          endif
        endif
#ifdef MPI
        t_enegrad=t_enegrad+MPI_Wtime()-tt0
#else
        t_enegrad=t_enegrad+tcpu()-tt0
#endif
c Second step of the velocity Verlet algorithm
        if (lang.eq.2) then
#ifndef LANG0
          call sd_verlet2
#endif
        else if (lang.eq.3) then
#ifndef LANG0
          call sd_verlet2_ciccotti
#endif
        else if (lang.eq.1) then
          call sddir_verlet2
        else
          call verlet2
        endif
        if (rattle) call rattle2
c Backup the coordinates, velocities, and accelerations
        do i=0,2*nres
          do j=1,3
            dc_old(j,i)=dc(j,i)
            d_t_old(j,i)=d_t(j,i)
            d_a_old(j,i)=d_a(j,i)
          enddo
        enddo 
      enddo

      enddo ! while scale

c Restore the time step
      d_time=d_time0
c Compute long-range forces
#ifdef MPI
      tt0 =MPI_Wtime()
#else
      tt0 = tcpu()
#endif
      call zerograd
      call etotal_long(energia_long)
#ifdef TIMING_ENE
#ifdef MPI
        t_elong=t_elong+MPI_Wtime()-tt0
#else
        t_elong=t_elong+tcpu()-tt0
#endif
#endif
      call cartgrad
      call lagrangian
#ifdef MPI
        t_enegrad=t_enegrad+MPI_Wtime()-tt0
#else
        t_enegrad=t_enegrad+tcpu()-tt0
#endif
c Compute accelerations from long-range forces
      if (ntwe.ne.0) then
      if (large.and. mod(itime,ntwe).eq.0) then
        write (iout,*) "energia_long",energia_long(0)
        write (iout,*) "Cartesian and internal coordinates: step 2"
c        call cartprint
        call pdbout(0.0d0,"cipiszcze",iout)
        call intout
        write (iout,*) "Accelerations from long-range forces"
        do i=0,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
     &      (d_a(j,i+nres),j=1,3)
        enddo
        write (iout,*) "Velocities, step 2"
        do i=0,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
     &      (d_t(j,i+nres),j=1,3)
        enddo
      endif
      endif
c Compute the final RESPA step (increment velocities)
c      write (iout,*) "*********************** RESPA fin"
      call RESPA_vel
c Compute the complete potential energy
      do i=0,n_ene
        potEcomp(i)=energia_short(i)+energia_long(i)
      enddo
      potE=potEcomp(0)-potEcomp(20)
c      potE=energia_short(0)+energia_long(0)
      totT=totT+d_time
c Calculate the kinetic and the total energy and the kinetic temperature
      call kinetic(EK)
      totE=EK+potE
c Couple the system to Berendsen bath if needed
      if (tbf .and. lang.eq.0) then
        call verlet_bath
      endif
      kinetic_T=2.0d0/(dimen3*Rb)*EK
c Backup the coordinates, velocities, and accelerations
      if (ntwe.ne.0) then
      if (mod(itime,ntwe).eq.0 .and. large) then
        write (iout,*) "Velocities, end"
        do i=0,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
     &      (d_t(j,i+nres),j=1,3)
        enddo
      endif
      endif
      return
      end