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

   1 #ifdef FIVEDIAG
   2        subroutine kinetic(KE_total)
   3 c----------------------------------------------------------------
   4 c   This subroutine calculates the total kinetic energy of the chain
   5 c-----------------------------------------------------------------
   6 c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups
   7 c   inside, implemented with five-diagonal inertia matrix
   8       implicit none
   9       include 'DIMENSIONS'
  10       include 'COMMON.VAR'
  11       include 'COMMON.CHAIN'
  12       include 'COMMON.DERIV'
  13       include 'COMMON.GEO'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.INTERACT'
  16       include 'COMMON.MD'
  17       include 'COMMON.LAGRANGE.5diag'
  18       include 'COMMON.IOUNITS'
  19       double precision KE_total
  20       integer i,j,k,iti
  21       double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),
  22      & mag1,mag2,v(3)
  23
  24       KEt_p=0.0d0
  25       KEt_sc=0.0d0
  26       KEr_p=0.0D0
  27       KEr_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 c Skip dummy peptide groups
  36         if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then
  37           do j=1,3
  38             v(j)=incr(j)+0.5d0*d_t(j,i)
  39           enddo
  40 c          write (iout,*) "Kinetic trp:",i,(v(j),j=1,3)
  41           vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
  42           KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
  43         endif
  44         do j=1,3
  45           incr(j)=incr(j)+d_t(j,i)
  46         enddo
  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 j=1,3
  53         incr(j)=d_t(j,0)
  54       enddo
  55       do i=nnt,nct
  56         iti=iabs(itype(i))
  57         if (itype(i).eq.10 .and. itype(i).ne.ntyp1) then
  58           do j=1,3
  59             v(j)=incr(j)
  60          enddo
  61         else
  62           do j=1,3
  63             v(j)=incr(j)+d_t(j,nres+i)
  64          enddo
  65         endif
  66 c        write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
  67 c        write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
  68         KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
  69         vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
  70         do j=1,3
  71           incr(j)=incr(j)+d_t(j,i)
  72         enddo
  73       enddo
  74 c      goto 111
  75 c      write(iout,*) 'KEt_sc', KEt_sc
  76 c  The part due to stretching and rotation of the peptide groups
  77        do i=nnt,nct-1
  78          if (itype(i).ne.ntyp1.and.itype(i+1).ne.ntyp1) then
  79 c        write (iout,*) "i",i
  80 c        write (iout,*) "i",i," mag1",mag1," mag2",mag2
  81          do j=1,3
  82            incr(j)=d_t(j,i)
  83          enddo
  84 c         write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
  85          KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2)
  86      &   +incr(3)*incr(3))
  87          endif
  88        enddo
  89 c      goto 111
  90 c       write(iout,*) 'KEr_p', KEr_p
  91 c  The rotational part of the side chain virtual bond
  92        do i=nnt,nct
  93         iti=iabs(itype(i))
  94         if (itype(i).ne.10.and.itype(i).ne.ntyp1) then
  95         do j=1,3
  96           incr(j)=d_t(j,nres+i)
  97         enddo
  98 c        write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
  99         KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+
 100      &     incr(3)*incr(3))
 101         endif
 102        enddo
 103 c The total kinetic energy
 104   111  continue
 105 c       write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
 106 c     &  ' KEr_sc', KEr_sc
 107        KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
 108 c       write (iout,*) "KE_total",KE_total
 109        return
 110        end
 111 #else
 112        subroutine kinetic(KE_total)
 113 c----------------------------------------------------------------
 114 c   This subroutine calculates the total kinetic energy of the chain
 115 c-----------------------------------------------------------------
 116       implicit none
 117       include 'DIMENSIONS'
 118       include 'COMMON.VAR'
 119       include 'COMMON.CHAIN'
 120       include 'COMMON.DERIV'
 121       include 'COMMON.GEO'
 122       include 'COMMON.LOCAL'
 123       include 'COMMON.INTERACT'
 124       include 'COMMON.MD'
 125       include 'COMMON.LAGRANGE'
 126       include 'COMMON.IOUNITS'
 127       double precision KE_total
 128       integer i,j,k,iti
 129       double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),
 130      & mag1,mag2,v(3)
 131
 132       KEt_p=0.0d0
 133       KEt_sc=0.0d0
 134 c      write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
 135 c   The translational part for peptide virtual bonds
 136       do j=1,3
 137         incr(j)=d_t(j,0)
 138       enddo
 139       do i=nnt,nct-1
 140 c        write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
 141         do j=1,3
 142           v(j)=incr(j)+0.5d0*d_t(j,i)
 143        enddo
 144         vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
 145         KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
 146         do j=1,3
 147           incr(j)=incr(j)+d_t(j,i)
 148         enddo
 149       enddo
 150 c      write(iout,*) 'KEt_p', KEt_p
 151 c The translational part for the side chain virtual bond
 152 c Only now we can initialize incr with zeros. It must be equal
 153 c to the velocities of the first Calpha.
 154       do j=1,3
 155         incr(j)=d_t(j,0)
 156       enddo
 157       do i=nnt,nct
 158         iti=iabs(itype(i))
 159         if (itype(i).eq.10) then
 160           do j=1,3
 161             v(j)=incr(j)
 162           enddo
 163         else
 164           do j=1,3
 165             v(j)=incr(j)+d_t(j,nres+i)
 166           enddo
 167         endif
 168 c        write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
 169 c        write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
 170         KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
 171         vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
 172         do j=1,3
 173           incr(j)=incr(j)+d_t(j,i)
 174         enddo
 175       enddo
 176 c      goto 111
 177 c      write(iout,*) 'KEt_sc', KEt_sc
 178 c  The part due to stretching and rotation of the peptide groups
 179        KEr_p=0.0D0
 180        do i=nnt,nct-1
 181 c        write (iout,*) "i",i
 182 c        write (iout,*) "i",i," mag1",mag1," mag2",mag2
 183          do j=1,3
 184            incr(j)=d_t(j,i)
 185          enddo
 186 c        write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
 187          KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2)
 188      &    +incr(3)*incr(3))
 189        enddo
 190 c      goto 111
 191 c       write(iout,*) 'KEr_p', KEr_p
 192 c  The rotational part of the side chain virtual bond
 193        KEr_sc=0.0D0
 194        do i=nnt,nct
 195         iti=iabs(itype(i))
 196         if (itype(i).ne.10) then
 197         do j=1,3
 198           incr(j)=d_t(j,nres+i)
 199         enddo
 200 c        write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
 201         KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+
 202      &    incr(3)*incr(3))
 203         endif
 204        enddo
 205 c The total kinetic energy
 206   111  continue
 207 c       write(iout,*) 'KEr_sc', KEr_sc
 208        KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
 209 c       write (iout,*) "KE_total",KE_total
 210        return
 211        end
 212 #endif