c-------------------------------------------------------------------------------
      subroutine RESPA_step(itime)
c-------------------------------------------------------------------------------
c  Perform a single RESPA step.
c-------------------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CONTROL'
      include 'COMMON.VAR'
      include 'COMMON.MD'
      include 'COMMON.LANGEVIN'
      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(0:n_ene),energia_short(0:n_ene),
     & energia_long(0:n_ene)
      double precision cm(3),L(3),vcm(3),incr(3)
      integer ilen,count,rstcount
      external ilen
      character*50 tytul
      integer maxcount_scale /10/
      common /gucio/ cm
      double precision stochforcvec(MAXRES6)
      double precision grad_tmp(3,0:maxres2)
      common /stochcalc/ stochforcvec
      integer itime
      logical scale
      double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
      if (large.and. mod(itime,ntwe).eq.0) then
        write (iout,*) "***************** RESPA itime",itime
        write (iout,*) "Cartesian and internal coordinates: step 0"
        call cartprint
        call intout
        write (iout,*) "Accelerations"
        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
c
c Perform the initial RESPA step (increment velocities)
c      write (iout,*) "*********************** RESPA ini"

      if (tnp1) then
creview          call tnp1_respa_step1
          call tnp_respa_step1
      else if (tnp) then
          call tnp_respa_step1
      else
          if (tnh.and..not.xiresp) then
            call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
            do i=0,2*nres
             do j=1,3
              d_t(j,i)=d_t(j,i)*scale_nh
             enddo
            enddo 
          endif
          call RESPA_vel
      endif

cd       if(tnp .or. tnp1) then
cd        write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
cd       endif

      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
c 
cr      if (tnp) then
cr          do i=0,nres
cr            do j=1,3
cr             grad_tmp(j,i)=gcart(j,i)
cr             grad_tmp(j,i+nres)=gxcart(j,i)
cr            enddo
cr          enddo
cr      endif
c
c Compute the short-range forces
      call zerograd
      call etotal_short(energia_short)
      if (tnp.or.tnp1) potE=energia_short(0)
      call cartgrad
      call lagrangian
      do i=0,2*nres
        do j=1,3
          dc_old(j,i)=dc(j,i)
          if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
          d_a_old(j,i)=d_a(j,i)
        enddo
      enddo 
      d_time0=d_time
c Split the time step
      d_time=d_time/ntime_split 
ctest      E_long2=E_long
c Perform the short-range RESPSA steps (velocity Verlet increments of
c positions and velocities using short-range forces)
c      write (iout,*) "*********************** RESPA split"
creview      E_old=potE
      do itsplit=1,ntime_split
        if (lang.eq.1) then
          call sddir_precalc
        else if (lang.eq.2 .or. lang.eq.3) then
          call stochastic_force(stochforcvec)
        endif
c First step of the velocity Verlet algorithm
        if (lang.eq.2) then
          call sd_verlet1
        else if (lang.eq.3) then
          call sd_verlet1_ciccotti
        else if (lang.eq.1) then
          call sddir_verlet1
        else if (tnp1) then
          call tnp1_respa_i_step1
