        subroutine HHO_ESA_init()
c       V. Engel and R. Schinke, J.C.P vol. 88, 129 (1988) 
        return
        end
c==================================
        subroutine HHO_ESA_pes(xx,pe)
	implicit none
	real*8 xx(3),x(3),pot,pe
c  transform coordinates from HHO to HOH
        x(1)=xx(2)
        x(2)=xx(1)
        x(3)=xx(3)

	x(1)=x(1)*0.5292d0
	x(2)=x(2)*0.5292d0
	x(3)=x(3)*0.5292d0
	pe=pot(x)/27.2116d0
	x(1)=x(1)/0.5292d0
        x(2)=x(2)/0.5292d0
        x(3)=x(3)/0.5292d0

	return
	end
c---------
      FUNCTION POTI(r)                                                          
C
C     H2O - FITPOTENTIAL                                                        
C
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      DIMENSION R(3),RE(3),ALPHA(3),C(50),PHI(100)                             
C                                                                               
      DATA DOH,BETOH,ROH/4.621D0 , 2.294D0 , 0.971D0/                           
      data norder/50/
      data re/1.0000d0,1.0000d0,1.0000d0/
      data alpha/1.2500d0,0.7500d0,1.2500d0/ 
      data gamma/0.25d0/
      data c/ -.93562626d0,   -.33992392d1,   -.37923375d1,
     $	       .35229926d1,    .20508566d1,    .33290770d0,
     $         .64682997d0,   -.51447712d2,   -.49046646d1,
     $         .24435891d2,    .76337272d0,    .10582765d2,
     $         .73669551d1,    .64329848d2,    .89923222d1,
     $        -.89690335d2,    .64529244d2,   -.53894246d2,
     $        -.59093483d0,   -.12138697d2,    .65697698d1,
     $         .59787268d2,   -.29842680d2,    .30342123d2,
     $        -.29963985d2,    .24344415d2,   -.62823352d2,
     $         .18190455d3,    .10358976d2,   -.39393237d2,
     $         .39698652d1,    .12165069d2,    .51422343d1,
     $        -.47491777d1,   -.15323205d1,    .92191236d1,
     $        -.74335345d1,    .15738989d3,    .33299372d1,
     $        -.35652609d1,   -.15416913d3,   -.17472317d1,
     $         .52101774d2,    .16650884d3,   -.23987576d2,
     $        -.69990914d2,    .19884456d2,    .40429481d2,
     $        -.97198548d1,    .18539893d1/
c---------------------------
C                                                                               
      RMORSE(X)=DOH*(DEXP(-BETOH*(X-ROH))-1.)**2 - DOH                          
      XSI(X,Y)=(1.-DEXP(-gamma*X**2)) * (1.-DEXP(-gamma*Y**2))                  
C                                                                               
C     EINLESEN DER DATEN                                                        
C                                                                               
c      READ(22,100)NORDER,(RE(I),I=1,3),(ALPHA(I),I=1,3),gamma                    
c     *,(C(I),I=1,NORDER)                                                        
c 100  format(i3,/,3f12.4,/,3f12.4,/,f12.4,50(/,e16.8))                          
      write(6,*) 'Nach Einlesen der Daten norder:',norder
      DUMMY=RINTI(R)                                                            
      RETURN                                                                    
C                                                                               
      ENTRY POT(r)                                                              
C                                                                               
      RMAX=DMAX1(R(1),R(3))/0.5292                                              
      IF(RMAX.GE.5.5)GOTO 10                                                    
      CALL FUNC(R,PHI,NORDER,RE,ALPHA)                                          
      POT=0.0                                                                   
      DO 1 I=1,NORDER                                                           
        POT=POT+C(I)*PHI(I)
  1   continue
      POT=POT+RMORSE(R(1))*XSI(R(2),R(3))                                       
     *       +RMORSE(R(3))*XSI(R(2),R(1))                                       
      IF(RMAX.LE.3.5)RETURN                                                     
      VINT=RINT(R)                                                              
      COSCOS=DCOS(1.5707963*(RMAX-3.5))                                         
      GW=0.5+0.5*COSCOS                                                         
      FW=0.5-0.5*COSCOS                                                         
      POT=POT*GW+VINT*FW                                                        
      RETURN                                                                    
