1 subroutine inertia_tensor
2 c Calculating the intertia tensor for the entire protein in order to
3 c remove the perpendicular components of velocity matrix which cause
4 c the molecule to rotate.
5 implicit real*8 (a-h,o-z)
7 include 'COMMON.CONTROL'
10 include 'COMMON.CHAIN'
11 include 'COMMON.DERIV'
13 include 'COMMON.LOCAL'
14 include 'COMMON.INTERACT'
15 include 'COMMON.IOUNITS'
16 include 'COMMON.NAMES'
18 double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
19 & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
20 & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
21 & pr2(3,3),pp(3),incr(3),v(3),mag,mag2
34 c calculating the center of the mass of the protein
37 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
46 M_SC=M_SC+msc(iabs(iti))
49 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
53 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
58 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
60 Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
61 Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
62 Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
63 Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
64 Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
65 Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
72 pr(j)=c(j,inres)-cm(j)
74 Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
75 Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
76 Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
77 Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
78 Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
79 Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
83 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*
84 & vbld(i+1)*vbld(i+1)*0.25d0
85 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
86 & vbld(i+1)*vbld(i+1)*0.25d0
87 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
88 & vbld(i+1)*vbld(i+1)*0.25d0
89 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
90 & vbld(i+1)*vbld(i+1)*0.25d0
91 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
92 & vbld(i+1)*vbld(i+1)*0.25d0
93 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
94 & vbld(i+1)*vbld(i+1)*0.25d0
99 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
102 Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
103 & dc_norm(1,inres))*vbld(inres)*vbld(inres)
104 Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
105 & dc_norm(2,inres))*vbld(inres)*vbld(inres)
106 Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
107 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
108 Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
109 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
110 Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
111 & dc_norm(2,inres))*vbld(inres)*vbld(inres)
112 Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
113 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
118 c write(iout,*) "The angular momentum before adjustment:"
119 c write(iout,*) (L(j),j=1,3)
125 c Copying the Im matrix for the djacob subroutine
133 c Finding the eigenvectors and eignvalues of the inertia tensor
134 c write (iout,*) "Calling djacob"
136 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
137 c write (iout,*) "Eigenvalues & Eigenvectors"
138 c write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
141 c write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
144 c Constructing the diagonalized matrix
146 if (dabs(eigval(i)).gt.1.0d-15) then
147 Id(i,i)=1.0d0/eigval(i)
154 Imcp(i,j)=eigvec(j,i)
160 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
167 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
171 c Calculating the total rotational velocity of the molecule
174 vrot(i)=vrot(i)+pr2(i,j)*L(j)
177 c Resetting the velocities
179 call vecpr(vrot(1),dc(1,i),vp)
181 d_t(j,i)=d_t(j,i)-vp(j)
185 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
187 call vecpr(vrot(1),dc(1,inres),vp)
189 d_t(j,inres)=d_t(j,inres)-vp(j)
194 c write(iout,*) "The angular momentum after adjustment:"
195 c write(iout,*) (L(j),j=1,3)
198 c----------------------------------------------------------------------------
199 subroutine angmom(cm,L)
200 implicit real*8 (a-h,o-z)
202 include 'COMMON.CONTROL'
205 include 'COMMON.CHAIN'
206 include 'COMMON.DERIV'
208 include 'COMMON.LOCAL'
209 include 'COMMON.INTERACT'
210 include 'COMMON.IOUNITS'
211 include 'COMMON.NAMES'
213 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
216 c Calculate the angular momentum
225 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
228 v(j)=incr(j)+0.5d0*d_t(j,i)
231 incr(j)=incr(j)+d_t(j,i)
233 call vecpr(pr(1),v(1),vp)
241 call vecpr(pr(1),pp(1),vp)
253 pr(j)=c(j,inres)-cm(j)
255 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
257 v(j)=incr(j)+d_t(j,inres)
264 call vecpr(pr(1),v(1),vp)
265 c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
266 c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
268 L(j)=L(j)+msc(iabs(iti))*vp(j)
270 c write (iout,*) "L",(l(j),j=1,3)
271 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
273 v(j)=incr(j)+d_t(j,inres)
275 call vecpr(dc(1,inres),d_t(1,inres),vp)
277 L(j)=L(j)+Isc(iti)*vp(j)
281 incr(j)=incr(j)+d_t(j,i)
286 c------------------------------------------------------------------------------
287 subroutine vcm_vel(vcm)
288 implicit real*8 (a-h,o-z)
292 include 'COMMON.CHAIN'
293 include 'COMMON.DERIV'
295 include 'COMMON.LOCAL'
296 include 'COMMON.INTERACT'
297 include 'COMMON.IOUNITS'
298 double precision vcm(3),vv(3),summas,amas
308 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
311 amas=msc(iabs(itype(i)))
313 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
315 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
319 vcm(j)=vcm(j)+amas*vv(j)
326 c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas