      program qh2mol
      INCLUDE 'PARAM.P8'
      REAL ENERGY(3,3)           !energy
      REAL WEIGHT(MAXENS)        !weight of ensemble members
      REAL ENSMBL(NCOORD,MAXENS) !ensemble of configurations
      REAL ESAVE(MAXCRR) 
      REAL CONFIG(NCOORD),csave(ncoord),x(100),y(100)
      delta=1.0
      dt=1.0e-5
      nensem=900
      dseed=76453456.d0
      ntherm=100
      HBM=7.6359               !hbar*2/(mass)
      E2=14.409                !charge of the electron
      ABOHR=HBM/E2
      alpha=2.*abohr
      beta=0.25
      n=0
      do sep=0.1,3.0,0.05
        n=n+1
        s=sep
        s2=0.5*sep
        a=apar(ABOHR,s)
        accpt=0.0
        CALL INTENS(ENSMBL,WEIGHT,CONFIG)
c       w1=phi(config)*phi(config)
        en=0.0
        DO 1 ITHERM=1,NTHERM
           CALL TSTEP(ENSMBL,WEIGHT,EPSILN)   ! thermalizing
 1      CONTINUE
        DO 2 IQUANT=1,3                       !zero total sums
           ENERGY(TOTAL,IQUANT)=0.
 2      CONTINUE
        en=0.0
        do i=1,1
            CALL TSTEP(ENSMBL,WEIGHT,EPSILN)
            en=en+epsiln
        end do
        en=en/real(1)
      write(*,*) 'separacia = ',sep,'  energia = ',en,'eV'
      x(n)=sep
      y(n)=en
      end do
      call metafl('cons')
      call disini()
      call height(50)
      call hname(50)
      call labdig(2,'X')
      call labdig(2,'Y')
      call name('Separation (A)','X')
      call name('Energy (eV)','Y')
      call graf(0.0,3.0,0.0,0.5,-32.5,-23.0,-32.5,1.0)
      call linwid(4)
      call color('red')
      call curve(x,y,n)
      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
c 
      SUBROUTINE TSTEP(ENSMBL,WEIGHT,EPSILN)
C take a time step using the Path Integral Monte Carlo method
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Input/Output variables
      REAL WEIGHT(MAXENS)       !weight of ensemble members (I/O)
      REAL ENSMBL(NCOORD,MAXENS)!ensemble of configurations (I/O)
      REAL EPSILN               !local energy of CONFIG (output)
C Local variables:
      REAL CONFIG(NCOORD)       !configuration
      REAL W                    !weight for single config
      REAL EBAR,WBAR            !ensmble average local energy and weight
      INTEGER IENSEM            !ensemble index
      INTEGER ICOORD            !coordinate index
      REAL NORM                 !normalization of weights
      REAL SHIFT(NCOORD)        !array containing drift vector
C Functions:
      REAL GAUSS                !gaussian random number   
      REAL ELOCAL               !local energy of the configuration
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      EBAR=0.                   !zero sums
      WBAR=0.
      DO 10 IENSEM=1,NENSEM                        !loop over ensemble
         DO 20 ICOORD=1,NCOORD
            CONFIG(ICOORD)=ENSMBL(ICOORD,IENSEM)   !get a configuration
20       CONTINUE
         CALL DRIFT(CONFIG,SHIFT)                  !calc shifts
         DO 30 ICOORD=1,NCOORD
            CONFIG(ICOORD)=CONFIG(ICOORD)+         !shift configuration
     +                     GAUSS(DSEED)*SQHBDT+SHIFT(ICOORD)
30       CONTINUE
C
         EPSILN=ELOCAL(CONFIG)                     !calculate energy
         WEIGHT(IENSEM)=WEIGHT(IENSEM)*EXP(-EPSILN*DT)  !calc weight
         EBAR=EBAR+WEIGHT(IENSEM)*EPSILN           !update sums
         WBAR=WBAR+WEIGHT(IENSEM)
C       
         DO 40 ICOORD=1,NCOORD
            ENSMBL(ICOORD,IENSEM)=CONFIG(ICOORD)   !save configuration
40       CONTINUE
10    CONTINUE
C
      EPSILN=EBAR/WBAR                    !weighted average energy
      NORM=NENSEM/WBAR
      DO 50 IENSEM=1,NENSEM               !renormalize weights
         WEIGHT(IENSEM)=NORM*WEIGHT(IENSEM)