C                                                                               
 10   CONTINUE                                                                  
      POT=RINT(R)                                                               
      RETURN                                                                    
      END                                                                       
c                                                                               
      FUNCTION RINTI(R)                                                         
C                                                                               
C     H2O            INTERPOLATIONS POTENTIAL                                   
C                                                                               
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      REAL RLGDRE                                                               
C                                                                               
      DIMENSION R(3),GAMMA(5),A(5,5),B(5,6),P(5),X(6)                           
C                                                                               
      DATA D0,ALPHA0,R0/4.621,2.294,0.971/                                      
      DATA GAMMA/45.,75.,104.52,135.,180./                                      
      DATA B/-0.69629,-0.48884,-0.45792,-0.44409,-0.41761,                      
     *       1.40886,1.39153,1.39691,1.46819,1.54303,                           
     *       0.12042,0.14350,0.17865,0.19215,0.19151,                           
     *       0.19193,0.15280,0.13400,0.09861,0.05140,                           
     *       0.62717,1.38163,1.43055,1.36016,1.22214,                           
     *       0.31021,0.12674,0.18001,0.38321,0.51142/                           
C                                                                               
C                                                                               
C     BESTIMMUNG DER LEGENDRE ENTWICKLUNGS KOEFF.                               
C                                                                               
      DO 1 I=1,5                                                                
      DO 1 J=1,5                                                                
 1    A(I,J)=RLGDRE(J-1,COS(1.7453292E-2*SNGL(GAMMA(I))))                       
      CALL MAT56(A,5,B,6,DET)                                                   
      RETURN                                                                    
C                                                                               
      ENTRY RINT(R)                                                             
C                                                                               
      RMAX=DMAX1(R(1),R(3))                                                     
      RMAXA0=RMAX/0.5292 - 3.5                                                  
      RMIN=DMIN1(R(1),R(3))                                                     
      COSGAM=-(R(2)**2-R(1)**2-R(3)**2)/(2*R(1)*R(3))                           
      DO 2 I=1,5                                                                
 2    P(I)=RLGDRE(I-1,SNGL(COSGAM))                                             
      DO 3 J=1,6                                                                
      X(J)=0.0                                                                  
      DO 3 I=1,5                                                                
      X(J)=X(J)+B(I,J)*P(I)                                                     
 3    CONTINUE                                                                  
      D=D0+X(1)*DEXP(-X(2)*RMAXA0-X(3)*RMAXA0**2)                               
      ALPHA=ALPHA0+X(4)*DEXP(-X(5)*RMAXA0-X(6)*RMAXA0**2)                       
      RINT=D*(DEXP(-ALPHA*(RMIN-R0))-1)**2-D                                    
      RETURN                                                                    
      END                                                                       
c                                                                               
      FUNCTION RLGDRE(N,X)                                                      
C     CALCULATES LEGENDRE POLYNOMIALS USING RECURRENCE  ABRAM./STEGUN           
C     PAGE 782                                                                  
      RLGDRE=1                                                                  
      R2=RLGDRE                                                                 
      IF(N.EQ.0)RETURN                                                          
      RLGDRE=X                                                                  
      R1=RLGDRE                                                                 
      IF(N.EQ.1)RETURN                                                          
      DO 1 I=2,N                                                                
      RXR1=X*R1                                                                 
      RLGDRE=RXR1+RXR1-R2+(R2-RXR1)/I                                           
      R2=R1                                                                     
 1    R1=RLGDRE                                                                 
      RETURN                                                                    
      END                                                                       
