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 .and. itype(i).ne.21) 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 .and. itype(i).ne.21) 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 .and. itype(i).ne.21) 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 .and. itype(i).ne.21) 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 .and. itype(i).ne.21) 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