50    CONTINUE
C
      RETURN
      END
c
      SUBROUTINE DRIFT(CONFIG,SHIFT)
C calculate the drift vector (SHIFT) for CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Input/Output variables:
      REAL CONFIG(NCOORD)        !configuration (input)
      REAL SHIFT(NCOORD)         !array containing drift vector (output)
C Local variables:
      INTEGER ICOORD              !coordinate index
      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 CHI,FDCHI,SDCHI,LAPCHI !atomic orbital
      REAL FEE,FDFEE,SDFEE,LAPFEE !elec-elec correlations
      REAL R                      !dummy variables
      REAL FACTA,FACTB,FACTE      !useful factors
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
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     get coordinates and radii
      CALL RADII(X1,X2,Y1,Y2,Z1,Z2,R1L,R2L,R1R,R2R,R12,CONFIG)
C
      CHI1=CHI(R1R)+CHI(R1L)           !pieces of the total wavefunction
      CHI2=CHI(R2R)+CHI(R2L)
      F=FEE(R12)
C
      FACTA=HBMDT*(FDCHI(R1L)/R1L+FDCHI(R1R)/R1R)/CHI1  !useful factors
      FACTB=HBMDT*(FDCHI(R1L)/R1L-FDCHI(R1R)/R1R)/CHI1
      FACTE=HBMDT*FDFEE(R12)/F/R12
C
      SHIFT(1)=FACTA*X1+FACTE*(X1-X2)  !shift for electron one
      SHIFT(2)=FACTA*Y1+FACTE*(Y1-Y2)
      SHIFT(3)=FACTA*Z1+FACTE*(Z1-Z2)+FACTB*S2
C
      FACTA=HBMDT*(FDCHI(R2L)/R2L+FDCHI(R2R)/R2R)/CHI2
      FACTB=HBMDT*(FDCHI(R2L)/R2L-FDCHI(R2R)/R2R)/CHI2
C
      SHIFT(4)=FACTA*X2-FACTE*(X1-X2)  !shift for electron two
      SHIFT(5)=FACTA*Y2-FACTE*(Y1-Y2)
      SHIFT(6)=FACTA*Z2-FACTE*(Z1-Z2)+FACTB*S2
C
      RETURN
      END
c
      REAL FUNCTION GAUSS(DSEED)
C returns a gaussian random number with zero mean and unit variance
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      INTEGER IGAUSS                 !sum index
      DOUBLE PRECISION DSEED         !random number seed
      REAL RANNOS                    !uniform random number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      GAUSS=0.                       !sum 12 uniform random numbers
      DO 10 IGAUSS=1,12
         GAUSS=GAUSS+RANNOS(DSEED)
10    CONTINUE
      GAUSS=GAUSS-6.                 !subtract six so that mean=0
      RETURN
      END
c
      SUBROUTINE INTENS(ENSMBL,WEIGHT,CONFIG)
C generate the ENSMBL at t=0 for PIMC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Output variables:
      REAL CONFIG(NCOORD)        !configuration
      REAL WEIGHT(MAXENS)        !weight of ensemble members
      REAL ENSMBL(NCOORD,MAXENS) !ensemble of configurations
C Local variables:
      INTEGER ISTEP              !step index
      INTEGER ICOORD             !coordinate index
      REAL W                     !weight for single config
      REAL ACCPT                 !acceptance ratio
      INTEGER IENSEM             !ensemble index
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CALL INTCFG(CONFIG,W)      !generate a single intial configuration
C
      DO 10 ISTEP=1,20                 !do 20 thermalization steps
         CALL METROP(CONFIG,W,ACCPT)   !using Metropolis algorithm
10    CONTINUE
C
      DO 30 ISTEP=1,10*NENSEM          !generate the ensemble
         CALL METROP(CONFIG,W,ACCPT)   !take a Metrop step
         IF (MOD(ISTEP,10) .EQ. 0) THEN
             IENSEM=ISTEP/10           !save every 10th config
             DO 20 ICOORD=1,NCOORD     
                ENSMBL(ICOORD,IENSEM)=CONFIG(ICOORD)
20           CONTINUE
             WEIGHT(IENSEM)=1.         !set all weights=1
         END IF
30    CONTINUE
      RETURN
      END