Merge branch 'adasko' into bartek with corrections
[unres.git] / source / unres / src_MD-M / moments.f
index 007c089..983ce36 100644 (file)
@@ -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