FUNCTION GROUND(RHH,RLIHA,RLIHB)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     
C     THIS ROUTINE COMPUTES THE POTENTIAL ENERGY VALUES
C     OF THE GROUND ELECTRONIC STATE OF THE LIH2+ SYSTEM
C     AS EXPLAINED IN DETAIL IN
C
C       R. MARTINAZZO, G.F. TANTARDINI, E. BODO, F.A. GIANTURCO
C       J. CHEM. PHYS. 119 (2003) 11241
C
C     ON INPUT 
C       RHH, RLIHA, RLIHB ARE THE H-H AND LI-H INTERNUCLEAR DISTANCES
C       IN ATOMIC UNITS
C     ON OUTPUT
C       GROUND IS THE TOTAL ENERGY (IN ATOMIC UNITS) REFERRED TO THE 
C       THREE-BODY BREAK-UP ASYMPTOTE
C
C     CALL IT AFTER AN INITIALIZATION CALL TO GROUND_I
C     USE YOUR OWN 2-BODY POTENTIAL, IF DESIRED, BY CHANGING 
C     H2_KW AND XLIH_FCI FUNCTIONS
C
C     LIST OF SUBROUTINES AND FUNCTIONS:
C     
C      SUBROUTINE GROUND_I   INIZIALIZATION 
C      FUNCTION V3BODY       3-BODY POTENTIAL 
C      FUNCTION H2_KW        H2   GROUND-STATE POTENTIAL
C      FUNCTION XLIH_FCI     LIH+ GROUND-STATE POTENTIAL 
C      FUNCTION H2QUAD       HH QUADRUPOLE
C      FUNCTION H2POL        HH SPHERICAL, '3-BODY' POLARIZABILITY
C      FUNCTION H2POL2       HH ANYSOTROPIC POLARIZABILITY 
C     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT NONE
      DOUBLE PRECISION RHH,RLIHA,RLIHB
      DOUBLE PRECISION POT,V3BODY,H2_KW,XLIH_FCI
      DOUBLE PRECISION GROUND
      
      POT=V3BODY(RHH,RLIHA,RLIHB)
      POT=POT+H2_KW(RHH)+XLIH_FCI(RLIHA)+XLIH_FCI(RLIHB)      
      GROUND=POT

      END
      
      SUBROUTINE GROUND_I
      INTEGER MMAX,MLONG,MA
      PARAMETER (MMAX=10,MLONG=1,MA=143)
      INTEGER NDIM,NEXP1,NEXP2,NEXP3
      COMMON/DIM/ NDIM(MMAX)
      COMMON/EXP1/ NEXP1(MA-2-MLONG)
      COMMON/EXP2/ NEXP2(MA-2-MLONG)
      COMMON/EXP3/ NEXP3(MA-2-MLONG)
      INTEGER MAL,I,J,K,M,JJ,KK,IM
      NDIM(1)=0
      JJ=0
      DO M=2,MMAX
         KK=0
         DO I=0,M
            DO J=0,M
               
