      REAL FUNCTION ELOCAL(CONFIG)
C calculate the local energy for CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variable:
      parameter(ncoord=6)
C Input variables:
      REAL CONFIG(NCOORD)        !configuration
C Local variables:
      REAL TPOP,VPOP             !kinetic and potential contributions
      REAL EECORR                !elec-elec correlation
      REAL CROSS,CROSS1,CROSS2   !cross terms
      REAL ONEE1,ONEE2           !one electron terms
      REAL X1,X2,Y1,Y2,Z1,Z2     !coordinates of 2 electrons
      REAL R1L,R1R,R2L,R2R,R12   !relative distances
      REAL CHI1,CHI2,F           !parts of the wave function
      REAL DOTR1L,DOTR2L,DOTR1R,DOTR2R !dot products with R12
      REAL R12DR1,SR12Z          !temp vars for dot products
      REAL CHI,FDCHI,SDCHI,LAPCHI !atomic orbitals
      REAL FEE,FDFEE,SDFEE,LAPFEE !elec-elec correlations
      REAL DIST                   !Euclidean distance
      REAL R,X,Y,Z                !dummy variables
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     define functions
      CHI(R)=EXP(-R/A)                        !atomic orbital
      FDCHI(R)=-CHI(R)/A                      !its first derivative,
      SDCHI(R)=CHI(R)/A/A                     !second derivative,
      LAPCHI(R)=SDCHI(R)+2*FDCHI(R)/R         !and Laplacian
C
      FEE(R)=EXP(R/(ALPHA*(1+BETA*R)))        !elec-elec correlation
      FDFEE(R)=FEE(R)/(ALPHA*(1.+BETA*R)**2)  !its first,second deriv,
      SDFEE(R)=FDFEE(R)**2/FEE(R)-2.*BETA*FEE(R)/ALPHA/(1+BETA*R)**3
      LAPFEE(R)=SDFEE(R)+2*FDFEE(R)/R         !and Laplacian
C
      DIST(X,Y,Z)=SQRT(X**2+Y**2+Z**2)        !Euclidean distance
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     get coordinates and radii
      CALL RADII(X1,X2,Y1,Y2,Z1,Z2,R1L,R2L,R1R,R2R,R12,CONFIG)
C                                                             
C     calculate dot products with R12
      R12DR1=X1*(X1-X2)+Y1*(Y1-Y2)+Z1*(Z1-Z2) !convenient starting place
      SR12Z=S*(Z1-Z2)/2                !useful constant
      DOTR1L=R12DR1+SR12Z              !dot products with R12
      DOTR1R=R12DR1-SR12Z
      DOTR2L=DOTR1L-R12**2
      DOTR2R=DOTR1R-R12**2
      DOTR1L=DOTR1L/R12/R1L            !dot products of unit vectors
      DOTR2L=DOTR2L/R12/R2L
      DOTR1R=DOTR1R/R12/R1R
      DOTR2R=DOTR2R/R12/R2R
C
      CHI1=CHI(R1R)+CHI(R1L)           !pieces of the total wavefunction
      CHI2=CHI(R2R)+CHI(R2L)
      F=FEE(R12)
C
      EECORR=2*LAPFEE(R12)/F           !correlation contribution
      ONEE1=(LAPCHI(R1L)+LAPCHI(R1R))/CHI1               !electron one
      ONEE2=(LAPCHI(R2L)+LAPCHI(R2R))/CHI2               !electron two
      CROSS1=(FDCHI(R1L)*DOTR1L+FDCHI(R1R)*DOTR1R)/CHI1  !cross terms
      CROSS2=(FDCHI(R2L)*DOTR2L+FDCHI(R2R)*DOTR2R)/CHI2
      CROSS=2*FDFEE(R12)*(CROSS1-CROSS2)/F
C
      TPOP=-HBM*(EECORR+ONEE1+ONEE2+CROSS)/2.0              !kinetic
      VPOP=-E2*(1./R1L + 1./R1R + 1./R2L + 1./R2R - 1./R12) !potential
      ELOCAL=TPOP+VPOP+E2/s                            !total
C
      RETURN
      END

