C***********************************************************************
C FILE: FDIRTR.FOR
C
C HISTORY: 23/2/90: COPIED FROM FILE FDIR16.FOR
C                   TRICLINIC GEOMETRY
C          25/2/90: OFF-DIAGONAL ELEMENTS IN STRESS TENSOR
C           8/3/90: SIGN OF STRESS-CONTRIB. CHANGED
C          18/5/90: DIRDAT() FOR USE WITH CONJUGATE GRADIENT
C***********************************************************************
      SUBROUTINE FDIRTR(NGX,NGY,NGZ,NPLWV,RECC,VOLC,CVD,CHDEN,
     &     SIGHA,ENHAC,LPCTX,LPCTY,LPCTZ,DIRDAT)
C=======================================================================
C
C{{{{{{{{{{{{{{{{{{{{{{{{ SUBROUTINE FDIRTR }}}}}}}}}}}}}}}}}}}}}}}}}}}}
C
C THIS SUBROUTINE CALCULATES THE HARTREE POTENTIAL FROM THE ELECTRONIC
C CHARGE DENSITY CALCULATED IN SUBROUTINE CHSP. THE CORRECTION TO THE
C TOTAL ENERGY DUE TO OVERCOUNTING THE HARTREE ENERGY ON SUMMING THE
C ELECTRONIC EIGENVALUES IS ALSO COMPUTED (HARTREE CONTRIBUTION TO THE
C TOTAL ENERGY = 0.5*SUM (VH(G)*RHO(-G)). THE SUM OF EIGENVALUES GIVES
C SUM (VH(G)*RHO(-G)) WHERE VH(G) IS THE HARTREE POTENTIAL AT WAVEVECTOR
C G AND RHO(G) IS THE CHARGE DENSITY AT WAVEVECTOR G)
C
C                                               G  G
C                   1        --   ! n(G) ! [     Ai Aj            ]
C     Sigma      = --- 4 PI  >    -------- [ 2 --------  - &      ]
C          Ai,Aj    2        --     G**2   [     G**2       Ai,Aj ]
C                           G,G<>0
C
C=======================================================================
      IMPLICIT COMPLEX (C)
      DIMENSION RECC(3,*)
      DIMENSION CVD(*)
      DIMENSION DIRDAT(*)
      DIMENSION CHDEN(*)
      DIMENSION SIGHA(*)
      DIMENSION LPCTX(*),LPCTY(*),LPCTZ(*)
C=======================================================================
C EDEPS = ELECTRON CHARGE DIVIDED BY THE PERMITTIVITY OF FREE SPACE
C=======================================================================
      DATA TPI,EDEPS / 6.28318530718,180.95283 /
      DO 8620 M=1,6
      SIGHA(M)=0.0
 8620 CONTINUE
      ENHAC=0.0
      SCALE=EDEPS/VOLC
C=======================================================================
C CALCULATE THE HARTREE POTENTIAL ON THE GRID OF RECIPROCAL LATTICE
C VECTORS, THE CORRECTION TO THE TOTAL ENERGY AND THE FORCE ON THE UNIT
C CELL
C=======================================================================
      NI=1
      DO 8630 N=1,NGZ
        GZX=RECC(3,1)*LPCTZ(N)
        GZY=RECC(3,2)*LPCTZ(N)
        GZZ=RECC(3,3)*LPCTZ(N)
        DO 8631 NN=1,NGY
          GYX=RECC(2,1)*LPCTY(NN)
          GYY=RECC(2,2)*LPCTY(NN)
          GYZ=RECC(2,3)*LPCTY(NN)
          DO 8632 NNN=1,NGX
            GXX=RECC(1,1)*LPCTX(NNN)
            GXY=RECC(1,2)*LPCTX(NNN)
            GXZ=RECC(1,3)*LPCTX(NNN)
            GISQ=(GXX+GYX+GZX)**2+(GXY+GYY+GZY)**2+(GXZ+GYZ+GZZ)**2
C=======================================================================
C SINCE THE G=0 COULOMB CONTRIBUTIONS TO THE HARTREE, EWALD AND
C ELECTRON-ION ENERGIES ARE INDIVIDUALLY DIVERGENT BUT TOGETHER SUM TO
C ZERO, SET THE HARTREE POTENTIAL AT G=0 TO ZERO.
C=======================================================================
            IF (GISQ.LT.0.001) THEN
              CVD(NI)=(0.0,0.0)
              DIRDAT(NI)=0.0
            ELSE
              GISQ=1.0/GISQ
              CVD(NI)=CHDEN(NI)*GISQ
              DIRDAT(NI)=SCALE*GISQ
              DUM=CHDEN(NI)*CONJG(CHDEN(NI))*GISQ
              ENHAC=ENHAC+DUM
              G1=GXX+GYX+GZX
              G2=GXY+GYY+GZY
              G3=GXZ+GYZ+GZZ
              SIGHA(1)=SIGHA(1)+DUM*(2*(G1**2)*GISQ-1)
              SIGHA(2)=SIGHA(2)+DUM*(2*(G2**2)*GISQ-1)
              SIGHA(3)=SIGHA(3)+DUM*(2*(G3**2)*GISQ-1)
              SIGHA(4)=SIGHA(4)+DUM*(2*(G2*G3)*GISQ)
              SIGHA(5)=SIGHA(5)+DUM*(2*(G3*G1)*GISQ)
              SIGHA(6)=SIGHA(6)+DUM*(2*(G1*G2)*GISQ)
            ENDIF
            NI=NI+1
 8632     CONTINUE
 8631   CONTINUE
 8630 CONTINUE
C=======================================================================
C SCALE THE HARTREE POTENTIAL BY EDEPS DIVIDED BY THE VOLUME OF THE UNIT
C CELL
C=======================================================================
      SCALE=EDEPS/VOLC
      DO 8640 N=1,NPLWV
        CVD(N)=CVD(N)*SCALE
 8640 CONTINUE
C=======================================================================
C SCALE FORCE ON THE UNIT CELL (INCLUDING THE FACTOR OF 1/2 FOR THE
C SELF-ENERGY FORM OF THE HARTREE ENERGY)
C=======================================================================
      SIGHA(1)=SIGHA(1)*SCALE/(VOLC*2)
      SIGHA(2)=SIGHA(2)*SCALE/(VOLC*2)
      SIGHA(3)=SIGHA(3)*SCALE/(VOLC*2)
      SIGHA(4)=SIGHA(4)*SCALE/(VOLC*2)
      SIGHA(5)=SIGHA(5)*SCALE/(VOLC*2)
      SIGHA(6)=SIGHA(6)*SCALE/(VOLC*2)
C=======================================================================
C SCALE THE CORRECTION TO THE TOTAL ENERGY WHICH IS REQUIRED TO TAKE
C ACCOUNT OF OVERCOUNTING THE HARTREE ENERGY IN THE ELECTRON ENERGIES
C=======================================================================
      ENHAC=-ENHAC*SCALE/2
      RETURN
      END