C     DO K=0,M
C     MODIFICATIONS FOR AB2 SYSTEMS ----> K.LE.J
C     FOR A3  SYSTEMS ----> K.LE.J.AND.I.LE.J
               
               DO K=0,J
                  IM=I+J+K
                  IF(IM.EQ.M)THEN
                     IF(I.NE.M.AND.J.NE.M.AND.K.NE.M)THEN
                        JJ=JJ+1
                        KK=KK+1
                        NDIM(M)=KK
                        NEXP1(JJ)=I
                        NEXP2(JJ)=J
                        NEXP3(JJ)=K
                     ENDIF
                  ENDIF
               ENDDO
            ENDDO
         ENDDO
         WRITE(*,*)'ORDER...',M
         WRITE(*,*)'NUMBER OF TERMS...',NDIM(M)
         WRITE(*,*)'EXPONENTS..'
         DO KK=1,NDIM(M)
            WRITE(*,*) NEXP1(JJ-NDIM(M)+KK),
     &           NEXP2(JJ-NDIM(M)+KK),
     &           NEXP3(JJ-NDIM(M)+KK)
         ENDDO
      ENDDO
      RETURN
      END

      FUNCTION V3BODY(R1,R2,R3)
      IMPLICIT NONE

      DOUBLE PRECISION R1,R2,R3,V3BODY

      INTEGER MMAX,MLONG,MA,NSW
      PARAMETER (MMAX=10,MLONG=1,MA=143,NSW=2)

      INTEGER NDIM,NEXP1,NEXP2,NEXP3
      COMMON/DIM/ NDIM(MMAX)
      COMMON/EXP1/ NEXP1(MA-2-MLONG)
      COMMON/EXP2/ NEXP2(MA-2-MLONG)
      COMMON/EXP3/ NEXP3(MA-2-MLONG)
      DOUBLE PRECISION A(MA)

      DATA A/     59.9041369487,
     &            -4.4190710855,
     &         -1793.3690652080,
     &           -92.7339536285,
     &           231.0047518887,
     &            23.6312294844,
     &         42386.6962109884,
     &         20887.4910336156,
     &          -101.8488429800,
     &         -4783.2820589398,
     &          5601.4299484618,
     &           647.9336802252,
     &         -2188.0964589478,
     &       -333391.4339867483,
     &        -50650.0836783211,
     &       -216249.6666120038,
     &        -10968.2992516742,
     &         14746.3468957760,
     &         64623.8750484771,
     &         -6585.4074148974,
     &       -121946.1513253672,
     &          2470.5168826501,
     &         39388.8512400949,
     &        674663.0178284654,
     &       -623209.8463494780,
     &       -625800.6376789842,
     &       5061315.2544497522,
     &       -499972.9861830364,
     &        738165.5689500275,
     &      -3062084.0673502581,
     &       -890975.0918916499,
     &         37100.2541750794,
     &        481176.7795094966,
     &         68335.2125605610,
     &        860916.1659883368,
     &       -182446.6634467566,
     &       -313469.8588537204,
     &      16296025.2245538831,
     &      24283257.8151341975,
     &       2226423.3841185221,
     &     -55670137.2963673249,
     &     -41861251.3156600073,
     &      12522078.8982328828,
     &     -10780744.6234461721,
     &      19920454.5483050309,
     &       3725852.5719358511,
     &       -888096.4455630139,
     &       5805044.5654820763,
     &      -1174676.0747827506,
     &       -152672.8060593364,
     &      -3582082.6546327011,
     &        252873.3936963925,
     &      -3291637.8143871599,
     &       1409200.5075722407,
     &       1356914.6186123572,
     &     -79156096.6236694306,
     &    -134289948.7898918688,
     &    -144168137.8314403296,
     &      18010795.4982191287,
     &     184494908.8457103074,
     &     201704765.4667312205,
     &     -99315107.5417414308,
     &      64030848.3463394046,
     &      62186512.0140229985,
     &     -69323914.7188201100,
     &     -18225109.7742339559,
     &       8767846.0950237475,
     &     -71368402.0258482397,
     &      20651551.6130857319,
     &       1535688.4530311411,
     &      12342388.8021743763,
     &       5825471.3843415780,
     &      -3734791.5687767509,
     &      10453889.2599439640,
     &       -704923.1248095465,
     &       7092463.0621291548,
     &      -5229921.5270726196,
     &      -3297509.8045103801,
     &     170414954.6685615480,
     &     442318060.4571669102,
     &     338310296.6504344344,
     &    -106295716.0848532170,
     &    -243177154.1358072460,
     &    -272160559.5405060053,
     &    -560001873.6517249346,
     &     345105141.7074456215,
     &    -179453069.6148003042,
     &    -205181566.6644286513,
     &     178625484.6556275785,
     &      80765193.6117105931,
     &     -34816331.7087118998,
     &     -22162585.9540813640,
     &     117201851.0785028338,
     &     -94646151.0332854241,
     &     -10872024.4385482874,
     &     134794077.4815447628,
     &     -23336856.3831176758,
     &      15174814.7355045378,
     &     -78927252.4897296727,
     &     -12871086.9000668079,
     &       4843936.1891660253,
     &     -11668605.0056117158,
     &        426732.0261371761,
     &      -8474699.3137837537,
     &       9738486.1953393687,
     &       4215023.2278999090,
     &      83709195.3309966177,
     &     -81955642.6349996924,
     &    -509178802.8625016212,
     &    -253164045.3674562573,
     &     145445180.5049311221,
     &    -261403144.5543840528,
     &    -214500767.8142659366,
     &     597722827.9000399113,
     &    -439796210.5928990245,
     &     195491463.8727613091,
     &     412858125.1965236664,
     &     989693212.4453485012,
     &     179744133.8334875405,
     &    -158573195.2249470055,
     &      47721187.7360739559,
     &    -309115588.9086966515,
     &    -865951061.6187733412,
     &     103449948.7291473746,
     &      25790340.4448547661,
     &      33253194.3017021790,
     &     415688777.1302601099,
     &     169521356.8767541051,
     &     -23916729.3615326211,
     &    -214129802.4606905282,
     &     -89171572.5754874349,
     &     -11048394.8646354862,
     &     122655134.3360205442,
     &      37156079.7304429114,
     &       1487230.4260974966,
     &      -1885306.1160342137,
     &       -502645.5527109082,
     &       5039172.7627611337,
     &      -7203742.8395571373,
     &      -2200156.4276824687,
     &             1.0277126653,
     &             1.4179477712,
     &             2.5479621480/

      DOUBLE PRECISION BETA1,BETA2,BETA3
      DOUBLE PRECISION RO1,RO2,RO3
      DOUBLE PRECISION RHO1,RHO2,RHO3,RHO_0,RHOU
      DOUBLE PRECISION YMOD,DYDA
      DOUBLE PRECISION RSCAT,COSANG
      DOUBLE PRECISION FACTOR,DENOM,VLONG,V3B
      DOUBLE PRECISION Q,DELALP_0,ALP_2,P2
      DOUBLE PRECISION H2POL,H2POL2,H2QUAD
      INTEGER MAL,I,J,K,M,JJ,KK,IM,N1,N2,N3

