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