X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fkinetic_lesyng.f;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fkinetic_lesyng.f;h=db959b3b3ecf1143ff08ee38a80bd3c1a4546e16;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git

diff --git a/source/unres/src_MD-M-newcorr/kinetic_lesyng.f b/source/unres/src_MD-M-newcorr/kinetic_lesyng.f
new file mode 100644
index 0000000..db959b3
--- /dev/null
+++ b/source/unres/src_MD-M-newcorr/kinetic_lesyng.f
@@ -0,0 +1,104 @@
+       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=iabs(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=iabs(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	
+	
+	
+	
+