X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD_DFA%2Fmoments.f;fp=source%2Funres%2Fsrc_MD_DFA%2Fmoments.f;h=0000000000000000000000000000000000000000;hb=664d495e70d14eed4e97f7b8efd2e107dee2fd4e;hp=5adbf21178c16d17fa6e2b949e40dc80cc0d28de;hpb=77276d5043cefaf512451c3ad9f669ed22b90d04;p=unres.git diff --git a/source/unres/src_MD_DFA/moments.f b/source/unres/src_MD_DFA/moments.f deleted file mode 100644 index 5adbf21..0000000 --- a/source/unres/src_MD_DFA/moments.f +++ /dev/null @@ -1,328 +0,0 @@ - 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