SUBROUTINE FCCE(nnorbs)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
      DATA BOLTZ /8.617335408D-5/
      SAVE BOLTZ
C

       
      BETAH = 1.0D0/(BOLTZ*TEMPK)
      nfull = 0
      nfroz = 0
        NPASS = 1
        DOEVAL = IPASS.EQ.1
          do K=1,1
       
C
C         DIAGONAL BLOCKS FIRST.
C
          NORBSK = IORBPT(K+1)-IORBPT(K)
          IIORB = 0
          IJLOC0 = 0
          DO 100 I=IATOM1(K),IATOM1(K+1)-1
            IATM = IATOMS(I)
            NORBSI = NATORB(IATNUM(IATM))
            IJGLB = IIMAT(IATM)
            DO 20 IF=1,NORBSI
              IJLOC = IJLOC0 + IF + IIORB
              DO 10 JF=1,IF
                FF(IJLOC) = FDIAG(IJGLB)
                IJGLB = IJGLB + 1
                IJLOC = IJLOC + NORBSK
 10           CONTINUE
 20         CONTINUE
            IIORB = IIORB + NORBSI
            IJLOC0 = IIORB*NORBSK
 100      CONTINUE
C
C         NOW OFF-DIAGONAL BLOCKS.
C
          NPAIRS = IP1(NATOMS+1)-1
          NATMS = IATOM1(K+1) - IATOM1(K)
          IF(NATMS.GT.1)THEN
            IATM0 = IATOM1(K) - 1
            IATM1 = IATOMS(IATOM1(K))
            II = NATORB(IATNUM(IATM1))
            DO 200 I=2,NATMS
              IATM = IATOMS(IATM0+I)
              NORBSI = NATORB(IATNUM(IATM))
              JJ = 0
              IJLOC0 = 0
              DO 180 J=1,I-1
                JATM = IATOMS(IATM0+J)
C
C               GET POSITION IN PAIRLIST OF (IATM,JATM) PAIR.
C
                CALL IJFIND(NPAIRS,IATM,JATM,IJADDR)
                NORBSJ = NATORB(IATNUM(JATM))
C
C               COPY GLOBAL DIATOMIC BLOCK FDIAT TO FF, BUT ONLY IF IATM
C               AND JATM ARE CLOSE ENOUGH TO BOND.  OTHERWISE, SET THE
C               BLOCK OF FF TO ZERO.
C
                IF(IJADDR.NE.0)THEN
                  IJGLB = IJMAT(IJADDR)
                  DO 130 IF=1,NORBSI
                    IJLOC = IJLOC0 + II + IF
                    DO 120 JF=1,NORBSJ
                      FF(IJLOC) = FDIAT(IJGLB)
                      IJLOC = IJLOC + NORBSK
                      IJGLB = IJGLB + 1
 120                CONTINUE
 130              CONTINUE
                ELSE
                  DO 150 IF=1,NORBSI
                    IJLOC = IJLOC0 + II + IF
                    DO 140 JF=1,NORBSJ
                      FF(IJLOC) = 0.0D0
                      IJLOC = IJLOC + NORBSK
 140                CONTINUE
 150              CONTINUE
                ENDIF
                JJ = JJ + NORBSJ
                IJLOC0 = JJ*NORBSK
 180          CONTINUE
              II = II + NORBSI
 200        CONTINUE
          ENDIF
          enddo
c             write(*,*)'FF',(ff(i),i=1,nnorbs*nnorbs)
c             write(*,*)'EVEC',(EVEC(i),i=1,nnorbs*nnorbs)
c             write(*,*)'eval',(eval(i),i=1,nnorbs)
          call MULTIMAT(nnorbs,FF,EVEC,eval)

           END

         SUBROUTINE MULTIMAT(ndim,AAA,BBB,CCC) 

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C#include "divcon.dim"
C#include "divcon.h"
       DIMENSION AAA(*),BBB(*),cCC(*)
       DIMENSION ediff1(ndim,ndim),ediff2(ndim,ndim),ediff(ndim,ndim)
       DIMENSION cccc(ndim,ndim),bbbb(ndim,ndim),aaaa(ndim,ndim)
       DIMENSION rnorm(ndim,ndim)

        do j=1,ndim
         do i=1,ndim
           if (i.eq.j) then
            cccc(j,i)=ccc(j)
           else
            cccc(j,i)=0.0d0
           endif
         enddo
        enddo
      ij=0
       do i=1,ndim
        do j=1,ndim
             ij=ij+1
          bbbb(i,j)=bbb(ij)
        enddo
       enddo  
       do i=1,ndim
        do j=1,i
           aaaa(i,j)=aaa(i+(j-1)*ndim) 
           aaaa(j,i)= aaaa(i,j)
        enddo
       enddo 
c      ij=0
c       do i=1,ndim
c        do j=1,i
c             ij=ij+1
c          bbbb(i,j)=bbb(ij)
c          bbbb(j,i)=bbb(ij)
c        enddo
c       enddo


       do i=1,NDIM
         do j=1,ndim
           ediff1(i,j)=0.0d0
             do k=1,ndim
               ediff1(i,j)=ediff1(i,j)+AAAA(j,k)*BBBB(i,k)
             enddo
        enddo
       enddo

       do i=1,NDIM
        do j=1,ndim
           ediff2(i,j)=0.0d0
             do k=1,ndim
               ediff2(i,j)=ediff2(i,j)+CCCC(j,k)*BBBB(k,i)
             enddo
        enddo
       enddo
         ediffmax=1.0d-16
       do i=1,NDIM
        do j=1,ndim
          ediff(i,j)=ediff1(i,j)-ediff2(j,i)
         if (dabs(ediff(i,j)).gt.ediffmax) ediffmax=dabs(ediff(i,j))
        enddo
       enddo
c           write(*,*)
            write(*,*)'max(FC-Ce)=',ediffmax
c          do i=1,NDIM
c           write(*,*)(ediff(j,i),j=1,ndim)
c          enddo

         do i=1,ndim
          do j=1,ndim
           rnorm(i,j)=0.0d0
           do k=1,ndim
            rnorm(i,j)=rnorm(i,j)+bbbb(k,i)*bbbb(k,j)
           enddo
          enddo
         enddo
 

         

           write(1,*)
c            write(1,*)'FC-Ce'
          do j=1,NDIM
           write(1,'(12f14.11)')(rnorm(j,i),i=1,ndim)
          enddo
c           write(2,*)
c            write(2,*)'FC-Ce'
cc          do i=1,NDIM
c           write(2,*)(ediff2(j,1),j=1,ndim)
cc          enddo
 99     format(36f10.7 ) 


         END