c-----------------------------------------------------------
      subroutine sd_verlet_p_setup
c Sets up the parameters of stochastic Verlet algorithm       
      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 emgdt(MAXRES6),
     & pterm,vterm,rho,rhoc,vsig,
     & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
     & afric_vec(MAXRES6),prand_vec(MAXRES6),
     & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
      logical lprn /.false./
      double precision zero /1.0d-8/, gdt_radius /0.05d0/ 
      double precision ktm
      tt0 = tcpu()
c
c AL 8/17/04 Code adapted from tinker
c
c Get the frictional and random terms for stochastic dynamics in the
c eigenspace of mass-scaled UNRES friction matrix
c
      do i = 1, dimen
            gdt = fricgam(i) * d_time
c
c Stochastic dynamics reduces to simple MD for zero friction
c
            if (gdt .le. zero) then
               pfric_vec(i) = 1.0d0
               vfric_vec(i) = d_time
               afric_vec(i) = 0.5d0 * d_time * d_time
               prand_vec(i) = 0.0d0
               vrand_vec1(i) = 0.0d0
               vrand_vec2(i) = 0.0d0
c
c Analytical expressions when friction coefficient is large
c
            else 
               if (gdt .ge. gdt_radius) then
                  egdt = dexp(-gdt)
                  pfric_vec(i) = egdt
                  vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
                  afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
                  pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
                  vterm = 1.0d0 - egdt**2
                  rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
c
c Use series expansions when friction coefficient is small
c
               else
                  gdt2 = gdt * gdt
                  gdt3 = gdt * gdt2
                  gdt4 = gdt2 * gdt2
                  gdt5 = gdt2 * gdt3
                  gdt6 = gdt3 * gdt3
                  gdt7 = gdt3 * gdt4
                  gdt8 = gdt4 * gdt4
                  gdt9 = gdt4 * gdt5
                  afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
     &                          - gdt5/120.0d0 + gdt6/720.0d0
     &                          - gdt7/5040.0d0 + gdt8/40320.0d0
     &                          - gdt9/362880.0d0) / fricgam(i)**2
                  vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
                  pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
                  pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
     &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
     &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
     &                       + 127.0d0*gdt9/90720.0d0
                  vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
     &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
     &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
     &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
                  rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
     &                       - 17.0d0*gdt2/1280.0d0
     &                       + 17.0d0*gdt3/6144.0d0
     &                       + 40967.0d0*gdt4/34406400.0d0
     &                       - 57203.0d0*gdt5/275251200.0d0
     &                       - 1429487.0d0*gdt6/13212057600.0d0)
               end if
c
c Compute the scaling factors of random terms for the nonzero friction case
c
               ktm = 0.5d0*d_time/fricgam(i)
               psig = dsqrt(ktm*pterm) / fricgam(i)
               vsig = dsqrt(ktm*vterm)
               rhoc = dsqrt(1.0d0 - rho*rho)
               prand_vec(i) = psig 
               vrand_vec1(i) = vsig * rho 
               vrand_vec2(i) = vsig * rhoc
            end if
      end do
      if (lprn) then
      write (iout,*) 
     &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
     &  " vrand_vec2"
      do i=1,dimen
        write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
     &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
      enddo
      endif
c
c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
c
      call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
      call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
      call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
      call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
      call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
      call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
c      call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
c      call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
c      call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
c      call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
c      call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
c      call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
      t_sdsetup=t_sdsetup+tcpu()-tt0
      return
      end