cr          if(itsplit.eq.1)then
cr           d_time_s12=d_time0*0.5*s12_np
cr           do j=1,3
cr            d_t_half(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
cr           enddo
cr           do i=nnt,nct-1
cr            do j=1,3
cr             d_t_half(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
cr            enddo
cr           enddo
cr           do i=nnt,nct
cr            if (itype(i).ne.10) then
cr             inres=i+nres
cr             do j=1,3
cr              d_t_half(j,inres)=d_t_old(j,inres)
cr     &                +d_a_old(j,inres)*d_time_s12
cr             enddo
cr            endif
cr           enddo
cr          endif
        else if (tnp) then
          call tnp_respa_i_step1
        else
          if (tnh.and.xiresp) then
            call kinetic(EK)
            call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
            do i=0,2*nres
             do j=1,3
              d_t_old(j,i)=d_t_old(j,i)*scale_nh
             enddo
            enddo 
cd            write(iout,*) "SSS1",itsplit,EK,scale_nh
          endif
          call verlet1
        endif	  
c Build the chain from the newly calculated coordinates	
        call chainbuild_cart
        if (rattle) call rattle1
        if (large.and. mod(itime,ntwe).eq.0) then
          write (iout,*) "Cartesian and internal coordinates: step 1"
          call cartprint
          call intout
          write (iout,*) "Accelerations"
          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 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
        tt0 = tcpu()
c Calculate energy and forces
c
c E_long aproximation
cr         if (tnp .or. tnp1) then
cr          dtmp=0.5*d_time*(1.0/s12_np+1.0/s_np)
cr          do i=0,2*nres
cr            do j=1,3
cr             E_long=E_long+d_t_new(j,i)*dtmp*grad_tmp(j,i)
cr            enddo
cr          enddo
cr         endif
c-------------------------------------
c test of reviewer's comment
cr        E_long=0
c-------------------------------------
 
         
c
ctest        call etotal_long(energia_long)
ctest        E_long=energia_long(0)
ctest
        call zerograd

        call etotal_short(energia_short)
        E_old=potE
        potE=energia_short(0)

c       if(tnp .or. tnp1) then
c        write (iout,*) "kkk",E_long2,E_long
c       endif

          
        call cartgrad
c Get the new accelerations
        call lagrangian
        t_enegrad=t_enegrad+tcpu()-tt0
c Second step of the velocity Verlet algorithm
        if (lang.eq.2) then	
          call sd_verlet2
        else if (lang.eq.3) then
          call sd_verlet2_ciccotti
        else if (lang.eq.1) then
          call sddir_verlet2
        else if (tnp1) then
            call tnp1_respa_i_step2
        else if (tnp) then
            call tnp_respa_i_step2
        else
          call verlet2
          if (tnh.and.xiresp) then
            call kinetic(EK)
            call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
            do i=0,2*nres
             do j=1,3
              d_t(j,i)=d_t(j,i)*scale_nh
             enddo
            enddo 
cd            write(iout,*) "SSS2",itsplit,EK,scale_nh
          endif

        endif 	    	    
        if (rattle) call rattle2
c Backup the coordinates, velocities, and accelerations
        if (tnp .or. tnp1) then 
         do i=0,2*nres
          do j=1,3
            d_t_old(j,i)=d_t(j,i)
            if (tnp) d_t(j,i)=d_t(j,i)/s_np
            if (tnp1) d_t(j,i)=d_t(j,i)/s_np
          enddo
         enddo 

        endif
        do i=0,2*nres
          do j=1,3
            dc_old(j,i)=dc(j,i)
            if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
            d_a_old(j,i)=d_a(j,i)
          enddo
        enddo 

cd       if(tnp .or. tnp1) then
cd        call kinetic(EK)
cd        HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
cd        H=(HNose1-H0)*s_np
cd        write (iout,*) "jjj",EK,potE
cd        write (iout,*) "iii H=",H,abs(HNose1-H0)/H0
cd        write (iout,'(a,3f)') 
cd     &             "III NP S, pi",totT+itsplit*d_time, s_np, pi_np
cd       endif


      enddo
c Restore the time step
      d_time=d_time0
c Compute long-range forces
      call zerograd
      call etotal_long(energia_long)
      E_long=energia_long(0)
      potE=energia_short(0)+energia_long(0)
      call cartgrad
c Compute accelerations from long-range forces
      call lagrangian
      if (large.and. mod(itime,ntwe).eq.0) then
        write (iout,*) "Cartesian and internal coordinates: step 2"
        call cartprint
        call intout
        write (iout,*) "Accelerations"
        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
c Compute the final RESPA step (increment velocities)
c      write (iout,*) "*********************** RESPA fin"
      if (tnp1) then
creview          call tnp1_respa_step2
          call tnp_respa_step2
      else if (tnp) then
          call tnp_respa_step2
      else
          call RESPA_vel
          if (tnh.and..not.xiresp) then
            call kinetic(EK)
            call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
            do i=0,2*nres
             do j=1,3
              d_t(j,i)=d_t(j,i)*scale_nh
             enddo
            enddo 
          endif
      endif

        if (tnp .or. tnp1) then 
         do i=0,2*nres
          do j=1,3
            d_t(j,i)=d_t_old(j,i)/s_np
          enddo
         enddo 
        endif

c Compute the complete potential energy
cc      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/(dimen*Rb)*EK
c Backup the coordinates, velocities, and accelerations
      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

c-----review
c       if(tnp .or. tnp1) then
c        HNose1=Hnose(EK,s_np,energia_short(0),pi_np,Q_np,t_bath,dimen)
c_new_var_csplit         Csplit=H0-E_long
c         Csplit=H0-energia_short(0)
c       endif
c----------


      if (mod(itime,ntwe).eq.0) then

       if(tnp .or. tnp1) then
        write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit,
     &                          E_long,energia_short(0)
        HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
        H=(HNose1-H0)*s_np
cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen*0.001986d0*t_bath*log(s_np)
        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
cd        write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
       endif

       if(tnh) then
        HNose1=Hnose_nh(EK,potE)
        H=HNose1-H0
        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
       endif


       if (large) then
       itnp=0
       do j=1,3
        itnp=itnp+1
        vtnp(itnp)=d_t(j,0)
       enddo
       do i=nnt,nct-1	
        do j=1,3    
          itnp=itnp+1
          vtnp(itnp)=d_t(j,i)
        enddo
       enddo
       do i=nnt,nct
        if (itype(i).ne.10) then
          inres=i+nres
          do j=1,3  
           itnp=itnp+1  
           vtnp(itnp)=d_t(j,inres)
          enddo
        endif      
       enddo 

c Transform velocities from UNRES coordinate space to cartesian and Gvec
c eigenvector space

        do i=1,dimen
          vtnp_(i)=0.0d0
          vtnp_a(i)=0.0d0
          do j=1,dimen
            vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
            vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
          enddo
          vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
        enddo

        do i=1,dimen
         write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
        enddo

       endif
      endif
      return
      end