*dk h2ogr
      subroutine h2ogr(pot,r)
C     function potabc(rb,rs,cosgam)
C
C      H2O-GRUNDZUSTANDSFLAECHE VON SORBOIE + MURRELL
C
C      ABSTAENDE IN ANGSTROEM , ENERGIE IN EV
C
      implicit real*8 (a-h,o-z)
      dimension r(3)
      dimension pin(22),pout(22),vdiatm(3)
C
      data pin/1.d0,1.5421d0,4.6539d0,1.5720d0,-4.4355d0,
     *17.9914d0,-1.9761d0,
     *6.5279d0,-7.5781d0,4.9240d0,11.3559d0,22.3568d0,
     *23.1978d0,0.6787d0,
     *-11.4851d0,1.5255d0,-5.7924d0,59.6057d0,
     *-22.2496d0,20.8998d0,-18.39d0,
     *93.565d0/
C
      data pout/1.d0,1.875d0,-0.3691d0,3.1767d0,-2.4168d0,
     *3.1803d0,1.2577d0,
     *7.1737d0,-5.22d0,8.1229d0,6.2908d0,-2.2456d0,-1.7899d0,
     *0.3677d0,-6.7255d0,
     *0.8266d0,-3.1387d0,9.1025d0,-15.1988d0,8.0856d0,
     *-13.1075d0,15.906d0/
C
      data ain,aout/-0.9418d0,-1.7381d0/
      data gam1,gam2,gam3/4.5348d0,4.5348d0,2.0d0/
      data roh,rh2/0.9572d0,1.5139d0/
C     data rmasse/15.872893d0/
      real*8 angtoau, autoev
      data angtoau/0.529177249d0/, autoev/27.2113961d0/
C
C     return
C
C     entry pot(rb,rs,cosgam)
c
c  clip, if necessary
      pot=108.d0
      if(dmin1(r(2),r(3)).lt.0.71d0) return
c
c  The following conventions exist for the internuclear distances:
c     r1 = r(OH),  r2 = r(OH),  r3 = r(HH)
c  Convert from aph to sorbie-murrell coordinate convention
c  Convert distances from au to angstroms
      r1 = r(2)*angtoau
      r2 = r(3)*angtoau
      r3 = r(1)*angtoau
C
C     r1=rs
C     x=r1/(1.0d0+rmasse)
C     x=0.0d0
c     modifizierung
C     y=r1-x
C
C     r2=sqrt(rb**2+x**2-2*rb*x*cosgam)
C     r3=sqrt(rb**2+y**2+2*rb*y*cosgam)
C
C     OH-POTENTIALE :
C
      do 1 i=1,2
      if(i.eq.1)rr=r1
      if(i.eq.2)rr=r2
      if(rr.lt.0.6758804d0)goto 2
      dr=rr-0.9706d0
      vdiatm(i)=-4.63160d0*(1.0d0+4.534735d0*dr+5.03069d0*dr**2
     *+3.92551d0*dr**3)*exp(-4.534735d0*dr)
      goto 1
 2    vdiatm(i)=0.02338676d0/rr**12-2.573544d0
 1    continue
C
C     H2-POTENTIAL :
C
      dr=r3-0.7417d0
      if(r3.gt.1.67067d0)goto 3
      vdiatm(3)=-4.75010d0*(1.0d0+3.87591d0*dr+3.74086d0*dr**2
     *+3.25460d0*dr**3)*exp(-3.87591d0*dr) + 1.96d0
      goto 4
 3    vdiatm(3)=5.85684d0*(1.0d0+1.47818d0*dr+1.78015d0*dr**2)
     **exp(-3.9022d0*dr)
 4    continue
C
C     DER DREI-ATOMIGE ANTEIL :
C
      s1=r1-roh
      s2=r2-roh
      s3=r3-rh2
C
      s1q=s1*s1
      s2q=s2*s2
      s3q=s3*s3
      s1c=s1q*s1
      s2c=s2q*s2
      s3c=s3q*s3
      s1qu=s1q*s1q
      s2qu=s2q*s2q
      s3qu=s3q*s3q
C
      polyno=0.0d0
      do 100 i=1,22
C
      if(i.eq.1)pol=1.
      if(i.eq.2)pol=s1+s2
      if(i.eq.3)pol=s3
      if(i.eq.4)pol=s1q+s2q
      if(i.eq.5)pol=s3q
      if(i.eq.6)pol=s3*(s1+s2)
      if(i.eq.7)pol=s1*s2
      if(i.eq.8)pol=s1c+s2c
      if(i.eq.9)pol=s3c
      if(i.eq.10)pol=s1*s2q+s2*s1q
      if(i.eq.11)pol=s3q*(s1+s2)
      if(i.eq.12)pol=s3*(s1q+s2q)
      if(i.eq.13)pol=s1*s2*s3
      if(i.eq.14)pol=s1qu+s2qu
      if(i.eq.15)pol=s3qu
      if(i.eq.16)pol=s1c*s2+s2c*s1
      if(i.eq.17)pol=s1q*s2q
      if(i.eq.18)pol=s3*(s1c+s2c)
      if(i.eq.19)pol=s3q*(s1q+s2q)
      if(i.eq.20)pol=s3c*(s1+s2)
      if(i.eq.21)pol=s1*s2*s3q
      if(i.eq.22)pol=s3*(s1q*s2+s1*s2q)
C
      if(s3.lt.0.1568d0)p=pin(i)
      if(s3.gt.0.1568d0)p=pout(i)
 100  polyno=polyno+p*pol
C
      if(s3.lt.0.1568d0)polyno=ain *polyno
      if(s3.gt.0.1568d0)polyno=aout*polyno
C
      polyno=polyno*(1.0d0-tanh(0.5d0*gam1*s1))
     *             *(1.0d0-tanh(0.5d0*gam2*s2))
     *             *(1.0d0-tanh(0.5d0*gam3*s3))
C
C
      pot=vdiatm(1)+vdiatm(2)+vdiatm(3)+polyno
c  Convert energy from eV to au
      pot = (pot +10.08d0)/autoev
C
C
      return
      end
