+#ifdef FIVEDIAG
+c Here accelerations due to friction forces are computed right after forces.
+ d_t_work=0.0d0
+ do j=1,3
+ v_work(j,1)=d_t(j,0)
+ v_work(j,nnt)=d_t(j,0)
+ enddo
+ do i=nnt+1,nct
+ do j=1,3
+ v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
+ do j=1,3
+ v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
+ enddo
+ endif
+ enddo
+#ifdef DEBUG
+ write (iout,*) "v_work"
+ do i=1,2*nres
+ write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
+ enddo
+#endif
+ do j=1,3
+ ind=0
+ do ichain=1,nchain
+ n=dimen_chain(ichain)
+ iposc=iposd_chain(ichain)
+c write (iout,*) "friction_force j",j," ichain",ichain,
+c & " n",n," iposc",iposc,iposc+n-1
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+c diagnostics
+c innt=chain_border(1,1)
+c inct=chain_border(2,1)
+ do i=innt,inct
+ vvec(ind+1)=v_work(j,i)
+ ind=ind+1
+ if (iabs(itype(i)).ne.10) then
+ vvec(ind+1)=v_work(j,i+nres)
+ ind=ind+1
+ endif
+ enddo
+#ifdef DEBUG
+ write (iout,*) "vvec ind",ind," n",n
+ write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
+#endif
+c write (iout,*) "chain",i," ind",ind," n",n
+ call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
+ & DU2fric(iposc),vvec(iposc),rs)
+#ifdef DEBUG
+ write (iout,*) "rs"
+ write (iout,'(f10.5)') (rs(i),i=1,n)
+#endif
+ do i=iposc,iposc+n-1
+c write (iout,*) "ichain",ichain," i",i," j",j,
+c & "index",3*(i-1)+j,"rs",rs(i-iposc+1)
+ fric_work(3*(i-1)+j)=-rs(i-iposc+1)
+ enddo
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "Vector fric_work dimen3",dimen3
+ write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
+#endif
+#else