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) cycle
187 call vecpr(vrot(1),dc(1,i),vp)
189 d_t(j,i)=d_t(j,i)-vp(j)
194 call vecpr(vrot(1),dc(1,i),vp)
196 d_t(j,i)=d_t(j,i)-vp(j)
201 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
203 call vecpr(vrot(1),dc(1,inres),vp)
205 d_t(j,inres)=d_t(j,inres)-vp(j)
210 c write(iout,*) "The angular momentum after adjustment:"
211 c write(iout,*) (L(j),j=1,3)
214 c----------------------------------------------------------------------------
215 subroutine angmom(cm,L)
218 include 'COMMON.CONTROL'
222 include 'COMMON.LAGRANGE.5diag'
224 include 'COMMON.LAGRANGE'
228 include 'COMMON.LANGEVIN.lang0.5diag'
230 include 'COMMON.LANGEVIN.lang0'
233 include 'COMMON.LANGEVIN'
235 include 'COMMON.CHAIN'
236 include 'COMMON.DERIV'
238 include 'COMMON.LOCAL'
239 include 'COMMON.INTERACT'
240 include 'COMMON.IOUNITS'
241 include 'COMMON.NAMES'
242 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
244 integer iti,inres,i,j
245 c Calculate the angular momentum
253 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
255 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
258 v(j)=incr(j)+0.5d0*d_t(j,i)
261 incr(j)=incr(j)+d_t(j,i)
263 call vecpr(pr(1),v(1),vp)
271 call vecpr(pr(1),pp(1),vp)
283 pr(j)=c(j,inres)-cm(j)
285 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
287 v(j)=incr(j)+d_t(j,inres)
294 call vecpr(pr(1),v(1),vp)
295 c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
296 c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
298 L(j)=L(j)+msc(iabs(iti))*vp(j)
300 c write (iout,*) "L",(l(j),j=1,3)
301 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
303 v(j)=incr(j)+d_t(j,inres)
305 call vecpr(dc(1,inres),d_t(1,inres),vp)
307 L(j)=L(j)+Isc(iti)*vp(j)
311 incr(j)=incr(j)+d_t(j,i)
316 c------------------------------------------------------------------------------
317 subroutine vcm_vel(vcm)
323 include 'COMMON.LAGRANGE.5diag'
325 include 'COMMON.LAGRANGE'
327 include 'COMMON.CHAIN'
328 include 'COMMON.DERIV'
330 include 'COMMON.LOCAL'
331 include 'COMMON.INTERACT'
332 include 'COMMON.IOUNITS'
333 double precision vcm(3),vv(3),summas,amas
344 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
347 amas=msc(iabs(itype(i)))
349 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
351 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
355 vcm(j)=vcm(j)+amas*vv(j)
362 c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
369 subroutine inertia_tensor
370 c Calculating the intertia tensor for the entire protein in order to
371 c remove the perpendicular components of velocity matrix which cause
372 c the molecule to rotate.
375 include 'COMMON.CONTROL'
379 include 'COMMON.LAGRANGE.5diag'
381 include 'COMMON.LAGRANGE'
383 include 'COMMON.CHAIN'
384 include 'COMMON.DERIV'
386 include 'COMMON.LOCAL'
387 include 'COMMON.INTERACT'
388 include 'COMMON.IOUNITS'
389 include 'COMMON.NAMES'
391 double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC,
392 & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),
393 & vpp(3,0:MAXRES),vs_p(3),pr1(3,3),
394 & pr2(3,3),pp(3),incr(3),v(3),mag,mag2
396 integer iti,inres,i,j,k
407 c calculating the center of the mass of the protein
410 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
419 M_SC=M_SC+msc(iabs(iti))
422 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
426 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
431 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
433 Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
434 Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
435 Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
436 Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
437 Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
438 Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
445 pr(j)=c(j,inres)-cm(j)
447 Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
448 Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
449 Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
450 Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
451 Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
452 Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
456 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))*
457 & vbld(i+1)*vbld(i+1)*0.25d0
458 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))*
459 & vbld(i+1)*vbld(i+1)*0.25d0
460 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))*
461 & vbld(i+1)*vbld(i+1)*0.25d0
462 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))*
463 & vbld(i+1)*vbld(i+1)*0.25d0
464 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))*
465 & vbld(i+1)*vbld(i+1)*0.25d0
466 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))*
467 & vbld(i+1)*vbld(i+1)*0.25d0
472 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
475 Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)*
476 & dc_norm(1,inres))*vbld(inres)*vbld(inres)
477 Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)*
478 & dc_norm(2,inres))*vbld(inres)*vbld(inres)
479 Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)*
480 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
481 Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)*
482 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
483 Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)*
484 & dc_norm(2,inres))*vbld(inres)*vbld(inres)
485 Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)*
486 & dc_norm(3,inres))*vbld(inres)*vbld(inres)
491 c write(iout,*) "The angular momentum before adjustment:"
492 c write(iout,*) (L(j),j=1,3)
498 c Copying the Im matrix for the djacob subroutine
506 c Finding the eigenvectors and eignvalues of the inertia tensor
507 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
508 c write (iout,*) "Eigenvalues & Eigenvectors"
509 c write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
512 c write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
514 c Constructing the diagonalized matrix
516 if (dabs(eigval(i)).gt.1.0d-15) then
517 Id(i,i)=1.0d0/eigval(i)
524 Imcp(i,j)=eigvec(j,i)
530 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
537 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
541 c Calculating the total rotational velocity of the molecule
544 vrot(i)=vrot(i)+pr2(i,j)*L(j)
547 c Resetting the velocities
549 call vecpr(vrot(1),dc(1,i),vp)
551 d_t(j,i)=d_t(j,i)-vp(j)
555 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
557 call vecpr(vrot(1),dc(1,inres),vp)
559 d_t(j,inres)=d_t(j,inres)-vp(j)
564 c write(iout,*) "The angular momentum after adjustment:"
565 c write(iout,*) (L(j),j=1,3)
568 c----------------------------------------------------------------------------
569 subroutine angmom(cm,L)
572 include 'COMMON.CONTROL'
576 include 'COMMON.LAGRANGE.5diag'
578 include 'COMMON.LAGRANGE'
582 include 'COMMON.LANGEVIN.lang0.5diag'
584 include 'COMMON.LANGEVIN.lang0'
587 include 'COMMON.LANGEVIN'
589 include 'COMMON.CHAIN'
590 include 'COMMON.DERIV'
592 include 'COMMON.LOCAL'
593 include 'COMMON.INTERACT'
594 include 'COMMON.IOUNITS'
595 include 'COMMON.NAMES'
597 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),
599 integer iti,inres,i,j
600 c Calculate the angular momentum
609 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
612 v(j)=incr(j)+0.5d0*d_t(j,i)
615 incr(j)=incr(j)+d_t(j,i)
617 call vecpr(pr(1),v(1),vp)
625 call vecpr(pr(1),pp(1),vp)
637 pr(j)=c(j,inres)-cm(j)
639 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
641 v(j)=incr(j)+d_t(j,inres)
648 call vecpr(pr(1),v(1),vp)
649 c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
650 c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
652 L(j)=L(j)+msc(iabs(iti))*vp(j)
654 c write (iout,*) "L",(l(j),j=1,3)
655 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
657 v(j)=incr(j)+d_t(j,inres)
659 call vecpr(dc(1,inres),d_t(1,inres),vp)
661 L(j)=L(j)+Isc(iti)*vp(j)
665 incr(j)=incr(j)+d_t(j,i)
670 c------------------------------------------------------------------------------
671 subroutine vcm_vel(vcm)
677 include 'COMMON.LAGRANGE.5diag'
679 include 'COMMON.LAGRANGE'
681 include 'COMMON.CHAIN'
682 include 'COMMON.DERIV'
684 include 'COMMON.LOCAL'
685 include 'COMMON.INTERACT'
686 include 'COMMON.IOUNITS'
687 double precision vcm(3),vv(3),summas,amas
698 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
701 amas=msc(iabs(itype(i)))
703 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
705 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
709 vcm(j)=vcm(j)+amas*vv(j)
716 c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas