+#ifdef FIVEDIAG
+ subroutine inertia_tensor
+ use comm_gucio
+ use energy_data
+ real(kind=8) Im(3,3),Imcp(3,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,M_PEP
+ integer iti,inres,i,j,k,mnum,mnum1
+ 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
+ M_PEP=0.0d0
+
+!c caulating the center of the mass of the protein
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=0.0d0
+ M_PEP=M_PEP+mp(mnum)
+
+ do j=1,3
+ cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
+ enddo
+ enddo
+! do j=1,3
+! cm(j)=mp*cm(j)
+! enddo
+ M_SC=0.0d0
+ do i=nnt,nct
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ iti=iabs(itype(i,mnum))
+ if (iti.eq.ntyp1_molec(mnum)) cycle
+ M_SC=M_SC+msc(iabs(iti),mnum)
+ inres=i+nres
+ do j=1,3
+ cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/(M_SC+M_PEP)
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=0.0d0
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ do j=1,3
+ pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
+ enddo
+ Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
+ Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
+ Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
+ Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
+ Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
+ Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
+ enddo
+
+ do i=nnt,nct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (iti.eq.ntyp1_molec(mnum)) cycle
+ inres=i+nres
+ do j=1,3
+ pr(j)=c(j,inres)-cm(j)
+ enddo
+ Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
+ Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
+ Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
+ Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
+ Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
+ Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ iti=iabs(itype(i,mnum))
+ if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then
+ inres=i+nres
+ Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*&
+ dc_norm(1,inres))*vbld(inres)*vbld(inres)
+ Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*&
+ dc_norm(2,inres))*vbld(inres)*vbld(inres)
+ Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*&
+ dc_norm(2,inres))*vbld(inres)*vbld(inres)
+ Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ endif
+ enddo
+
+ call angmom(cm,L)
+ Im(2,1)=Im(1,2)
+ Im(3,1)=Im(1,3)
+ Im(3,2)=Im(2,3)
+
+!c Copng 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)
+ 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
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+! write (iout,*) itype(i+1,mnum1),itype(i,mnum)
+ if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) &
+ .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.&
+ itype(i,mnum).ne.ntyp1_molec(mnum) &
+ .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ 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
+ mnum=molnum(i)
+ if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) 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)
+ return
+ end subroutine inertia_tensor
+!------------------------------------------------------------
+ subroutine angmom(cm,L)
+ implicit none
+ double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),&
+ pp(3),mscab
+ integer iti,inres,i,j,mnum,mnum1
+!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
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=0.0d0
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ 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(mnum)*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(mnum)*vp(j)
+ enddo
+ enddo
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ inres=i+nres
+ do j=1,3
+ pr(j)=c(j,inres)-cm(j)
+ enddo
+ if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) 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)
+! if (mnum.gt.4) then
+! mscab=0.0d0
+! else
+ mscab=msc(iti,mnum)
+! endif
+ do j=1,3
+ L(j)=L(j)+mscab*vp(j)
+ enddo
+!c write (iout,*) "L",(l(j),j=1,3)
+ if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) 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,mnum)*vp(j)
+ enddo
+ endif
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+ return
+ end subroutine angmom
+!---------------------------------------------------------------
+ subroutine vcm_vel(vcm)
+ double precision vcm(3),vv(3),summas,amas
+ integer i,j,mnum,mnum1
+ do j=1,3
+ vcm(j)=0.0d0
+ vv(j)=d_t(j,0)
+ enddo
+ summas=0.0d0
+ do i=nnt,nct
+ mnum=molnum(i)
+ if ((mnum.gt.5).or.(mnum.eq.3))&
+ 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.4) then
+ amas=msc(iabs(itype(i,mnum)),mnum)
+! else
+! amas=0.0d0
+! endif
+! amas=msc(iabs(itype(i)))
+ summas=summas+amas
+! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.lt.4)) 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
+ write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
+ do j=1,3
+ vcm(j)=vcm(j)/summas
+ enddo
+ return
+ end subroutine vcm_vel
+#else