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