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)
49 cm(j)=cm(j)+msc(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(iti)*(pr(2)*pr(2)+pr(3)*pr(3))
75 Im(1,2)=Im(1,2)-msc(iti)*pr(1)*pr(2)
76 Im(1,3)=Im(1,3)-msc(iti)*pr(1)*pr(3)
77 Im(2,3)=Im(2,3)-msc(iti)*pr(2)*pr(3)
78 Im(2,2)=Im(2,2)+msc(iti)*(pr(3)*pr(3)+pr(1)*pr(1))
79 Im(3,3)=Im(3,3)+msc(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) 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 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
135 c write (iout,*) "Eigenvalues & Eigenvectors"
136 c write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
139 c write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
141 c Constructing the diagonalized matrix
143 if (dabs(eigval(i)).gt.1.0d-15) then
144 Id(i,i)=1.0d0/eigval(i)
151 Imcp(i,j)=eigvec(j,i)
157 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
164 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
168 c Calculating the total rotational velocity of the molecule
171 vrot(i)=vrot(i)+pr2(i,j)*L(j)
174 c Resetting the velocities
176 call vecpr(vrot(1),dc(1,i),vp)
178 d_t(j,i)=d_t(j,i)-vp(j)
182 if(itype(i).ne.10) then
184 call vecpr(vrot(1),dc(1,inres),vp)
186 d_t(j,inres)=d_t(j,inres)-vp(j)
191 c write(iout,*) "The angular momentum after adjustment:"
192 c write(iout,*) (L(j),j=1,3)
195 c----------------------------------------------------------------------------
196 subroutine angmom(cm,L)
197 implicit real*8 (a-h,o-z)
199 include 'COMMON.CONTROL'
202 include 'COMMON.CHAIN'
203 include 'COMMON.DERIV'
205 include 'COMMON.LOCAL'
206 include 'COMMON.INTERACT'
207 include 'COMMON.IOUNITS'
208 include 'COMMON.NAMES'
210 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
213 c Calculate the angular momentum
222 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
225 v(j)=incr(j)+0.5d0*d_t(j,i)
228 incr(j)=incr(j)+d_t(j,i)
230 call vecpr(pr(1),v(1),vp)
238 call vecpr(pr(1),pp(1),vp)
250 pr(j)=c(j,inres)-cm(j)
252 if (itype(i).ne.10) then
254 v(j)=incr(j)+d_t(j,inres)
261 call vecpr(pr(1),v(1),vp)
262 c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
263 c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
265 L(j)=L(j)+msc(iti)*vp(j)
267 c write (iout,*) "L",(l(j),j=1,3)
268 if (itype(i).ne.10) then
270 v(j)=incr(j)+d_t(j,inres)
272 call vecpr(dc(1,inres),d_t(1,inres),vp)
274 L(j)=L(j)+Isc(iti)*vp(j)
278 incr(j)=incr(j)+d_t(j,i)
283 c------------------------------------------------------------------------------
284 subroutine vcm_vel(vcm)
285 implicit real*8 (a-h,o-z)
289 include 'COMMON.CHAIN'
290 include 'COMMON.DERIV'
292 include 'COMMON.LOCAL'
293 include 'COMMON.INTERACT'
294 include 'COMMON.IOUNITS'
295 double precision vcm(3),vv(3),summas,amas
305 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
310 if (itype(i).ne.10) then
312 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
316 vcm(j)=vcm(j)+amas*vv(j)
323 c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas