      program h2dist
      INCLUDE 'PARAM.P8'
      REAL CONFIG(NCOORD),csave(ncoord)
      real x1(1000),y1(1000),x2(1000),y2(1000)
      delta=1.0
      dseed=76453456.d0
      HBM=7.6359               !hbar*2/(mass)
      E2=14.409                !charge of the electron
      ABOHR=HBM/E2
      alpha=2.*abohr
      beta=0.25
      write(*,*) 'Zadaj separaciu s'
      read(*,*) sep
      s=sep
      s2=0.5*sep
      a=apar(ABOHR,s)
      accpt=0.0
      call intcfg(config,w1)
c      w1=phi(config)*phi(config)
      en=0.0
      n=0
30    continue      
      DO 10 ICOORD=1,NCOORD
       CSAVE(ICOORD)=CONFIG(ICOORD)        !save previous values
       CONFIG(ICOORD)=CONFIG(ICOORD)+DELTA*(RANNOS(dseed)-.5)!trial step
10    CONTINUE
      w2=phi(config)*phi(config)
      IF (w2/w1 .GT. RANNOS(DSEED)) THEN  !sometimes accept the step
         W1=W2                            !save new weight
         ACCPT=ACCPT+1.0                     !update accpt ratio
         n=n+1
         x1(n)=config(1)
         y1(n)=config(3)
         x2(n)=config(4)
         y2(n)=config(6)
         eloc=elocal(config)
         en=en+eloc
         if (accpt.EQ.990) goto 40
      ELSE
        DO 20 ICOORD=1,NCOORD
          CONFIG(ICOORD)=CSAVE(ICOORD)     !or else restore old config
20      CONTINUE
      END IF
c      write(*,*) 'separacia = ',sep,'  energia = ',en/accpt,'eV'
      goto 30
40    call metafl('cons')
      call disini()
      call height(50)
      call hname(50)
c      call labdig(2,'X')
c      call labdig(2,'Y')
c      call name('Separation (A)','X')
c      call name('Energy (eV)','Y')
      call graf(-3.0,3.0,-3.0,1.0,-3.0,3.0,-3.0,1.0)
      CALL HSYMBL(8)
      call color('red')
      do i=1,999
         call rlsymb(21,y1(i),x1(i))
      end do
      call color('yellow')
      do i=1,999
         call rlsymb(21,y2(i),x2(i))
      end do
      CALL HSYMBL(50)
      call color('blue')
      call rlsymb(21,-s2,0.0)
      call rlsymb(21,s2,0.0)
      call disfin
      end
c  
      SUBROUTINE INTCFG(CONFIG,W)
C generate a configuration (CONFIG) and calculate its weight (W)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variable:
      INCLUDE 'PARAM.P8'
C Output variables:
      REAL CONFIG(NCOORD)        !configuration
      REAL W                     !weight for single config
C Local variables:
      INTEGER ICOORD             !coordinate index
C Function:
      REAL PHI                   !total wavefunction
      REAL RANNOS                !uniform random number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 10 ICOORD=1,NCOORD      !pick configuration at random
         CONFIG(ICOORD)=A*(RANNOS(DSEED)-.5)
10    CONTINUE
      CONFIG(3)=CONFIG(3)+S2    !center elec 1. at right
      CONFIG(6)=CONFIG(6)-S2    !center elec 2. at left 
      W=PHI(CONFIG)**2          !weight=phi**2
      RETURN
      END
c
      SUBROUTINE METROP(CONFIG,W,ACCPT)
C take a Metropolis step
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Input/Output variables:
      REAL CONFIG(NCOORD)        !configuration 
      REAL W                     !weight for single config 
      REAL ACCPT                 !acceptance ratio 
C Local variables:
      INTEGER ICOORD             !coordinate index
      REAL CSAVE(NCOORD)         !temp storage for last config
      REAL WTRY                  !weight for trial config
C Function:
      REAL PHI                   !total wavefunction
      REAL RANNOS                !uniform random number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 10 ICOORD=1,NCOORD
       CSAVE(ICOORD)=CONFIG(ICOORD)        !save previous values
       CONFIG(ICOORD)=CONFIG(ICOORD)+DELTA*(RANNOS(DSEED)-.5)!trial step
