source/unres/src-HCD-5D/kinetic_lesyng.F.safe

   1        subroutine kinetic(KE_total)
   2 c----------------------------------------------------------------
   3 c   This subroutine calculates the total kinetic energy of the chain
   4 c-----------------------------------------------------------------
   5       implicit none
   6       include 'DIMENSIONS'
   7       include 'COMMON.VAR'
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.DERIV'
  10       include 'COMMON.GEO'
  11       include 'COMMON.LOCAL'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.MD'
  14 #ifdef FIVEDIAG
  15       include 'COMMON.LAGRANGE.5diag'
  16 #else
  17       include 'COMMON.LAGRANGE'
  18 #endif
  19       include 'COMMON.IOUNITS'
  20       double precision KE_total
  21
  22       integer i,j,k,iti
  23       double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),
  24      & mag1,mag2,v(3)
  25
  26       KEt_p=0.0d0
  27       KEt_sc=0.0d0
  28 c      write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
  29 c   The translational part for peptide virtual bonds
  30       do j=1,3
  31         incr(j)=d_t(j,0)
  32       enddo
  33       do i=nnt,nct-1
  34 c        write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
  35         do j=1,3
  36           v(j)=incr(j)+0.5d0*d_t(j,i)
  37         enddo
  38         vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
  39         KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
  40         do j=1,3
  41           incr(j)=incr(j)+d_t(j,i)
  42         enddo
  43       enddo
  44 c      write(iout,*) 'KEt_p', KEt_p
  45 c The translational part for the side chain virtual bond
  46 c Only now we can initialize incr with zeros. It must be equal
  47 c to the velocities of the first Calpha.
  48       do j=1,3
  49         incr(j)=d_t(j,0)
  50       enddo
  51       do i=nnt,nct
  52         iti=iabs(itype(i))
  53         if (itype(i).eq.10) then
  54           do j=1,3
  55             v(j)=incr(j)
  56           enddo
  57         else
  58           do j=1,3
  59             v(j)=incr(j)+d_t(j,nres+i)
  60           enddo
  61         endif
  62 c        write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
  63 c        write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
  64         KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
  65         vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
  66         do j=1,3
  67           incr(j)=incr(j)+d_t(j,i)
  68         enddo
  69       enddo
  70 c      goto 111
  71 c      write(iout,*) 'KEt_sc', KEt_sc
  72 c  The part due to stretching and rotation of the peptide groups
  73        KEr_p=0.0D0
  74        do i=nnt,nct-1
  75 c        write (iout,*) "i",i
  76 c        write (iout,*) "i",i," mag1",mag1," mag2",mag2
  77         do j=1,3
  78           incr(j)=d_t(j,i)
  79         enddo
  80 c        write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
  81           KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2)
  82      &    +incr(3)*incr(3))
  83        enddo
  84 c      goto 111
  85 c       write(iout,*) 'KEr_p', KEr_p
  86 c  The rotational part of the side chain virtual bond
  87        KEr_sc=0.0D0
  88        do i=nnt,nct
  89         iti=iabs(itype(i))
  90         if (itype(i).ne.10) then
  91         do j=1,3
  92           incr(j)=d_t(j,nres+i)
  93         enddo
  94 c        write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
  95         KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+
  96      &    incr(3)*incr(3))
  97         endif
  98        enddo
  99 c The total kinetic energy
 100   111  continue
 101 c       write(iout,*) 'KEr_sc', KEr_sc
 102        KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
 103 c       write (iout,*) "KE_total",KE_total
 104        return
 105        end
 106
 107
 108
 109