c                                                                               
c                                                                               
c                                                                               
      SUBROUTINE MAT56(A,N,B,M,DETERM)                                          
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      DIMENSION INDEX(5,2),IPIVOT(5),A(5,5),B(5,6),PIVOT(5)                     
      EQUIVALENCE (IROW,JROW), (ICOLUM,JCOLUM), (AMAX, T, SWAP)      
C     INITIALIZATION                                                
   10 DETERM=1.0D0                                                 
   15 DO 20 J=1,N                                                 
   20 IPIVOT(J)=0                                                
   30 DO 550 I=1,N                                              
C     SEARCH FOR PIVOT ELEMENT                                 
   40 AMAX=0.0D0                                              
   45 DO 105 J=1,N                                           
   50 IF (IPIVOT(J)-1) 60, 105, 60                          
   60 DO 100 K=1,N                                         
   70 IF (IPIVOT(K)-1) 80, 100, 740                       
   80 IF( DABS(AMAX)- DABS(A(J,K)))85,100,100                                   
   85 IROW=J                                             
   90 ICOLUM=K                                          
   95 AMAX=A(J,K)                                      
  100 CONTINUE                                        
  105 CONTINUE                                       
  110 IPIVOT(ICOLUM)=IPIVOT(ICOLUM)+1               
C     INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL       
  130 IF (IROW-ICOLUM) 140, 260, 140                         
  140 DETERM=-DETERM                                        
  150 DO 200 L=1,N                                         
  160 SWAP=A(IROW,L)                                      
  170 A(IROW,L)=A(ICOLUM,L)                              
  200 A(ICOLUM,L)=SWAP                                  
  205 IF(M) 260, 260, 210                              
  210 DO 250 L=1, M                                   
  220 SWAP=B(IROW,L)                                 
  230 B(IROW,L)=B(ICOLUM,L)                         
  250 B(ICOLUM,L)=SWAP                             
  260 INDEX(I,1)=IROW                             
  270 INDEX(I,2)=ICOLUM                                                         
  310 PIVOT(I)=A(ICOLUM,ICOLUM)                  
  320 DETERM=DETERM*PIVOT(I)                    
C     DIVIDE PIVOT ROW BY PIVOT ELEMENT        
  330 A(ICOLUM,ICOLUM)=1.0D0                  
  340 DO 350 L=1,N                           
  350 A(ICOLUM,L)=A(ICOLUM,L)/PIVOT(I)      
  355 IF(M) 380, 380, 360                  
  360 DO 370 L=1,M                        
  370 B(ICOLUM,L)=B(ICOLUM,L)/PIVOT(I)   
C     REDUCE NON-PIVOT ROWS             
  380 DO 550 L1=1,N                    
  390 IF(L1-ICOLUM) 400, 550, 400     
  400 T=A(L1,ICOLUM)                 
  420 A(L1,ICOLUM)=0.0D0            
  430 DO 450 L=1,N                 
  450 A(L1,L)=A(L1,L)-A(ICOLUM,L)*T                    
  455 IF(M) 550, 550, 460                             
  460 DO 500 L=1,M                                   
  500 B(L1,L)=B(L1,L)-B(ICOLUM,L)*T                 
  550 CONTINUE                                     
C     INTERCHANGE COLUMNS                         
  600 DO 710 I=1,N                               
  610 L=N+1-I                                   
  620 IF (INDEX(L,1)-INDEX(L,2)) 630, 710, 630 
  630 JROW=INDEX(L,1)                         
  640 JCOLUM=INDEX(L,2)                      
  650 DO 705 K=1,N                          
  660 SWAP=A(K,JROW)                       
  670 A(K,JROW)=A(K,JCOLUM)               
  700 A(K,JCOLUM)=SWAP                   
  705 CONTINUE                          
  710 CONTINUE                         
  740 RETURN                          
      END                            
c                                                                               
c                                                                               
c                                                                               
      SUBROUTINE FUNC(R,PHI,NORDER,X,ALPHA)                                     
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      DIMENSION R(3),S(3),S2(3),S3(3),S4(3),S5(3)                               
      DIMENSION PHI(100),S6(3)                                                  
      DIMENSION ALPHA(3),X(3)                                                   
      COMMON BETA(3)                                                            
