2 subroutine inertia_tensor
3 c Calculating the intertia tensor for the entire protein in order to
4 c remove the perpendicular components of velocity matrix which cause
5 c the molecule to rotate.
8 include 'COMMON.CONTROL'
12 include 'COMMON.LAGRANGE.5diag'
14 include 'COMMON.LAGRANGE'
16 include 'COMMON.CHAIN'
17 include 'COMMON.DERIV'
19 include 'COMMON.LOCAL'
20 include 'COMMON.INTERACT'
21 include 'COMMON.IOUNITS'
22 include 'COMMON.NAMES'
24 double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
25 & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
26 & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
27 & pr2(3,3),pp(3),incr(3),v(3),mag,mag2
29 integer iti,inres,i,j,k
40 c caulating the center of the mass of the protein
42 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
44 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
53 if (iti.eq.ntyp1) cycle
54 M_SC=M_SC+msc(iabs(iti))
57 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
61 cm(j)=cm(j)/(M_SC+dimenp*mp)
65 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
67 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
69 Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
70 Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
71 Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
72 Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
73 Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
74 Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
79 if (iti.eq.ntyp1) cycle
82 pr(j)=c(j,inres)-cm(j)
84 Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
85 Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
86 Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
87 Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
88 Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
89 Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
93 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
94 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*
95 & vbld(i+1)*vbld(i+1)*0.25d0
96 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
97 & vbld(i+1)*vbld(i+1)*0.25d0
98 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
99 & vbld(i+1)*vbld(i+1)*0.25d0
100 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
101 & vbld(i+1)*vbld(i+1)*0.25d0
102 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
103 & vbld(i+1)*vbld(i+1)*0.25d0
104 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
105 & vbld(i+1)*vbld(i+1)*0.25d0
109 if (iti.ne.10 .and. iti.ne.ntyp1) then
111 Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
112 & dc_norm(1,inres))*vbld(inres)*vbld(inres)
113 Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
114 & dc_norm(2,inres))*vbld(inres)*vbld(inres)
115 Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
116 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
117 Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
118 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
119 Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
120 & dc_norm(2,inres))*vbld(inres)*vbld(inres)
121 Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
122 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
127 c write(iout,*) "The angular momentum before adjustment:"
128 c write(iout,*) (L(j),j=1,3)
134 c Copng the Im matrix for the djacob subroutine
141 c Finding the eigenvectors and eignvalues of the inertia tensor
142 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
143 c write (iout,*) "Eigenvalues & Eigenvectors"
144 c write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
147 c write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
149 c Constructing the diagonalized matrix
151 if (dabs(eigval(i)).gt.1.0d-15) then
152 Id(i,i)=1.0d0/eigval(i)
159 Imcp(i,j)=eigvec(j,i)
165 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
172 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
176 c Calculating the total rotational velocity of the molecule
179 vrot(i)=vrot(i)+pr2(i,j)*L(j)
182 c Resetting the velocities
185 write (iout,*) itype(i+1),itype(i)
186 if (itype(i+1).ne.ntyp1 .and. itype(i).eq.ntyp1 .or.
187 & itype(i).ne.ntyp1 .and. itype(i+1).eq.ntyp1) cycle
188 call vecpr(vrot(1),dc(1,i),vp)
190 d_t(j,i)=d_t(j,i)-vp(j)
195 call vecpr(vrot(1),dc(1,i),vp)
197 d_t(j,i)=d_t(j,i)-vp(j)
202 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
204 call vecpr(vrot(1),dc(1,inres),vp)
206 d_t(j,inres)=d_t(j,inres)-vp(j)
211 c write(iout,*) "The angular momentum after adjustment:"
212 c write(iout,*) (L(j),j=1,3)
215 c----------------------------------------------------------------------------
216 subroutine angmom(cm,L)
219 include 'COMMON.CONTROL'
223 include 'COMMON.LAGRANGE.5diag'
225 include 'COMMON.LAGRANGE'
229 include 'COMMON.LANGEVIN.lang0.5diag'
231 include 'COMMON.LANGEVIN.lang0'
234 include 'COMMON.LANGEVIN'
236 include 'COMMON.CHAIN'
237 include 'COMMON.DERIV'
239 include 'COMMON.LOCAL'
240 include 'COMMON.INTERACT'
241 include 'COMMON.IOUNITS'
242 include 'COMMON.NAMES'
243 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
245 integer iti,inres,i,j
246 c Calculate the angular momentum
254 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
256 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
259 v(j)=incr(j)+0.5d0*d_t(j,i)
262 incr(j)=incr(j)+d_t(j,i)
264 call vecpr(pr(1),v(1),vp)
272 call vecpr(pr(1),pp(1),vp)
284 pr(j)=c(j,inres)-cm(j)
286 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
288 v(j)=incr(j)+d_t(j,inres)
295 call vecpr(pr(1),v(1),vp)
296 c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
297 c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
299 L(j)=L(j)+msc(iabs(iti))*vp(j)
301 c write (iout,*) "L",(l(j),j=1,3)
302 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
304 v(j)=incr(j)+d_t(j,inres)
306 call vecpr(dc(1,inres),d_t(1,inres),vp)
308 L(j)=L(j)+Isc(iti)*vp(j)
312 incr(j)=incr(j)+d_t(j,i)
317 c------------------------------------------------------------------------------
318 subroutine vcm_vel(vcm)
324 include 'COMMON.LAGRANGE.5diag'
326 include 'COMMON.LAGRANGE'
328 include 'COMMON.CHAIN'
329 include 'COMMON.DERIV'
331 include 'COMMON.LOCAL'
332 include 'COMMON.INTERACT'
333 include 'COMMON.IOUNITS'
334 double precision vcm(3),vv(3),summas,amas
345 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
348 amas=msc(iabs(itype(i)))
350 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
352 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
356 vcm(j)=vcm(j)+amas*vv(j)
363 c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
370 subroutine inertia_tensor
371 c Calculating the intertia tensor for the entire protein in order to
372 c remove the perpendicular components of velocity matrix which cause
373 c the molecule to rotate.
376 include 'COMMON.CONTROL'
380 include 'COMMON.LAGRANGE.5diag'
382 include 'COMMON.LAGRANGE'
384 include 'COMMON.CHAIN'
385 include 'COMMON.DERIV'
387 include 'COMMON.LOCAL'
388 include 'COMMON.INTERACT'
389 include 'COMMON.IOUNITS'
390 include 'COMMON.NAMES'
392 double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
393 & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
394 & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
395 & pr2(3,3),pp(3),incr(3),v(3),mag,mag2
397 integer iti,inres,i,j,k
408 c calculating the center of the mass of the protein
411 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
420 M_SC=M_SC+msc(iabs(iti))
423 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
427 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
432 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
434 Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
435 Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
436 Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
437 Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
438 Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
439 Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
446 pr(j)=c(j,inres)-cm(j)
448 Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
449 Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
450 Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
451 Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
452 Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
453 Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
457 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*
458 & vbld(i+1)*vbld(i+1)*0.25d0
459 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
460 & vbld(i+1)*vbld(i+1)*0.25d0
461 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
462 & vbld(i+1)*vbld(i+1)*0.25d0
463 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
464 & vbld(i+1)*vbld(i+1)*0.25d0
465 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
466 & vbld(i+1)*vbld(i+1)*0.25d0
467 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
468 & vbld(i+1)*vbld(i+1)*0.25d0
473 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
476 Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
477 & dc_norm(1,inres))*vbld(inres)*vbld(inres)
478 Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
479 & dc_norm(2,inres))*vbld(inres)*vbld(inres)
480 Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
481 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
482 Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
483 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
484 Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
485 & dc_norm(2,inres))*vbld(inres)*vbld(inres)
486 Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
487 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
492 c write(iout,*) "The angular momentum before adjustment:"
493 c write(iout,*) (L(j),j=1,3)
499 c Copying the Im matrix for the djacob subroutine
507 c Finding the eigenvectors and eignvalues of the inertia tensor
508 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
509 c write (iout,*) "Eigenvalues & Eigenvectors"
510 c write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
513 c write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
515 c Constructing the diagonalized matrix
517 if (dabs(eigval(i)).gt.1.0d-15) then
518 Id(i,i)=1.0d0/eigval(i)
525 Imcp(i,j)=eigvec(j,i)
531 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
538 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
542 c Calculating the total rotational velocity of the molecule
545 vrot(i)=vrot(i)+pr2(i,j)*L(j)
548 c Resetting the velocities
550 call vecpr(vrot(1),dc(1,i),vp)
552 d_t(j,i)=d_t(j,i)-vp(j)
556 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
558 call vecpr(vrot(1),dc(1,inres),vp)
560 d_t(j,inres)=d_t(j,inres)-vp(j)
565 c write(iout,*) "The angular momentum after adjustment:"
566 c write(iout,*) (L(j),j=1,3)
569 c----------------------------------------------------------------------------
570 subroutine angmom(cm,L)
573 include 'COMMON.CONTROL'
577 include 'COMMON.LAGRANGE.5diag'
579 include 'COMMON.LAGRANGE'
583 include 'COMMON.LANGEVIN.lang0.5diag'
585 include 'COMMON.LANGEVIN.lang0'
588 include 'COMMON.LANGEVIN'
590 include 'COMMON.CHAIN'
591 include 'COMMON.DERIV'
593 include 'COMMON.LOCAL'
594 include 'COMMON.INTERACT'
595 include 'COMMON.IOUNITS'
596 include 'COMMON.NAMES'
598 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
600 integer iti,inres,i,j
601 c Calculate the angular momentum
610 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
613 v(j)=incr(j)+0.5d0*d_t(j,i)
616 incr(j)=incr(j)+d_t(j,i)
618 call vecpr(pr(1),v(1),vp)
626 call vecpr(pr(1),pp(1),vp)
638 pr(j)=c(j,inres)-cm(j)
640 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
642 v(j)=incr(j)+d_t(j,inres)
649 call vecpr(pr(1),v(1),vp)
650 c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
651 c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
653 L(j)=L(j)+msc(iabs(iti))*vp(j)
655 c write (iout,*) "L",(l(j),j=1,3)
656 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
658 v(j)=incr(j)+d_t(j,inres)
660 call vecpr(dc(1,inres),d_t(1,inres),vp)
662 L(j)=L(j)+Isc(iti)*vp(j)
666 incr(j)=incr(j)+d_t(j,i)
671 c------------------------------------------------------------------------------
672 subroutine vcm_vel(vcm)
678 include 'COMMON.LAGRANGE.5diag'
680 include 'COMMON.LAGRANGE'
682 include 'COMMON.CHAIN'
683 include 'COMMON.DERIV'
685 include 'COMMON.LOCAL'
686 include 'COMMON.INTERACT'
687 include 'COMMON.IOUNITS'
688 double precision vcm(3),vv(3),summas,amas
699 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
702 amas=msc(iabs(itype(i)))
704 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
706 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
710 vcm(j)=vcm(j)+amas*vv(j)
717 c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas