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'
      character*16 form
      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.REMD'
      real*8 energia_long(0:n_ene),
     &  energia_short(0:n_ene),vcm(3),incr(3)
      double precision cm(3),L(3),xv,sigv,lowb,highb
      double precision varia(maxvar)
      character*256 qstr
      integer ilen
      external ilen
      character*50 tytul
      logical file_exist
      common /gucio/ cm
      d_time0=d_time
c      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(mdpdb) then
        if (ilen(tmpdir).gt.0) 
     &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
     &      liczba(:ilen(liczba))//".pdb")
        open(ipdb,
     &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
     &  //".pdb")
      else
#ifdef NOXDR
        if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) 
     &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
     &      liczba(:ilen(liczba))//".x")
        cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
     &  //".x"
#else
        if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) 
     &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
     &      liczba(:ilen(liczba))//".cx")
        cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
     &  //".cx"
#endif
      endif
#else
      if(mdpdb) then
         if (ilen(tmpdir).gt.0) 
     &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
         open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
      else
         if (ilen(tmpdir).gt.0) 
     &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
         cartname=prefix(:ilen(prefix))//"_MD.cx"
      endif
#endif
      if (usampl) then
        write (qstr,'(256(1h ))')
        ipos=1
        do i=1,nfrag
          iq = qinfrag(i,iset)*10
          iw = wfrag(i,iset)/100
          if (iw.gt.0) then
            if(me.eq.king.or..not.out1file)
     &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),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,iset)*10
          iw = wpair(i,iset)/100
          if (iw.gt.0) then
            if(me.eq.king.or..not.out1file)
     &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
            write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
            ipos=ipos+7
          endif
        enddo
c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
#ifdef NOXDR
c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
#else
c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
#endif
c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
      endif
      icg=1
      if (rest) then
       if (restart1file) then
         if (me.eq.king)
     &     inquire(file=mremd_rst_name,exist=file_exist)
           write (*,*) me," Before broadcast: file_exist",file_exist
         call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
     &          IERR)
         write (*,*) me," After broadcast: file_exist",file_exist
c        inquire(file=mremd_rst_name,exist=file_exist)
        if(me.eq.king.or..not.out1file)
     &   write(iout,*) "Initial state read by master and distributed"
       else
         if (ilen(tmpdir).gt.0)
     &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
     &      //liczba(:ilen(liczba))//'.rst')
        inquire(file=rest2name,exist=file_exist)
       endif
       if(file_exist) then
         if(.not.restart1file) then
           if(me.eq.king.or..not.out1file)
     &      write(iout,*) "Initial state will be read from file ",
     &      rest2name(:ilen(rest2name))
           call readrst
         endif  
         call rescale_weights(t_bath)
       else
        if(me.eq.king.or..not.out1file)then
         if (restart1file) then
          write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
     &       " does not exist"
         else
          write(iout,*) "File ",rest2name(:ilen(rest2name)),
     &       " does not exist"
         endif
         write(iout,*) "Initial velocities randomly generated"
        endif
        call random_vel
        totT=0.0d0
       endif
      else
c Generate initial velocities
        if(me.eq.king.or..not.out1file)
     &   write(iout,*) "Initial velocities randomly generated"
        call random_vel
        totT=0.0d0
      endif
c      rest2name = prefix(:ilen(prefix))//'.rst'
      if(me.eq.king.or..not.out1file)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 		 	 
c  Zeroing the total angular momentum of the system
       write(iout,*) "Calling the zero-angular 
     & momentum subroutine"
      endif
      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)
      if(me.eq.king.or..not.out1file)then
       write (iout,*) "vcm right after adjustment:"
       write (iout,*) (vcm(j),j=1,3) 
      endif
      if (.not.rest) then	  	
         call chainbuild
         if(iranconf.ne.0) then
          if (overlapsc) then 
           print *, 'Calling OVERLAP_SC'
           call overlap_sc(fail)
          endif 

          if (searchsc) then 
           call sc_move(2,nres-1,10,1d10,nft_sc,etot)
           print *,'SC_move',nft_sc,etot
           if(me.eq.king.or..not.out1file)
     &      write(iout,*) 'SC_move',nft_sc,etot
          endif 

          if(dccart)then
           print *, 'Calling MINIM_DC'
           call minim_dc(etot,iretcode,nfun)
          else
           call geom_to_var(nvar,varia)
           print *,'Calling MINIMIZE.'
           call minimize(etot,varia,iretcode,nfun)
           call var_to_geom(nvar,varia)
          endif
          if(me.eq.king.or..not.out1file)
     &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
         endif
      endif 	  
      call chainbuild_cart
      call kinetic(EK)
      if (tbf) then
        call verlet_bath(EK)
      endif	  
      kinetic_T=2.0d0/(dimen3*Rb)*EK
      if(me.eq.king.or..not.out1file)then
       call cartprint
       call intout
      endif
#ifdef MPI
      tt0=MPI_Wtime()
#else
      tt0=tcpu()
#endif
      call zerograd
      call etotal(potEcomp)
#ifdef TIMING_ENE
#ifdef MPI
      t_etotal=t_etotal+MPI_Wtime()-tt0
#else
      t_etotal=t_etotal+tcpu()-tt0
#endif
#endif
      potE=potEcomp(0)
      call cartgrad
      call lagrangian
      call max_accel
      if (amax*d_time .gt. dvmax) then
        d_time=d_time*dvmax/amax
        if(me.eq.king.or..not.out1file) write (iout,*) 
     &   "Time step reduced to",d_time,
     &   " because of too large initial acceleration."
      endif
      if(me.eq.king.or..not.out1file)then 
       write(iout,*) "Potential energy and its components"
       call enerprint(potEcomp)
c       write(iout,*) (potEcomp(i),i=0,n_ene)
      endif
      potE=potEcomp(0)-potEcomp(20)
      totE=EK+potE
      itime=0
      if (ntwe.ne.0) call statout(itime)
      if(me.eq.king.or..not.out1file)
     &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
     &   " Kinetic energy",EK," potential energy",potE, 
     &   " total energy",totE," maximum acceleration ",
     &   amax
      if (large) then
        write (iout,*) "Initial coordinates"
        do i=1,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
     &    (c(j,i+nres),j=1,3)
        enddo
        write (iout,*) "Initial dC"
        do i=0,nres
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
     &    (dc(j,i+nres),j=1,3)
        enddo
        write (iout,*) "Initial velocities"
        write (iout,"(13x,' backbone ',23x,' side chain')")
        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
c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
          write (iout,'(i3,3f15.10,3x,3f15.10)') 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
#ifdef MPI
      tt0 =MPI_Wtime()
#else
      tt0 = tcpu()
#endif
        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
        call lagrangian
        if(.not.out1file .and. large) then
          write (iout,*) "energia_long",energia_long(0),
     &     " energia_short",energia_short(0),
     &     " total",energia_long(0)+energia_short(0)
          write (iout,*) "Initial fast-force 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
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
#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
        if(.not.out1file .and. large) then
          write (iout,*) "energia_long",energia_long(0)
          write (iout,*) "Initial slow-force 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
#ifdef MPI
        t_enegrad=t_enegrad+MPI_Wtime()-tt0
#else
        t_enegrad=t_enegrad+tcpu()-tt0
#endif
      endif
      return
      end