*dk   psurfac
      subroutine psurfac(v,rij)
c     written by b.j. archer
c     compute  the potential for a proton (#2), an electron (#3), and
c     an positron (#1)
      implicit none
      
      real*8 amass, bmass, cmass, udiatm, usys, vzero, autoev, 
     >       amutoau, angtoau, radtodeg, cmm1toau, rij, v
      integer mass, mass1, mass2, mass3 
      
      common /masses/ amass,bmass,cmass,udiatm,usys
      common / potzero / vzero
      common / convrsns / autoev,amutoau,angtoau,radtodeg,cmm1toau
      dimension rij(3), mass(2)
      logical first
      save first
      data first /.true./
      data mass/1000,1836/
c     on first time thru, check the masses for correctness
      if(first)then
         first = .false.
         mass1=amass*1.0e3+0.0001
         mass2=bmass*1.0e0+0.0001
         mass3=cmass*1.0e3+0.0001
         if(mass1.ne.mass(1).or.mass2.ne.mass(2)
     &        .or.mass3.ne.mass(1))then
          write(6,*)'*** error ***: masses not as expected '
          write(6,*)'for H+positron potential; expected the first four'
          write(6,*)'digits, in order, to be:',mass(1),mass(2),mass(1)
          write(6,*)'received, in order,    :',mass1,mass2,mass3
            stop 'psurfac'
         end if
      end if
c     compute the potential in hartrees, rij should be in bohr
c     rij(1) is the distance between particles 2 and 3, etc.
c     set the potential to a large number (2000eV) if the particles coalesce
      if(rij(1).eq.0.0.or.rij(2).eq.0.0)then
         v=-2000.0/autoev
      else if(rij(3).eq.0.0)then
         v=2000.0/autoev
      else
         v=-1.0/rij(1) - 1.0/rij(2) + 1.0/rij(3)
c     add an almost random factor to make sure the energies are
c     positive (0.57347952 hartree approx. 15.6053 ev)
         v=v+vzero
      end if
      return
      end