source/unres/src_MD/src/kinetic_lesyng.f

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