source/unres/src_MD/md-diff/np/random_vel.f

   1 c-----------------------------------------------------------
   2       subroutine random_vel
   3       implicit real*8 (a-h,o-z)
   4       include 'DIMENSIONS'
   5       include 'COMMON.CONTROL'
   6       include 'COMMON.VAR'
   7       include 'COMMON.MD'
   8       include 'COMMON.LANGEVIN'
   9       include 'COMMON.CHAIN'
  10       include 'COMMON.DERIV'
  11       include 'COMMON.GEO'
  12       include 'COMMON.LOCAL'
  13       include 'COMMON.INTERACT'
  14       include 'COMMON.IOUNITS'
  15       include 'COMMON.NAMES'
  16       include 'COMMON.TIME1'
  17       double precision xv,sigv,lowb,highb
  18 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
  19 c First generate velocities in the eigenspace of the G matrix
  20 c      write (iout,*) "Calling random_vel"
  21       xv=0.0d0
  22       do i=1,dimen
  23         sigv=dsqrt((Rb*t_bath)/geigen(i))
  24         lowb=-5*sigv
  25         highb=5*sigv
  26         d_t_work_new(i)=anorm_distr(xv,sigv,lowb,highb)
  27       enddo
  28 c      Ek1=0.0d0
  29 c      do i=1,dimen
  30 c        Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(i)**2
  31 c      enddo
  32 c Transform velocities to UNRES coordinate space
  33       do i=1,dimen
  34         d_t_work(i)=0.0d0
  35         do j=1,dimen
  36           d_t_work(i)=d_t_work(i)+Gvec(i,j)*d_t_work_new(j)
  37         enddo
  38       enddo
  39 c Transfer to the d_t vector
  40       do j=1,3
  41         d_t(j,0)=d_t_work(j)
  42       enddo
  43       ind=3
  44       do i=nnt,nct-1
  45         do j=1,3
  46           ind=ind+1
  47           d_t(j,i)=d_t_work(ind)
  48         enddo
  49       enddo
  50       do i=nnt,nct
  51         if (itype(i).ne.10) then
  52           do j=1,3
  53             ind=ind+1
  54             d_t(j,i+nres)=d_t_work(ind)
  55           enddo
  56         endif
  57       enddo
  58 c      call kinetic(EK)
  59 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
  60 c     &  2.0d0/(dimen*Rb)*EK,2.0d0/(dimen*Rb)*EK1
  61       return
  62       end