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
11 include 'COMMON.CHAIN'
12 include 'COMMON.DERIV'
14 include 'COMMON.LOCAL'
15 include 'COMMON.INTERACT'
17 include 'COMMON.LAGRANGE.5diag'
18 include 'COMMON.IOUNITS'
19 double precision KE_total
21 double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),
28 c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
29 c The translational part for peptide virtual bonds
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
38 v(j)=incr(j)+0.5d0*d_t(j,i)
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))
45 incr(j)=incr(j)+d_t(j,i)
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.
57 if (itype(i).eq.10 .and. itype(i).ne.ntyp1) then
63 v(j)=incr(j)+d_t(j,nres+i)
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)
71 incr(j)=incr(j)+d_t(j,i)
75 c write(iout,*) 'KEt_sc', KEt_sc
76 c The part due to stretching and rotation of the peptide groups
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
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)
90 c write(iout,*) 'KEr_p', KEr_p
91 c The rotational part of the side chain virtual bond
94 if (itype(i).ne.10.and.itype(i).ne.ntyp1) then
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)+
103 c The total kinetic energy
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
112 subroutine kinetic(KE_total)
113 c----------------------------------------------------------------
114 c This subroutine calculates the total kinetic energy of the chain
115 c-----------------------------------------------------------------
119 include 'COMMON.CHAIN'
120 include 'COMMON.DERIV'
122 include 'COMMON.LOCAL'
123 include 'COMMON.INTERACT'
125 include 'COMMON.LAGRANGE'
126 include 'COMMON.IOUNITS'
127 double precision KE_total
129 double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),
134 c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
135 c The translational part for peptide virtual bonds
140 c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
142 v(j)=incr(j)+0.5d0*d_t(j,i)
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))
147 incr(j)=incr(j)+d_t(j,i)
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.
159 if (itype(i).eq.10) then
165 v(j)=incr(j)+d_t(j,nres+i)
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)
173 incr(j)=incr(j)+d_t(j,i)
177 c write(iout,*) 'KEt_sc', KEt_sc
178 c The part due to stretching and rotation of the peptide groups
181 c write (iout,*) "i",i
182 c write (iout,*) "i",i," mag1",mag1," mag2",mag2
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)
191 c write(iout,*) 'KEr_p', KEr_p
192 c The rotational part of the side chain virtual bond
196 if (itype(i).ne.10) then
198 incr(j)=d_t(j,nres+i)
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)+
205 c The total kinetic energy
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