C     MAL=MA-3
C     AB2 SYSTEM

      MAL=MA-2-MLONG
     
      BETA1=A(MAL+1)
C     AB2 SYSTEM
      BETA2=A(MAL+2)
      BETA3=A(MAL+2)

C     LONG-RANGE PARAMETERS

      RHO_0=ABS(A(MA-MLONG+1))
      
      RO1=R1
      RO2=R2
      RO3=R3
      RHO1=RO1*DEXP(-BETA1*RO1)
      RHO2=RO2*DEXP(-BETA2*RO2)
      RHO3=RO3*DEXP(-BETA3*RO3)
      
      YMOD=0.0D0
      JJ=0

      DO M=2,MMAX
         DO KK=1,NDIM(M)
            JJ=JJ+1
            N1=NEXP1(JJ)
            N2=NEXP2(JJ)
            N3=NEXP3(JJ)
            IF(N2.EQ.N3)THEN
               DYDA=RHO1**N1*RHO2**N2*RHO3**N3
            ELSE
               DYDA=RHO1**N1*
     &              (RHO2**N2*RHO3**N3+
     &               RHO2**N3*RHO3**N2)
            ENDIF
            YMOD=YMOD+A(JJ)*DYDA
         ENDDO
      ENDDO

C     -----------------------------------------------------------
C     LONG-RANGE TERM EVALUATION
C

C     LONG-RANGE CONTRIBUTIONS FOR ARRANGEMENT 1, HH--LI
C     RSCAT IS THE SCATTERING COORDINATE IN THIS
C     ARRANGEMENT AND COSANG IS THE COSINE OF
C     THE JACOBI ANGLE
C
C                    LI
C                   /|
C                  / |
C              R2 /  | R3
C                /   |
C               /  THETA 
C             H------H
C                R1

      RSCAT=DSQRT(0.50D0*R2**2+0.50D0*R3**2+
     &     -0.25D0*R1**2)
      COSANG=(R2**2-R3**2)/R1/RSCAT/2.0D0

