C
C $Id: escf.F,v 1.8 1999/04/06 20:43:27 arjan Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE ESCF(EELECT1)
C
C     DETERMINES THE ELECTRONIC ENERGY FOR THE CURRENT SCF ITERATION.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"


      data tnow /0.0d0/
      save tnow

      if (wrtscr) then
         tlast = tnow
         call etimer(tnow)
         tcycle = tnow-tlast
C-RDC          write(iscr,'(" time for cycle=",F11.3)')tcycle
         if(icall.ne.0)then
            tcalls(icall) = tcalls(icall) + tcycle
            ncalls(icall) = ncalls(icall) + 1
         endif
      endif

C
      EELECT1 = 0.0D0
      IIMAX = IIMAT(NATOMS+1)-1
      IJMAX = IJMAT(IP1(NATOMS+1))-1
      DO 400 I=1,IIMAX
        EELECT1 = EELECT1 + (HDIAG(I) + FDIAG(I))*PDIAG(I)
 400  CONTINUE
      EII = 0.0D0
      II = 0
      DO 500 I=1,NATOMS
        IAI = IATNUM(I)
        NORBSI = NATORB(IAI)
C
C       SKIP ATOM IF IT'S A DUMMY OR A SPARKLE.
C
        IF(NORBSI.EQ.0) GO TO 500
        DO 450 IORB=1,NORBSI
          II = II + IORB
          EII = EII + (HDIAG(II) + FDIAG(II))*PDIAG(II)
 450    CONTINUE
 500  CONTINUE
      EELECT1 = EELECT1 - 0.5D0*EII
      DO 600 IJ=1,IJMAX
        EELECT1 = EELECT1 + (HDIAT(IJ) + FDIAT(IJ))*PDIAT(IJ)
 600  CONTINUE

      if (pmeqm) then
         ijphi=0
         epmeqm=0.0d0

         do i=1,natoms
         IAI = IATNUM(I)
                 
            DO 234 J=1,NATOMS
               ijphi=ijphi+1
               if (chewald) then
                  epmeqm = epmeqm + PHIPME(IJPHI)*chgpme(j)*chgpme(i)
               elseif (mullewald) then
                  epmeqm = epmeqm + PHIPME(IJPHI)*atchg(j)*atchg(i)
               elseif (cmewald) then
                  epmeqm = epmeqm + PHIPME(IJPHI)*atchg2(j)*atchg2(i)
               else 
                  write(*,*)'definition of the set of charges (escf1)'
                  stop
               endif
c            epmeqm = epmeqm + PHIPME(IJPHI)*atchg(j)*zchg(iai)
 234        CONTINUE
c                write(*,*)atchg(i)
!         write(*,*) 'epmeqm', epmeqm, 'atchg(1)', atchg(1), 'atchg2(1)',
!     $              atchg2(1)

         enddo
         EELECT1 = EELECT1 + 0.5d0*epmeqm
        
         virlrqm = 0.5d0*epmeqm
         ijphi=0
         do i=1,natoms
                 
            DO 235 J=1,NATOMS
               ijphi=ijphi+1
               if (chewald) then
                  virlrqm = virlrqm-0.5d0*PHISR(IJPHI)*chgpme(j)*
     &                chgpme(i)
               elseif (mullewald) then
                  virlrqm = virlrqm-0.5d0*PHISR(IJPHI)*atchg(j)*atchg(i)
               elseif (cmewald) then
                  virlrqm = virlrqm - 0.5d0*PHISR(IJPHI)*atchg2(j)*
     &                atchg2(i)
               else 
                  write(*,*)'definition of the set of charges (escf2)'
                  stop
               endif
c            epmeqm = epmeqm + PHIPME(IJPHI)*atchg(j)*zchg(iai)
 235        CONTINUE
c                write(*,*)atchg(i)

         enddo
      endif

      RETURN
      END