--- /dev/null
+ subroutine kinetic(KE_total)
+c----------------------------------------------------------------
+c This subroutine calculates the total kinetic energy of the chain
+c-----------------------------------------------------------------
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.MD'
+ include 'COMMON.IOUNITS'
+ double precision KE_total
+
+ integer i,j,k
+ double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),
+ & mag1,mag2,v(3)
+
+ KEt_p=0.0d0
+ KEt_sc=0.0d0
+c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
+c The translational part for peptide virtual bonds
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct-1
+c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
+ do j=1,3
+ v(j)=incr(j)+0.5d0*d_t(j,i)
+ enddo
+ vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
+ KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+c write(iout,*) 'KEt_p', KEt_p
+c The translational part for the side chain virtual bond
+c Only now we can initialize incr with zeros. It must be equal
+c to the velocities of the first Calpha.
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct
+ iti=itype(i)
+ if (itype(i).eq.10) then
+ do j=1,3
+ v(j)=incr(j)
+ enddo
+ else
+ do j=1,3
+ v(j)=incr(j)+d_t(j,nres+i)
+ enddo
+ endif
+c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
+c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
+ KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+c goto 111
+c write(iout,*) 'KEt_sc', KEt_sc
+c The part due to stretching and rotation of the peptide groups
+ KEr_p=0.0D0
+ do i=nnt,nct-1
+c write (iout,*) "i",i
+c write (iout,*) "i",i," mag1",mag1," mag2",mag2
+ do j=1,3
+ incr(j)=d_t(j,i)
+ enddo
+c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
+ KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2)
+ & +incr(3)*incr(3))
+ enddo
+c goto 111
+c write(iout,*) 'KEr_p', KEr_p
+c The rotational part of the side chain virtual bond
+ KEr_sc=0.0D0
+ do i=nnt,nct
+ iti=itype(i)
+ if (itype(i).ne.10) then
+ do j=1,3
+ incr(j)=d_t(j,nres+i)
+ enddo
+c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
+ KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+
+ & incr(3)*incr(3))
+ endif
+ enddo
+c The total kinetic energy
+ 111 continue
+c write(iout,*) 'KEr_sc', KEr_sc
+ KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
+c write (iout,*) "KE_total",KE_total
+ return
+ end
+
+
+
+