subroutine MD
c------------------------------------------------
c  The driver for molecular dynamics subroutines
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 cm(3),L(3),vcm(3)
      double precision energia(0:n_ene)
      integer ilen,rstcount
      external ilen
      character*50 tytul
      common /gucio/ cm,energia
      integer itime
c
      t_MDsetup=0.0d0
      t_langsetup=0.0d0
      t_MD=0.0d0
      t_enegrad=0.0d0
      t_sdsetup=0.0d0
      write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
      tt0 = tcpu()
c Determine the inverse of the inertia matrix.
      call setup_MD_matrices
c Initialize MD
      call init_MD
      t_MDsetup = tcpu()-tt0
      rstcount=0 
c   Entering the MD loop       
      tt0 = tcpu()
      if (lang.eq.2 .or. lang.eq.3) then
        call setup_fricmat
        if (lang.eq.2) then
          call sd_verlet_p_setup	
        else
          call sd_verlet_ciccotti_setup
        endif
        do i=1,dimen
          do j=1,dimen
            pfric0_mat(i,j,0)=pfric_mat(i,j)
            afric0_mat(i,j,0)=afric_mat(i,j)
            vfric0_mat(i,j,0)=vfric_mat(i,j)
            prand0_mat(i,j,0)=prand_mat(i,j)
            vrand0_mat1(i,j,0)=vrand_mat1(i,j)
            vrand0_mat2(i,j,0)=vrand_mat2(i,j)
          enddo
        enddo
        flag_stoch(0)=.true.
        do i=1,maxflag_stoch
          flag_stoch(i)=.false.
        enddo  
      else if (lang.eq.1 .or. lang.eq.4) then
        call setup_fricmat
      endif
      t_langsetup=tcpu()-tt0
      tt0=tcpu()
      do itime=1,n_timestep
        rstcount=rstcount+1
        if (lang.gt.0 .and. surfarea .and. 
     &      mod(itime,reset_fricmat).eq.0) then
          if (lang.eq.2 .or. lang.eq.3) then
            call setup_fricmat
            if (lang.eq.2) then
              call sd_verlet_p_setup
            else
              call sd_verlet_ciccotti_setup
            endif
            do i=1,dimen
              do j=1,dimen
                pfric0_mat(i,j,0)=pfric_mat(i,j)
                afric0_mat(i,j,0)=afric_mat(i,j)
                vfric0_mat(i,j,0)=vfric_mat(i,j)
                prand0_mat(i,j,0)=prand_mat(i,j)
                vrand0_mat1(i,j,0)=vrand_mat1(i,j)
                vrand0_mat2(i,j,0)=vrand_mat2(i,j)
              enddo
            enddo
            flag_stoch(0)=.true.
            do i=1,maxflag_stoch
              flag_stoch(i)=.false.
            enddo   
          else if (lang.eq.1 .or. lang.eq.4) then
            call setup_fricmat
          endif
          write (iout,'(a,i10)') 
     &      "Friction matrix reset based on surface area, itime",itime
        endif
        if (reset_vel .and. tbf .and. lang.eq.0 
     &      .and. mod(itime,count_reset_vel).eq.0) then
          call random_vel
          write(iout,'(a,f20.2)') 
     &     "Velocities reset to random values, time",totT	
          do i=0,2*nres
            do j=1,3
              d_t_old(j,i)=d_t(j,i)
            enddo
          enddo
        endif
       	if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
          call inertia_tensor  
          call vcm_vel(vcm)
          do j=1,3
             d_t(j,0)=d_t(j,0)-vcm(j)
          enddo
          call kinetic(EK)
          kinetic_T=2.0d0/(dimen*Rb)*EK
          scalfac=dsqrt(T_bath/kinetic_T)
          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT	
          do i=0,2*nres
            do j=1,3
              d_t_old(j,i)=scalfac*d_t(j,i)
            enddo
          enddo
        endif  
        if (lang.ne.4) then
          if (RESPA) then
c Time-reversible RESPA algorithm 
c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
            call RESPA_step(itime)
          else
c Variable time step algorithm.
            call velverlet_step(itime)
          endif
        else
          call brown_step(itime)
        endif
        if (mod(itime,ntwe).eq.0) call statout(itime)
        if (mod(itime,ntwx).eq.0) then
          write (tytul,'("time",f8.2)') totT
          if(mdpdb) then
             call pdbout(potE,tytul,ipdb)
          else 
             call cartout(totT)
          endif
        endif
        if (rstcount.eq.1000.or.itime.eq.n_timestep) then
           open(irest2,file=rest2name,status='unknown')
           write(irest2,*) totT,EK,potE,totE
           do i=1,2*nres
            write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
           enddo
           do i=1,2*nres
            write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
           enddo
          close(irest2)
          rstcount=0
        endif 
      enddo
      t_MD=tcpu()-tt0
      write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') 
     &  '  Timing  ',
     & 'MD calculations setup:',t_MDsetup,
     & 'Energy & gradient evaluation:',t_enegrad,
     & 'Stochastic MD setup:',t_langsetup,
     & 'Stochastic MD step setup:',t_sdsetup,
     & 'MD steps:',t_MD
      write (iout,'(/28(1h=),a25,27(1h=))') 
     & '  End of MD calculation  '
      return
      end