X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fmoments.f;h=983ce36bf713984a56e45e52f9dcc091260933a7;hb=d174d7fb88aee4a9122f02ebcda537b868723498;hp=007c0891a8769e8bcc5e13986ee94448e5b575df;hpb=478a9d9a1c99eb3f4bc4ca676ff3162bdd01d633;p=unres.git diff --git a/source/unres/src_MD-M/moments.f b/source/unres/src_MD-M/moments.f index 007c089..983ce36 100644 --- a/source/unres/src_MD-M/moments.f +++ b/source/unres/src_MD-M/moments.f @@ -42,11 +42,11 @@ c calculating the center of the mass of the protein enddo M_SC=0.0d0 do i=nnt,nct - iti=itype(i) - M_SC=M_SC+msc(iti) + iti=iabs(itype(i)) + M_SC=M_SC+msc(iabs(iti)) inres=i+nres do j=1,3 - cm(j)=cm(j)+msc(iti)*c(j,inres) + cm(j)=cm(j)+msc(iabs(iti))*c(j,inres) enddo enddo do j=1,3 @@ -66,17 +66,17 @@ c calculating the center of the mass of the protein enddo do i=nnt,nct - iti=itype(i) + 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(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)) + 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 @@ -96,8 +96,8 @@ c calculating the center of the mass of the protein do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then - iti=itype(i) + 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) @@ -179,7 +179,7 @@ c Resetting the velocities enddo enddo do i=nnt,nct - if(itype(i).ne.10 .and. itype(i).ne.21) then + 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 @@ -244,12 +244,12 @@ c Calculate the angular momentum incr(j)=d_t(j,0) enddo do i=nnt,nct - iti=itype(i) + 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.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo @@ -262,10 +262,10 @@ c Calculate the angular momentum 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) + 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 + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo @@ -305,9 +305,9 @@ c------------------------------------------------------------------------------ vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i)) enddo endif - amas=msc(itype(i)) + amas=msc(iabs(itype(i))) summas=summas+amas - if (itype(i).ne.10 .and. itype(i).ne.21) then + 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