+++ /dev/null
- 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=itype(i)
- M_SC=M_SC+msc(iti)
- inres=i+nres
- do j=1,3
- cm(j)=cm(j)+msc(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(iti)*(pr(2)*pr(2)+pr(3)*pr(3))
- Im(1,2)=Im(1,2)-msc(iti)*pr(1)*pr(2)
- Im(1,3)=Im(1,3)-msc(iti)*pr(1)*pr(3)
- Im(2,3)=Im(2,3)-msc(iti)*pr(2)*pr(3)
- Im(2,2)=Im(2,2)+msc(iti)*(pr(3)*pr(3)+pr(1)*pr(1))
- Im(3,3)=Im(3,3)+msc(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) 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) 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) 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(iti)*vp(j)
- enddo
-c write (iout,*) "L",(l(j),j=1,3)
- if (itype(i).ne.10) 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) 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
projects / unres.git / blobdiff