source/unres/src-HCD-5D/kinetic_CASC.F

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