c---------------------------------------------------------
      subroutine init_MD
c  Set up the initial conditions of a MD simulation
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MP
      include 'mpif.h'
      include 'COMMON.INFO'
      include 'COMMON.SETUP'
      character*4 liczba
      character*16 form
#endif
      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'
      real*8 energia(0:n_ene),energia_long(0:n_ene),
     &  energia_short(0:n_ene),vcm(3),incr(3)
      double precision cm(3),L(3),xv,sigv,lowb,highb
      character*256 qstr
      integer ilen
      external ilen
      character*50 tytul
      common /gucio/ cm
      d_time0=d_time
      write(iout,*) "d_time", d_time
c Compute the standard deviations of stochastic forces for Langevin dynamics
c if the friction coefficients do not depend on surface area
      if (lang.gt.0 .and. .not.surfarea) then
        do i=nnt,nct-1
          stdforcp(i)=stdfp*dsqrt(gamp)
        enddo
        do i=nnt,nct
          stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
        enddo
      endif
c Open the pdb file for snapshotshots
#ifdef MPI
      if (nprocs.eq.1) then
        npos=3
      else
        npos = dlog10(dfloat(nprocs-1))+1
      endif
      if (npos.lt.3) npos=3
      write (liczba,'(i1)') npos
      form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
     &  //')'
      write (liczba,form) myrank
      if(mdpdb) then
        open(ipdb,
     &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
     &  //".pdb")
      else
        cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
     &  //".cx"
      endif
#else
      if(mdpdb) then
         open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
      else
         cartname=prefix(:ilen(prefix))//"_MD.cx"
      endif
#endif
      if (usampl) then
        write (qstr,'(256(1h ))')
        ipos=1
        do i=1,nfrag
          iq = qinfrag(i)*10
          iw = wfrag(i)/100
          if (iw.gt.0) then
            write (iout,*) "Frag",qinfrag(i),wfrag(i),iq,iw
            write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
            ipos=ipos+7
          endif
        enddo
        do i=1,npair
          iq = qinpair(i)*10
          iw = wpair(i)/100
          if (iw.gt.0) then
            write (iout,*) "Pair",i,qinpair(i),wpair(i),iq,iw
            write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
            ipos=ipos+7
          endif
        enddo
        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
      endif
      icg=1
      if (rest) then
        write(iout,*) "Initial state will be read from file ",
     &    rest2name(:ilen(rest2name))
        call readrst
      else
c Generate initial velocities
        write(iout,*) "Initial velocities randomly generated"
        call random_vel
        totT=0.0d0
      endif
c      rest2name = prefix(:ilen(prefix))//'.rst'
      write(iout,*) "Initial backbone velocities"
      do i=nnt,nct-1
        write(iout,*) (d_t(j,i),j=1,3)
      enddo
      write(iout,*) "Initial side-chain velocities"
      do i=nnt,nct
        write(iout,*) (d_t(j,i+nres),j=1,3)
      enddo 		 	 
c  Zeroing the total angular momentum of the system
      write(iout,*) "Calling the zero-angular 
     & momentum subroutine"
      call inertia_tensor  
c  Getting the potential energy and forces and velocities and accelerations
      call vcm_vel(vcm)
c      write (iout,*) "velocity of the center of the mass:"
c      write (iout,*) (vcm(j),j=1,3)
      do j=1,3
        d_t(j,0)=d_t(j,0)-vcm(j)
      enddo
c Removing the velocity of the center of mass
      call vcm_vel(vcm)
      write (iout,*) "vcm right after adjustment:"
      write (iout,*) (vcm(j),j=1,3) 
      if (.not.rest) then	  	
         call chainbuild
      endif 	  
      call chainbuild_cart
      call kinetic(EK)
      if (tbf) then
        call verlet_bath(EK)
      endif	  
      kinetic_T=2.0d0/(dimen*Rb)*EK
      call cartprint
      call intout
      call zerograd
      call etotal(energia)
      potE=energia(0)
      call cartgrad
      call lagrangian
      call max_accel
      if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
      write(iout,*) "Potential energy"
      write(iout,*) (energia(i),i=0,n_ene)
      potE=energia(0)-energia(20)
      totE=EK+potE
      call statout(itime)
      write (iout,*) "Initial:",
     &   " Kinetic energy",EK," potential energy",potE, 
     &   " total energy",totE," maximum acceleration ",
     &   amax
      if (large) then
        write (iout,*) "Initial velocities"
        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
        write (iout,*) "Initial 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
      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
c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
      enddo 
      if (RESPA) then
        call zerograd
        call etotal_long(energia_long)
        call cartgrad
        call lagrangian
c        call etotal_short(energia_short)
c        write (iout,*) "energia_long",energia_long(0),
c     &   " energia_short",energia_short(0),
c     &   " total",energia_long(0)+energia_short(0)
      endif
      return
      end