C                                                                               
      FAK=1.0                                                                   
      DO 1 I=1,3                                                                
      S(I)=R(I)-X(I)                                                            
      S2(I)=S(I)**2                                                             
      S3(I)=S2(I)*S(I)                                                          
      S4(I)=S2(I)*S2(I)                                                         
      S5(I)=S3(I)*S2(I)                                                         
      S6(I)=S3(I)**2                                                            
      FAK=FAK*(1.-DTANH(ALPHA(I)*S(I)))                                         
 1    CONTINUE                                                                  
C                                                                               
C                                                                               
      PHI(1)=1.                                                                 
C                                                                               
      PHI(2)=S(1)+S(3)                                                          
      PHI(3)=S(2)                                                               
C                                                                               
      PHI(4)=S2(1)+S2(3)                                                        
      PHI(5)=S2(2)                                                              
      PHI(6)=S(2)*(S(1)+S(3))                                                   
      PHI(7)=S(1)*S(3)                                                          
C                                                                               
      PHI(8)=S3(1)+S3(3)                                                        
      PHI(9)=S3(2)                                                              
      PHI(10)=S(1)*S2(3)+S2(1)*S(3)                                             
      PHI(11)=S2(2)*(S(1)+S(3))                                                 
      PHI(12)=S(2)*(S2(1)+S2(3))                                                
      PHI(13)=S(1)*S(2)*S(3)                                                    
C                                                                               
      PHI(14)=S4(1)+S4(3)                                                       
      PHI(15)=S4(2)                                                             
      PHI(16)=S3(1)*S(3)+S(1)*S3(3)                                             
      PHI(17)=S2(1)*S2(3)                                                       
      PHI(18)=S(2)*(S3(1)+S3(3))                                                
      PHI(19)=S2(2)*(S2(1)+S2(3))                                               
      PHI(20)=S3(2)*(S(1)+S(3))                                                 
      PHI(21)=S2(2)*S(1)*S(3)                                                   
      PHI(22)=S(2)*(S2(1)*S(3)+S(1)*S2(3))                                      
C                                                                               
C                                                                               
      PHI(23)=S5(1)+S5(3)                                                       
      PHI(24)=S4(1)*S(3)+S(1)*S4(3)                                             
      PHI(25)=S3(1)*S2(3)+S2(1)*S3(3)                                           
      PHI(26)=S(2)*(S4(1)+S4(3))                                                
      PHI(27)=S(2)*(S3(1)*S(3)+S(1)*S3(3))                                      
      PHI(28)=S(2)*S2(1)*S2(3)                                                  
      PHI(29)=S2(2)*(S3(1)+S3(3))                                               
      PHI(30)=S2(2)*(S2(1)*S(3)+S(1)*S2(3))                                     
      PHI(31)=S3(2)*(S2(1)+S2(3))                                               
      PHI(32)=S3(2)*S(1)*S(3)                                                   
      PHI(33)=S4(2)*(S(1)+S(3))                                                 
      PHI(34)=S5(2)                                                             
C                                                                               
      PHI(35)=S6(1)+S6(3)                                                       
      PHI(36)=S5(1)*S(3)+S(1)*S5(3)                                             
      PHI(37)=S4(1)*S2(3)+S2(1)*S4(3)                                           
      PHI(38)=S3(1)*S3(3)                                                       
      PHI(39)=S(2)*(S5(1)+S5(3))                                                
      PHI(40)=S(2)*(S4(1)*S(3)+S(1)*S4(3))                                      
      PHI(41)=S(2)*(S3(1)*S2(3)+S2(1)*S3(3))                                    
      PHI(42)=S2(2)*(S4(1)+S4(3))                                               
      PHI(43)=S2(2)*(S3(1)*S(3)+S(1)*S3(3))                                     
      PHI(44)=S2(1)*S2(2)*S2(3)                                                 
      PHI(45)=S3(2)*(S3(1)+S3(3))                                               
      PHI(46)=S3(2)*(S2(1)*S(3)+S(1)*S2(3))                                     
      PHI(47)=S4(2)*(S2(1)+S2(3))                                               
      PHI(48)=S4(2)*S(1)*S(3)                                                   
      PHI(49)=S5(2)*(S(1)+S(3))                                                 
      PHI(50)=S6(2)                                                             
