critical bug in set_matrices

[unres4.git] / source / unres / MD.f90

diff --git a/source/unres/MD.f90 b/source/unres/MD.f90
index 82ecbe8..1af6b18 100644 (file)
--- a/source/unres/MD.f90
+++ b/source/unres/MD.f90
@@ -652,7 +652,7 @@
 !      include 'COMMON.INTERACT'
 !      include 'COMMON.MD'
 !      include 'COMMON.IOUNITS'
-      real(kind=8) :: KE_total
+      real(kind=8) :: KE_total,mscab
                                                              
       integer :: i,j,k,iti,mnum
       real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
@@ -673,6 +673,7 @@
       enddo
       do i=nnt,nct-1
        mnum=molnum(i)
+       if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 !        write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3) 
         do j=1,3
           v(j)=incr(j)+0.5d0*d_t(j,i)
@@ -693,6 +694,11 @@
       do i=nnt,nct
          mnum=molnum(i)
         iti=iabs(itype(i,mnum))
+        if (mnum.eq.5) then
+         mscab=0.0d0
+        else
+         mscab=msc(iti,mnum)
+        endif
 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
           .or.(mnum.eq.5)) then
@@ -706,7 +712,7 @@
         endif
 !        write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) 
 !        write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3) 
-        KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))            
+        KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))            
         vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
         do j=1,3
           incr(j)=incr(j)+d_t(j,i)
@@ -977,7 +983,8 @@
           enddo
         enddo
         do i=nnt,nct
-          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+          mnum=molnum(i)
+          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
             do j=1,3
               ind=ind+1
               v_work(ind)=d_t(j,i+nres)
@@ -2209,7 +2216,8 @@
           enddo
         endif
 ! Side chains
-        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
+        molnum(i).ne.5) then
           do j=1,3 
             epdriftij= &
              dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
@@ -2911,6 +2919,7 @@
       do i=nnt,nct
 !        if (itype(i,1).eq.10) then
          mnum=molnum(i)
+          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
         iti=iabs(itype(i,mnum))
 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
@@ -3742,6 +3751,7 @@
         M_PEP=0.0d0
         do i=nnt,nct-1
           mnum=molnum(i)
+          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
           if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
           M_PEP=M_PEP+mp(mnum)
           do j=1,3
@@ -3754,6 +3764,7 @@
         M_SC=0.0d0                             
         do i=nnt,nct
            mnum=molnum(i)
+           if (mnum.eq.5) cycle
            iti=iabs(itype(i,mnum))              
           M_SC=M_SC+msc(iabs(iti),mnum)
            inres=i+nres
@@ -3767,6 +3778,7 @@
        
         do i=nnt,nct-1
            mnum=molnum(i)
+          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
           do j=1,3
             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
           enddo
@@ -3781,6 +3793,7 @@
        do i=nnt,nct    
            mnum=molnum(i)
            iti=iabs(itype(i,mnum))
+           if (mnum.eq.5) cycle
            inres=i+nres
            do j=1,3
              pr(j)=c(j,inres)-cm(j)        
@@ -3931,7 +3944,7 @@
 !       include 'COMMON.INTERACT'
 !       include 'COMMON.IOUNITS'
 !       include 'COMMON.NAMES'
-      
+      real(kind=8) :: mscab
       real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
       integer :: iti,inres,i,j,mnum
 !  Calculate the angular momentum
@@ -3943,6 +3956,7 @@
        enddo                  
        do i=nnt,nct-1
           mnum=molnum(i)
+          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
           do j=1,3
             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
           enddo
@@ -3972,6 +3986,11 @@
           mnum=molnum(i)
          iti=iabs(itype(i,mnum))
          inres=i+nres
+        if (mnum.eq.5) then
+         mscab=0.0d0
+        else
+         mscab=msc(iti,mnum)
+        endif
          do j=1,3
            pr(j)=c(j,inres)-cm(j)          
          enddo
@@ -3989,7 +4008,7 @@
 !         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 !     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
          do j=1,3
-            L(j)=L(j)+msc(iabs(iti),mnum)*vp(j)
+            L(j)=L(j)+mscab*vp(j)
          enddo
 !         write (iout,*) "L",(l(j),j=1,3)
          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
@@ -4033,13 +4052,18 @@
        summas=0.0d0
        do i=nnt,nct
          mnum=molnum(i)
+         if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
          if (i.lt.nct) then
            summas=summas+mp(mnum)
            do j=1,3
              vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
            enddo
          endif
+         if (mnum.ne.5) then 
          amas=msc(iabs(itype(i,mnum)),mnum)
+         else
+         amas=0.0d0
+         endif
          summas=summas+amas                     
          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
           .and.(mnum.ne.5)) then