subroutine inertia_tensor c Calculating the intertia tensor for the entire protein in order to c remove the perpendicular components of velocity matrix which cause c the molecule to rotate. implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision Im(3,3),Imcp(3,3),cm(3),pr(3),M_SC, & eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3), & vpp(3,0:MAXRES),vs_p(3),pr1(3,3), & pr2(3,3),pp(3),incr(3),v(3),mag,mag2 common /gucio/ cm integer iti,inres do i=1,3 do j=1,3 Im(i,j)=0.0d0 pr1(i,j)=0.0d0 pr2(i,j)=0.0d0 enddo L(i)=0.0d0 cm(i)=0.0d0 vrot(i)=0.0d0 enddo c calculating the center of the mass of the protein do i=nnt,nct-1 do j=1,3 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i) enddo enddo do j=1,3 cm(j)=mp*cm(j) enddo M_SC=0.0d0 do i=nnt,nct iti=iabs(itype(i)) M_SC=M_SC+msc(iabs(iti)) inres=i+nres do j=1,3 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres) enddo enddo do j=1,3 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp) enddo do i=nnt,nct-1 do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3)) Im(1,2)=Im(1,2)-mp*pr(1)*pr(2) Im(1,3)=Im(1,3)-mp*pr(1)*pr(3) Im(2,3)=Im(2,3)-mp*pr(2)*pr(3) Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct iti=itype(i) inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) enddo Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3)) Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2) Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3) Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3) Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct-1 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 enddo do i=nnt,nct if (itype(i).ne.10 .and. itype(i).ne.21) then iti=itype(i) inres=i+nres Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* & dc_norm(1,inres))*vbld(inres)*vbld(inres) Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) endif enddo call angmom(cm,L) c write(iout,*) "The angular momentum before adjustment:" c write(iout,*) (L(j),j=1,3) Im(2,1)=Im(1,2) Im(3,1)=Im(1,3) Im(3,2)=Im(2,3) c Copying the Im matrix for the djacob subroutine do i=1,3 do j=1,3 Imcp(i,j)=Im(i,j) Id(i,j)=0.0d0 enddo enddo c Finding the eigenvectors and eignvalues of the inertia tensor call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval) c write (iout,*) "Eigenvalues & Eigenvectors" c write (iout,'(5x,3f10.5)') (eigval(i),i=1,3) c write (iout,*) c do i=1,3 c write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3) c enddo c Constructing the diagonalized matrix do i=1,3 if (dabs(eigval(i)).gt.1.0d-15) then Id(i,i)=1.0d0/eigval(i) else Id(i,i)=0.0d0 endif enddo do i=1,3 do j=1,3 Imcp(i,j)=eigvec(j,i) enddo enddo do i=1,3 do j=1,3 do k=1,3 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j) enddo enddo enddo do i=1,3 do j=1,3 do k=1,3 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j) enddo enddo enddo c Calculating the total rotational velocity of the molecule do i=1,3 do j=1,3 vrot(i)=vrot(i)+pr2(i,j)*L(j) enddo enddo c Resetting the velocities do i=nnt,nct-1 call vecpr(vrot(1),dc(1,i),vp) do j=1,3 d_t(j,i)=d_t(j,i)-vp(j) enddo enddo do i=nnt,nct if(itype(i).ne.10 .and. itype(i).ne.21) then inres=i+nres call vecpr(vrot(1),dc(1,inres),vp) do j=1,3 d_t(j,inres)=d_t(j,inres)-vp(j) enddo endif enddo call angmom(cm,L) c write(iout,*) "The angular momentum after adjustment:" c write(iout,*) (L(j),j=1,3) return end c---------------------------------------------------------------------------- subroutine angmom(cm,L) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3), & pp(3) integer iti,inres c Calculate the angular momentum do j=1,3 L(j)=0.0d0 enddo do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct-1 do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo call vecpr(pr(1),v(1),vp) do j=1,3 L(j)=L(j)+mp*vp(j) enddo do j=1,3 pr(j)=0.5d0*dc(j,i) pp(j)=0.5d0*d_t(j,i) enddo call vecpr(pr(1),pp(1),vp) do j=1,3 L(j)=L(j)+Ip*vp(j) enddo enddo do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct iti=itype(i) inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) enddo if (itype(i).ne.10 .and. itype(i).ne.21) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo else do j=1,3 v(j)=incr(j) enddo endif call vecpr(pr(1),v(1),vp) c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3), c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3) do j=1,3 L(j)=L(j)+msc(iabs(iti))*vp(j) enddo c write (iout,*) "L",(l(j),j=1,3) if (itype(i).ne.10 .and. itype(i).ne.21) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo call vecpr(dc(1,inres),d_t(1,inres),vp) do j=1,3 L(j)=L(j)+Isc(iti)*vp(j) enddo endif do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo return end c------------------------------------------------------------------------------ subroutine vcm_vel(vcm) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' double precision vcm(3),vv(3),summas,amas do j=1,3 vcm(j)=0.0d0 vv(j)=d_t(j,0) enddo summas=0.0d0 do i=nnt,nct if (i.lt.nct) then summas=summas+mp do j=1,3 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i)) enddo endif amas=msc(itype(i)) summas=summas+amas if (itype(i).ne.10 .and. itype(i).ne.21) then do j=1,3 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres)) enddo else do j=1,3 vcm(j)=vcm(j)+amas*vv(j) enddo endif do j=1,3 vv(j)=vv(j)+d_t(j,i) enddo enddo c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas do j=1,3 vcm(j)=vcm(j)/summas enddo return end