--- /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