adding ebend_nucl to UCGM+some further reading
[unres4.git] / source / unres / MD.f90
index 3d3dc39..1cc216e 100644 (file)
         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)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))         
+        KEt_sc=KEt_sc+msc(iti,1)*(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)
          incr(j)=d_t(j,nres+i)
        enddo
 !        write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) 
-       KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
+       KEr_sc=KEr_sc+Isc(iti,1)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
          incr(3)*incr(3))
         endif
        enddo
 ! The total kinetic energy     
   111  continue
 !       write(iout,*) 'KEr_sc', KEr_sc 
-       KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)         
+       KE_total=0.5d0*(mp(1)*KEt_p+KEt_sc+0.25d0*Ip(1)*KEr_p+KEr_sc)           
 !       write (iout,*) "KE_total",KE_total 
       return
       end subroutine kinetic
         if (i.lt.nct) then
           do k=1,3
 !            write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
-            Ek1=Ek1+0.5d0*mp*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
-            0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
+            Ek1=Ek1+0.5d0*mp(1)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
+            0.5d0*0.25d0*IP(1)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
           enddo
         endif
         if (itype(i,1).ne.10) ii=ii+3
-        write (iout,*) "i",i," itype",itype(i,1)," mass",msc(itype(i,1))
+        write (iout,*) "i",i," itype",itype(i,1)," mass",msc(itype(i,1),1)
         write (iout,*) "ii",ii
         do k=1,3
           ii=ii+1
           write (iout,*) "k",k," ii",ii,"EK1",EK1
-          if (iabs(itype(i,1)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i,1)))*(d_t_work(ii)-d_t_work(ii-3))**2
-          Ek1=Ek1+0.5d0*msc(iabs(itype(i,1)))*d_t_work(ii)**2
+          if (iabs(itype(i,1)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i,1),1))*(d_t_work(ii)-d_t_work(ii-3))**2
+          Ek1=Ek1+0.5d0*msc(iabs(itype(i,1)),1)*d_t_work(ii)**2
         enddo
         write (iout,*) "i",i," ii",ii
       enddo
           enddo
         enddo
         do j=1,3
-         cm(j)=mp*cm(j)
+         cm(j)=mp(1)*cm(j)
         enddo
         M_SC=0.0d0                             
         do i=nnt,nct
            iti=iabs(itype(i,1))                 
-          M_SC=M_SC+msc(iabs(iti))
+          M_SC=M_SC+msc(iabs(iti),1)
            inres=i+nres
            do j=1,3
-            cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)          
+            cm(j)=cm(j)+msc(iabs(iti),1)*c(j,inres)        
            enddo
         enddo
         do j=1,3
-          cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
+          cm(j)=cm(j)/(M_SC+(nct-nnt)*mp(1))
         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))
+          Im(1,1)=Im(1,1)+mp(1)*(pr(2)*pr(2)+pr(3)*pr(3))
+          Im(1,2)=Im(1,2)-mp(1)*pr(1)*pr(2)
+          Im(1,3)=Im(1,3)-mp(1)*pr(1)*pr(3)
+          Im(2,3)=Im(2,3)-mp(1)*pr(2)*pr(3)    
+          Im(2,2)=Im(2,2)+mp(1)*(pr(3)*pr(3)+pr(1)*pr(1))
+          Im(3,3)=Im(3,3)+mp(1)*(pr(1)*pr(1)+pr(2)*pr(2))
         enddo                  
         
        do i=nnt,nct    
            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))                
+          Im(1,1)=Im(1,1)+msc(iabs(iti),1)*(pr(2)*pr(2)+pr(3)*pr(3))
+          Im(1,2)=Im(1,2)-msc(iabs(iti),1)*pr(1)*pr(2)
+          Im(1,3)=Im(1,3)-msc(iabs(iti),1)*pr(1)*pr(3)
+          Im(2,3)=Im(2,3)-msc(iabs(iti),1)*pr(2)*pr(3) 
+          Im(2,2)=Im(2,2)+msc(iabs(iti),1)*(pr(3)*pr(3)+pr(1)*pr(1))
+          Im(3,3)=Im(3,3)+msc(iabs(iti),1)*(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))* &    
+          Im(1,1)=Im(1,1)+Ip(1)*(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))* &
+         Im(1,2)=Im(1,2)+Ip(1)*(-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))* &
+          Im(1,3)=Im(1,3)+Ip(1)*(-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))* &
+          Im(2,3)=Im(2,3)+Ip(1)*(-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))* &
+          Im(2,2)=Im(2,2)+Ip(1)*(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))* &
+          Im(3,3)=Im(3,3)+Ip(1)*(1-dc_norm(3,i)*dc_norm(3,i))* &
           vbld(i+1)*vbld(i+1)*0.25d0
         enddo
         
          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
            iti=iabs(itype(i,1))                 
            inres=i+nres
-          Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* &
+          Im(1,1)=Im(1,1)+Isc(iti,1)*(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)* &
+          Im(1,2)=Im(1,2)-Isc(iti,1)*(dc_norm(1,inres)* &
          dc_norm(2,inres))*vbld(inres)*vbld(inres)
-          Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* &
+          Im(1,3)=Im(1,3)-Isc(iti,1)*(dc_norm(1,inres)* &
          dc_norm(3,inres))*vbld(inres)*vbld(inres)
-          Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* &
+          Im(2,3)=Im(2,3)-Isc(iti,1)*(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)* &
+          Im(2,2)=Im(2,2)+Isc(iti,1)*(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)* &
+          Im(3,3)=Im(3,3)+Isc(iti,1)*(1-dc_norm(3,inres)* &
                  dc_norm(3,inres))*vbld(inres)*vbld(inres)
          endif
         enddo
           enddo                
           call vecpr(pr(1),v(1),vp)
           do j=1,3
-            L(j)=L(j)+mp*vp(j)
+            L(j)=L(j)+mp(1)*vp(j)
           enddo
           do j=1,3
              pr(j)=0.5d0*dc(j,i)
           enddo
          call vecpr(pr(1),pp(1),vp)
          do j=1,3               
-             L(j)=L(j)+Ip*vp(j)         
+             L(j)=L(j)+Ip(1)*vp(j)      
           enddo
         enddo
         do j=1,3
 !         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))*vp(j)
+            L(j)=L(j)+msc(iabs(iti),1)*vp(j)
          enddo
 !         write (iout,*) "L",(l(j),j=1,3)
          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
            enddo
            call vecpr(dc(1,inres),d_t(1,inres),vp)
            do j=1,3                               
-             L(j)=L(j)+Isc(iti)*vp(j)   
+             L(j)=L(j)+Isc(iti,1)*vp(j)         
           enddo                           
          endif
         do j=1,3
        summas=0.0d0
        do i=nnt,nct
          if (i.lt.nct) then
-           summas=summas+mp
+           summas=summas+mp(1)
            do j=1,3
-             vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
+             vcm(j)=vcm(j)+mp(1)*(vv(j)+0.5d0*d_t(j,i))
            enddo
          endif
-         amas=msc(iabs(itype(i,1)))
+         amas=msc(iabs(itype(i,1)),1)
          summas=summas+amas                     
          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
            do j=1,3
       enddo
 !  Load peptide group radii
       do i=nnt,nct-1
-        radius(i)=pstok
+        radius(i)=pstok(1)
       enddo
 !  Load side chain radii
       do i=nnt,nct
         iti=itype(i,1)
-        radius(i+nres)=restok(iti)
+        radius(i+nres)=restok(iti,1)
       enddo
 !      do i=1,2*nres
 !        write (iout,*) "i",i," radius",radius(i)