X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fmoments.f;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fmoments.f;h=983ce36bf713984a56e45e52f9dcc091260933a7;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/moments.f b/source/unres/src_MD-M-newcorr/moments.f new file mode 100644 index 0000000..983ce36 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/moments.f @@ -0,0 +1,328 @@ + 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=iabs(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.ntyp1) then + iti=iabs(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.ntyp1) 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=iabs(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.ntyp1) 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.ntyp1) 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(iabs(itype(i))) + summas=summas+amas + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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