C     COMPUTE SWITCHING FUNCTION
      IF(RSCAT.LT.RHO_0)THEN
         FACTOR=0.0D0
         DENOM=1.0D0
         VLONG=0.0D0
         V3B=YMOD
      ELSE
         RHOU=RHO_0/RSCAT
         DENOM=1.0D0/(1.0D0-(RHOU)**NSW)
         FACTOR=EXP(1.0D0)*EXP(-DENOM)
C
C     COMPUTE HH QUADRUPOLE AND POLARIZABILITY FOR R=RHH
C     Q IS THE QZZ COMPONENT OF THE QUADRUPOLE TENSOR
C
C          ALP_0=1/3 *(ALPHA_XX+ALPHA_YY+ALPHA_ZZ)
C     AND 
C           POL2=(ALP_ZZ-ALPHA_0)
C
C           POL = ALP_0-2.0*ALPHA_H
C
C     P_0(X)=1, P_2(X)=1/2*(3*X**2-1)
C
         Q=H2QUAD(R1)
         DELALP_0=H2POL(R1)
         ALP_2= H2POL2(R1)
         P2=0.5*(3*COSANG**2-1.0)
         VLONG= Q*P2/RSCAT**3 +
     &        -DELALP_0/2.0D0/RSCAT**4            +
     &        -ALP_2*P2/2.0D0/RSCAT**4
         V3B=YMOD*(1.0D0-FACTOR) + VLONG*FACTOR
      ENDIF   
