+ write (iout,*) "nnt",nnt," nct",nct," nside",nside
+#ifdef FIVEDIAG
+#ifdef CRYST_BOND
+ ip4=ip/4
+ DM(1)=mp/4+ip4
+ if (iabs(itype(nnt)).eq.10) then
+ DM(1)=DM(1)+msc(10)
+ ind=2
+ nind=1
+ else
+ DM(1)=DM(1)+isc(iabs(itype(nnt)))
+ DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt)))
+ ind=3
+ nind=2
+ endif
+ do i=nnt+1,nct-1
+! if (iabs(itype(i,1)).eq.ntyp1) cycle
+ DM(ind)=2*ip4+mp/2
+ if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
+ if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
+ ind=ind+1
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
+ DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
+ ind=ind+2
+ endif
+ enddo
+ if (nct.gt.nnt) then
+ DM(ind)=ip4+mp/4
+ if (iabs(itype(nct)).eq.10) then
+ DM(ind)=DM(ind)+msc(10)
+ nind=ind
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(nct)))
+ DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct)))
+ nind=ind+1
+ endif
+ endif
+
+
+ ind=1
+ do i=nnt,nct
+ mnum=molnum(i)
+ if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
+ .and.mnum.eq.5) then
+ DU1(ind)=-isc(iabs(itype(i,1)))
+ DU1(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU1(ind)=mp/4-ip4
+ ind=ind+1
+ endif
+ enddo
+
+ ind=1
+ do i=nnt,nct-1
+ mnum=molnum(i)
+! if (iabs(itype(i,1)).eq.ntyp1) cycle
+ write (iout,*) "i",i," itype",itype(i,1),ntyp1
+ if (iabs(itype(i,1)).ne.10 .and. &
+ iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
+ DU2(ind)=mp/4-ip4
+ DU2(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU2(ind)=0.0d0
+ DU2(ind+1)=0.0d0
+ ind=ind+1
+ endif
+ enddo
+#else
+ ip4=ip(1)/4
+ DM(1)=mp(1)/4+ip4
+ if (iabs(itype(nnt,1)).eq.10) then
+ DM(1)=DM(1)+msc(10,1)
+ ind=2
+ nind=1
+ else
+ DM(1)=DM(1)+isc(iabs(itype(nnt,1)),1)
+ DM(2)=msc(iabs(itype(nnt,1)),1)+isc(iabs(itype(nnt,1)),1)
+ ind=3
+ nind=2
+ endif
+ do i=nnt+1,nct-1
+! if (iabs(itype(i,1)).eq.ntyp1) cycle
+ DM(ind)=2*ip4+mp(1)/2
+ if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
+ if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10,1)
+ ind=ind+1
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(i,1)),1)
+ DM(ind+1)=msc(iabs(itype(i,1)),1)+isc(iabs(itype(i,1)),1)
+ ind=ind+2
+ endif
+ enddo
+ if (nct.gt.nnt) then
+ DM(ind)=ip4+mp(1)/4
+ if (iabs(itype(nct,1)).eq.10) then
+ DM(ind)=DM(ind)+msc(10,1)
+ nind=ind
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(nct,1)),1)
+ DM(ind+1)=msc(iabs(itype(nct,1)),1)+isc(iabs(itype(nct,1)),1)
+ nind=ind+1
+ endif
+ endif
+
+
+ ind=1
+ do i=nnt,nct
+ mnum=molnum(i)
+ if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1 &
+ .and.mnum.eq.5) then
+ DU1(ind)=-isc(iabs(itype(i,1)),1)
+ DU1(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU1(ind)=mp(1)/4-ip4
+ ind=ind+1
+ endif
+ enddo
+
+ ind=1
+ do i=nnt,nct-1
+ mnum=molnum(i)
+! if (iabs(itype(i,1)).eq.ntyp1) cycle
+ write (iout,*) "i",i," itype",itype(i,1),ntyp1
+ if (iabs(itype(i,1)).ne.10 .and. &
+ iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
+ DU2(ind)=mp(1)/4-ip4
+ DU2(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU2(ind)=0.0d0
+ DU2(ind+1)=0.0d0
+ ind=ind+1
+ endif
+ enddo
+#endif
+ DMorig=DM
+ DU1orig=DU1
+ DU2orig=DU2
+ write (iout,*) "nind",nind," dimen",dimen
+ if (nind.ne.dimen) then
+ write (iout,*) "ERROR, dimensions differ"
+#ifdef MPI
+ call MPI_Finalize(ierr)
+#endif
+ stop
+ endif
+ write (iout,*) "The upper part of the five-diagonal inertia matrix"
+ do i=1,dimen
+ if (i.lt.dimen-1) then
+ write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
+ else if (i.eq.dimen-1) then
+ write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
+ else
+ write (iout,'(i3,3f10.5)') i,DM(i)
+ endif
+ enddo
+
+ call FDISYP (dimen, DM, DU1, DU2, MARK)
+
+ if (mark.eq.-1) then
+ write (iout,*) "ERROR: the inertia matrix is not positive definite"
+#ifdef MPI
+ call MPI_Finalize(ierr)
+#endif
+ stop
+ else if (mark.eq.0) then
+ write (iout,*) "ERROR: the inertia matrix is singular"
+#ifdef MPI
+ call MPI_Finalize(ierr)
+#endif
+ else if (mark.eq.1) then
+ write (iout,*) "The transformed inertia matrix"
+ do i=1,dimen
+ if (i.lt.dimen-1) then
+ write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
+ else if (i.eq.dimen-1) then
+ write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
+ else
+ write (iout,'(i3,3f10.5)') i,DM(i)
+ endif
+ enddo
+ endif
+! Diagonalization of the pentadiagonal matrix
+ allocate(DDU1(2*nres))
+ allocate(DDU2(2*nres))
+ allocate(DDM(2*nres))
+ do i=1,dimen-1
+ DDU1(i+1)=DU1orig(i)
+ enddo
+ do i=1,dimen-2
+ DDU2(i+2)=DU2orig(i)
+ enddo
+ DDM=DMorig
+ call quindibisect2(DDM,DDU1,DDU2,dimen,1,dimen,eps1,relfeh,eps2,iz,geigen)
+ if (lprn) then
+ write (iout,*) &
+ "Eigenvalues of the five-diagonal inertia matrix"
+ do i=1,dimen
+ write (iout,'(i5,f10.5)') i,geigen(i)
+ enddo
+ endif
+ deallocate(DDU1)
+ deallocate(DDU2)
+ deallocate(DDM)
+#else
+! Old Gmatrix