10    CONTINUE
      WTRY=PHI(CONFIG)**2                  !trial weight
C
      IF (WTRY/W .GT. RANNOS(DSEED)) THEN  !sometimes accept the step
         W=WTRY                            !save new weight
         ACCPT=ACCPT+1                     !update accpt ratio
      ELSE
        DO 20 ICOORD=1,NCOORD
          CONFIG(ICOORD)=CSAVE(ICOORD)     !or else restore old config
20      CONTINUE
      END IF
      RETURN
      END     
c      
      REAL FUNCTION APAR(ABOHR,s)
      AOLD=0.
      apar=10.0
10    IF (ABS(APAR-AOLD) .GT. 1.E-6) THEN
         AOLD=APAR
         APAR=ABOHR/(1+EXP(-S/AOLD))
      GOTO 10
      END IF
      RETURN
      END
      REAL FUNCTION ELOCAL(CONFIG)
C calculate the local energy for CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variable:
      INCLUDE 'PARAM.P8'
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 
c
      SUBROUTINE RADII(X1,X2,Y1,Y2,Z1,Z2,R1L,R2L,R1R,R2R,R12,CONFIG)
C calculates cartesian coordinates and radii given CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variable:
      INCLUDE 'PARAM.P8'
C Input variables:
      REAL CONFIG(NCOORD)        !configuration
C Output variables:
      REAL X1,X2,Y1,Y2,Z1,Z2     !coordinates of 2 electrons
      REAL R1L,R1R,R2L,R2R,R12   !relative distances
      REAL DIST                  !Euclidean distance
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DIST(X1,Y1,Z1)=SQRT(X1**2+Y1**2+Z1**2)        !Euclidean distance
C
      X1=CONFIG(1)            !give the CONFIG elements their real names
      X2=CONFIG(4)
      Y1=CONFIG(2)
      Y2=CONFIG(5)
      Z1=CONFIG(3)
      Z2=CONFIG(6)
C
      R1L=DIST(X1,Y1,Z1+S2)   !calculate separations
      R1R=DIST(X1,Y1,Z1-S2)
      R2L=DIST(X2,Y2,Z2+S2)
      R2R=DIST(X2,Y2,Z2-S2)
      R12=DIST(X1-X2,Y1-Y2,Z1-Z2)
C
      RETURN
      END    
c
      REAL FUNCTION PHI(CONFIG)
C calculates the total variational wavefunction for CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Input variables:
      REAL CONFIG(NCOORD)        !configuration
C Local variables:                                               
      REAL X1,X2,Y1,Y2,Z1,Z2     !coordinates of 2 electrons
      REAL R1L,R1R,R2L,R2R,R12   !relative distances
      REAL CHI1R,CHI1L,CHI2R,CHI2L,F  !parts of the wave function
      REAL CHI,FEE               !terms in the wavefunction
      REAL R                     !radius
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CHI(R)=EXP(-R/A)                   !atomic orbital
      FEE(R)=EXP(R/(ALPHA*(1+BETA*R)))   !electron-electron correlation
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     calculate the radii
      CALL RADII(X1,X2,Y1,Y2,Z1,Z2,R1L,R2L,R1R,R2R,R12,CONFIG)
C      
      CHI1R=CHI(R1R)          !pieces of the total wavefunction
      CHI1L=CHI(R1L)
      CHI2R=CHI(R2R)
      CHI2L=CHI(R2L)
      F=FEE(R12)
      PHI=(CHI1L +CHI1R)*(CHI2L+CHI2R)*F   !the whole thing
C
      RETURN
      END  
c
      REAL FUNCTION RANNOS(DSEED)
C returns a uniformly distributed random number between 0 and 1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DOUBLE PRECISION         DSEED
      DOUBLE PRECISION         D2P31M,D2P31
      DATA             D2P31M/2147483647.D0/
      DATA             D2P31 /2147483711.D0/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DSEED = MOD(16807.D0*DSEED,D2P31M)
      RANNOS = DSEED / D2P31
      RETURN
      END