+#ifdef FIVEDIAG
+ if (lprn) then
+ write (iout,*) "Potential forces backbone"
+ do i=nnt,nct
+ write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
+ enddo
+ write (iout,*) "Potential forces sidechain"
+ do i=nnt,nct
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
+ write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
+ enddo
+ endif
+ do j=1,3
+ ind=1
+ do i=nnt,nct
+ if (itype(i).eq.10 .or. itype(i).eq.ntyp1)then
+ rs(ind)=-gcart(j,i)-gxcart(j,i)
+ ind=ind+1
+ else
+ rs(ind)=-gcart(j,i)
+ rs(ind+1)=-gxcart(j,i)
+ ind=ind+2
+ end if
+ enddo
+#ifdef CHECK5DSOL
+ rsold=rs
+#endif
+ if (lprn) then
+ write(iout,*) "RHS of the 5-diag equations system"
+ do i=1,dimen
+ write(iout,*) i,rs(i)
+ enddo
+ endif
+
+ call FDISYS (dimen,DM,DU1,DU2,rs,xsolv)
+ if (lprn) then
+ write (iout,*) "Solution of the 5-diagonal equations system"
+ do i=1,dimen
+ write (iout,'(i5,f10.5)') i,xsolv(i)
+ enddo
+ endif
+#ifdef CHECK5DSOL
+! Check the solution
+ rscheck(1)=DMorig(1)*xsolv(1)+DU1orig(1)*xsolv(2)+&
+ DU2orig(1)*xsolv(3)
+ rscheck(2)=DU1orig(1)*xsolv(1)+DMorig(2)*xsolv(2)+&
+ DU1orig(2)*xsolv(3)+DU2orig(2)*xsolv(4)
+
+ do i=3,dimen-2
+ rscheck(i)=DU2orig(i-2)*xsolv(i-2)+DU1orig(i-1)*&
+ xsolv(i-1)+DMorig(i)*xsolv(i)+DU1orig(i)*&
+ xsolv(i+1)+DU2orig(i)*xsolv(i+2)
+ enddo
+ rscheck(dimen-1)=DU2orig(dimen-3)*xsolv(dimen-3)+&
+ DU1orig(dimen-2)*xsolv(dimen-2)+DMorig(dimen-1)*&
+ xsolv(dimen-1)+DU1orig(dimen-1)*&
+ xsolv(dimen)
+ rscheck(dimen)=DU2orig(dimen-2)*xsolv(dimen-2)+DU1orig(dimen-1)*&
+ xsolv(dimen-1)+DMorig(dimen)*xsolv(dimen)
+
+ do i=1,dimen
+ write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),&
+ "ratio",rscheck(i)/rsold(i)
+ enddo
+! end check
+#endif
+ ind=1
+ do i=nnt,nct
+ if (itype(i).eq.10 .or. itype(i).eq.ntyp1)then
+ d_a(j,i)=xsolv(ind)
+ ind=ind+1
+ else
+ d_a(j,i)=xsolv(ind)
+ d_a(j,i+nres)=xsolv(ind+1)
+ ind=ind+2
+ end if
+ enddo
+ enddo
+ if (lprn) then
+ do i=nnt,nct
+ write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
+ enddo
+ do i=nnt,nct
+ write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
+ enddo
+ endif
+ do j=1,3
+ d_a(j,0)=d_a(j,nnt)
+ enddo
+ do i=nnt,nct
+ if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
+ do j=1,3
+ d_a(j,i)=d_a(j,i+1)-d_a(j,i)
+ enddo
+ else
+ do j=1,3
+ d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
+ d_a(j,i)=d_a(j,i+1)-d_a(j,i)
+ enddo
+
+
+ end if
+ enddo
+
+ if (lprn) then
+ write(iout,*) 'acceleration 3D FIVEDIAG'
+ write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
+ do i=nnt,nct-1
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
+ enddo
+ do i=nnt,nct
+ write (iout,'(i3,3f10.5,3x,3f10.5)') &
+ i+nres,(d_a(j,i+nres),j=1,3)
+ enddo
+ endif
+#else
+! Old procedure