C                                                                               
C                                                                               
      DO 2 I=1,NORDER                                                           
      PHI(I)=FAK*PHI(I)                                                         
 2    CONTINUE                                                                  
C                                                                               
      RETURN                                                                    
      END                                                                       
C
C
      FUNCTION SPHARM(L,M,X,Y)                                                  
C     CALCULATES SPAERICAL HARMONICS  EDMONDS 2.5.29                            
      COMPLEX SPHARM                                                            
      REAL*8  FAK,ASLEG                                                         
      IF(M.NE.0)GOTO 1                                                          
      SPHARM=0.2820947917*SQRT(L+L+1.)*RLGDRE(L,COS(X))                         
     **CMPLX(1.,0.)                                                             
      RETURN                                                                    
 1    CONTINUE                                                                  
      MABS=IABS(M)                                                              
      IF(MABS.LE.L)GOTO 3                                                       
          spharm=0.0                                                            
          return                                                                
 3    CONTINUE                                                                  
      IO=MABS+MABS                                                              
      FAK=1.D0                                                                  
      DO 2 I=1,IO                                                               
 2    FAK=FAK*SQRT(FLOAT(L-MABS+I))                                             
      YM=MABS*Y                                                                 
      SPHAR=0.2820947917*SQRT(L+L+1.)/SNGL(FAK)*SNGL(ASLEG(L,MABS,X))*          
     *(-1.)**MABS                                                               
      SPHARM=SPHAR*CMPLX(COS(YM),SIN(YM))                                       
      IF(M.GE.0)RETURN                                                          
      SPHARM=(-1)**M*CONJG(SPHARM)                                              
      RETURN                                                                    
      END                                                                       
C
C
C
      FUNCTION ASLEG(L,M,BETA)                                                  
      REAL*8 ASLEG,FAK,PLM1,PLM2                                                
      X=COS(BETA)                                                               
      MABS=IABS(M)                                                              
      IF(MABS.GT.0)GOTO 1                                                       
      ASLEG=RLGDRE(L,X)                                                         
      RETURN                                                                    
 1    FAK=1.D0                                                                  
      DO 2 I=1,MABS                                                             
 2    FAK=FAK*(MABS+I)                                                          
      ARG=0.5*BETA                                                              
      PLM2=FAK*(COS(ARG)*SIN(ARG))**MABS                                        
      IF(L.GT.MABS)GOTO 3                                                       
      ASLEG=PLM2                                                                
      GOTO 100                                                                  
 3    PLM1=(MABS+MABS+1)*X*PLM2                                                 
      IF(L.GT.(MABS+1))GOTO 4                                                   
      ASLEG=PLM1                                                                
      GOTO 100                                                                  
 4    IU=MABS+2                                                                 
      DO 5 I=IU,L                                                               
      ASLEG=((I+I-1)*X*PLM1-(I+MABS-1)*PLM2)/(I-MABS)                           
      PLM2=PLM1                                                                 
      PLM1=ASLEG                                                                
 5    CONTINUE                                                                  
 100  IF(M.GE.0)RETURN                                                          
      FAK=1.D0                                                                  
      IO=MABS+MABS                                                              
      DO 110 I=1,IO                                                             
 110  FAK=FAK*(L-MABS+I)                                                        
      ASLEG=(-1)**MABS*ASLEG/FAK                                                
      RETURN                                                                    
      END                                                                       
c**************END OF PROGRAM FOR A-STATE********************





