source/unres/src_MD/md-diff/mts/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 #ifndef LANG0
   9       include 'COMMON.LANGEVIN'
  10 #else
  11       include 'COMMON.LANGEVIN.lang0'
  12 #endif
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.DERIV'
  15       include 'COMMON.GEO'
  16       include 'COMMON.LOCAL'
  17       include 'COMMON.INTERACT'
  18       include 'COMMON.IOUNITS'
  19       include 'COMMON.NAMES'
  20       include 'COMMON.TIME1'
  21       double precision xv,sigv,lowb,highb
  22 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
  23 c First generate velocities in the eigenspace of the G matrix
  24 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
  25 c      call flush(iout)
  26       xv=0.0d0
  27       ii=0
  28       do i=1,dimen
  29         do k=1,3
  30           ii=ii+1
  31           sigv=dsqrt((Rb*t_bath)/geigen(i))
  32           lowb=-5*sigv
  33           highb=5*sigv
  34           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
  35 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
  36 c     &      " d_t_work_new",d_t_work_new(ii)
  37         enddo
  38       enddo
  39 c diagnostics
  40 c      Ek1=0.0d0
  41 c      ii=0
  42 c      do i=1,dimen
  43 c        do k=1,3
  44 c          ii=ii+1
  45 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
  46 c        enddo
  47 c      enddo
  48 c      write (iout,*) "Ek from eigenvectors",Ek1
  49 c end diagnostics
  50 c Transform velocities to UNRES coordinate space
  51       do k=0,2
  52         do i=1,dimen
  53           ind=(i-1)*3+k+1
  54           d_t_work(ind)=0.0d0
  55           do j=1,dimen
  56             d_t_work(ind)=d_t_work(ind)
  57      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
  58           enddo
  59 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
  60 c          call flush(iout)
  61         enddo
  62       enddo
  63 c Transfer to the d_t vector
  64       do j=1,3
  65         d_t(j,0)=d_t_work(j)
  66       enddo
  67       ind=3
  68       do i=nnt,nct-1
  69         do j=1,3
  70           ind=ind+1
  71           d_t(j,i)=d_t_work(ind)
  72         enddo
  73       enddo
  74       do i=nnt,nct
  75         if (itype(i).ne.10) then
  76           do j=1,3
  77             ind=ind+1
  78             d_t(j,i+nres)=d_t_work(ind)
  79           enddo
  80         endif
  81       enddo
  82 c      call kinetic(EK)
  83 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
  84 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
  85 c      call flush(iout)
  86       return
  87       end