C     
C     END OF TOTAL LONG-RANGE POTENTIAL EVALUATION 
C     --------------------------------------------------------------

      V3BODY=V3B

      RETURN
      END


      FUNCTION H2QUAD(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION H2QUAD
C
C     HYDROGEN QUADRUPOLE (A.U.) 
C     AS A FUNCTION OF THE INTERNUCLEAR 
C     DISTANCE (A.U.)
C     STANDARD DEVIATION:  0.000325557608   
C
      PARAMETER(NA= 4)
      DIMENSION A(NA),FACT(2)
      DATA A /
     &     -0.0126271402,
     &     +0.341987071,
     &     -0.0249680149,
     &     +0.000251740527/
      DATA BETA/  0.0673276303/
      DATA FACT/  1.3899166,  3.30599025/
      DESP=DEXP(-BETA*R)
      RHO=R*DESP
      DUXPO=DEXP(FACT(1)*(R-FACT(2)))
      FACTOR=1.0D0/(1.0D0+DUXPO)
      YMOD=0.0D0
      DO M=1,NA
         YMOD=YMOD+A(M)*RHO**M
      ENDDO
      H2QUAD=YMOD*FACTOR
      RETURN
      END

      FUNCTION H2POL(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION H2POL
C
C     FULLCI H(Liu-augaug)
C     HYDROGEN 3-BODY
C     SPHERICAL POLARIZABILITY (A.U.)
C         DALPHA_0=ALPHA_0-2*ALPHA_H
C     AS A FUNCTION OF INTERNUCLEAR 
C     DISTANCE (A.U.)
C     STANDARD DEVIATION = 0.00741821229
C     
      PARAMETER(NA= 4)
      DIMENSION A(NA),FACT(2)
      DATA A /
     &     +0.0682048466,
     &     +1.69575813,
     &     -0.0379060701,
     &     +0.0149673329/
      DATA BETA/  0.0567659808/
      DATA FACT/  1.55657199,  3.19559801/
      DESP=DEXP(-BETA*R)
      RHO=R*DESP
      DESP1=DEXP(-  0.0724081327*R)
      DUXPO=DEXP(FACT(1)*(R-FACT(2)))
      FACTOR=1.0D0/(1.0D0+DUXPO)
      YMOD=0.0D0
      DO M=1,NA
         YMOD=YMOD+A(M)*RHO**M
      ENDDO
      YMOD=YMOD+ -7.4989518*DESP1
      H2POL=YMOD*FACTOR
      RETURN
      END

      FUNCTION H2POL2(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION H2POL2
C
C     FULLCI H(Liu-augaug)
C     HYDROGEN ANISOTROPIC 
C     POLARIZABILITY (A.U.)
C         ALPHA_2=ALPHA_ZZ-ALPHA_0
C     AS A FUNCTION OF INTERNUCLEAR
C     DISTANCE (A.U.)
C
C     STANDARD DEVIATION= 0.00863041785
C
      PARAMETER(NA= 5)
      DIMENSION A(NA),FACT(2)
      DATA A /
     &         0.0266251827,
     &         0.139081715,
     &         0.324064812,
     &        -0.0511107403,
     &         0.00279102422/
       DATA BETA/ -0.0729582362/
       DATA FACT/  1.59140975,  3.1416315/
       DESP=DEXP(-BETA*R)
       RHO=R*DESP
       DUXPO=DEXP(FACT(1)*(R-FACT(2)))
       FACTOR=1.0D0/(1.0D0+DUXPO)
       YMOD=0.0D0
       DO M=1,NA
             YMOD=YMOD+A(M)*RHO**M
       ENDDO
       H2POL2=YMOD*FACTOR
       RETURN
       END


C
C     ----  2-BODY POTENTIALS FOLLOW -------
C

      FUNCTION H2_KW(R)
C
C     FITTING OF THE KOLOS & WOLNIEWIC
C     PEC FOR THE H2 MOLECULE
C     FINAL STANDARD DEVIATION: 3.27320215E-05 A.U.
C     (i.e. 7.2 CM-1) 
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NA= 8)
      DIMENSION A(NA),FACT(2)
      DATA A /
     &        -0.142978719,
     &        -17.4496708,
     &         215.8658,
     &        -1763.38375,
     &         9338.55361,
     &        -30639.411,
     &         56147.2305,
     &        -43929.2685/
       DATA BETA/  1.19999843/
       DESP=DEXP(-BETA*R)
       RHO=R*DESP
       FACTOR=1.0D0
       YMOD=0.0D0
       DO M=1,NA
             YMOD=YMOD+A(M)*RHO**M
       ENDDO
       DESP1=DEXP(-  1.2*R)/R
       YMOD=YMOD+  0.637669117*DESP1
       H2_KW=YMOD*FACTOR
       RETURN
       END

      FUNCTION XLiH_FCI(R)
C
C     Fitting of the LiH+ FullCI curve 
C     with the Li(Y)+H(Liu-augaug)
C     basis-set. All values are in a.u..
C
C     STANDARD DEVIATION = 3.17247272E-05 	
C     (i.e. 6.96 cm^-1).
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NA= 6)
      DIMENSION A(NA),FACT(2)
      DATA A /
     &        -0.436998324,
     &        -10.9156475,
     &         55.0379023,
     &        -437.23156,
     &         1443.14556,
     &        -2340.82341/
       DATA BETA/  1.00000006/
       DATA FACT/  5.85874466,  1./
       DESP=DEXP(-BETA*R)
       RHO=R*DESP
       FACTOR=1.0D0
       factor1=0.0d0
       rhou=  5.85874466/R
       if(rhou.lt.1.0d0)then
         factor1=dexp(-1.0d0/(1.0d0-rhou** 2))
         factor1=factor1*dexp(1.0d0)
       endif
       YMOD=0.0D0
       DO M=1,NA
             YMOD=YMOD+A(M)*RHO**M
       ENDDO
       DESP1=DEXP(-  1.*R)/R
       YMOD=YMOD+  16.1279348*DESP1
       VLONG=-  4.5/2.0D0/R**4
       YMOD=YMOD*(1.0d0-FACTOR1)+FACTOR1*VLONG
       XLIH_FCI=YMOD*FACTOR
       RETURN
       END