      FUNCTION EXCITED(RHH,RLIHA,RLIHB)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     
C     THIS ROUTINE COMPUTES THE POTENTIAL ENERGY VALUES
C     OF THE FIRST EXCITED 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       EXCITED 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 EXCITED_I
C     USE YOUR OWN 2-BODY POTENTIAL, IF DESIRED, BY CHANGING 
C     H2P, H2P, HLIP FUNCTIONS
C
C     LIST OF SUBROUTINES AND FUNCTIONS:
C     
C      SUBROUTINE EXCITED_I  INIZIALIZATION 
C      FUNCTION EX3BODY      3-BODY POTENTIAL 
C      FUNCTION H2P          H2+ GROUND-STATE POTENTIAL
C      FUNCTION HLI          LIH GROUND-STATE POTENTIAL
C      FUNCTION HLIP         LIH+ POTENTIAL (A STATE)
C      FUNCTION HLIDIP       LIH DIPOLE
C      FUNCTION HLIQUAD      LIH QUADRUPOLE
C      FUNCTION HLIPOL       LIH SPHERICAL, '3-BODY' POLARIZABILITY
C      FUNCTION HLIPOL2      LIH ANYSOTROPIC POLARIZABILITY 
C     
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT NONE
      DOUBLE PRECISION RHH,RLIHA,RLIHB,EXCITED
      DOUBLE PRECISION POT,EX3BODY,H2P,HLI,HLIP
      DOUBLE PRECISION RMIN,RMAX

      POT=EX3BODY(RHH,RLIHA,RLIHB)
      RMIN=MIN(RLIHA,RLIHB)
      RMAX=MAX(RLIHA,RLIHB)
      POT=POT + H2P(RHH)+ HLI(RMIN)+ HLIP(RMAX)
      EXCITED=POT

      RETURN
      STOP
      END

      SUBROUTINE EXCITED_I
      IMPLICIT NONE
      DOUBLE PRECISION HLIPOL,HLIPOL2,HLIQUAD
      INTEGER MMAX,MA,MLONG,NSW
      PARAMETER (MMAX=12,MA=230,MLONG=1,NSW=4)
      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
      STOP
      END
 

      FUNCTION EX3BODY(R1,R2,R3)
      IMPLICIT NONE
      DOUBLE PRECISION R1,R2,R3,EX3BODY
      INTEGER MMAX,MA,MLONG,NSW
      PARAMETER (MMAX=12,MA=230,MLONG=1,NSW=4)
      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/      7.0139563146,
     &             3.0003108967,
     &           -76.4867842937,
     &           581.0183707108,
     &          -447.5847281554,
     &           -99.4864290612,
     &           982.7398204780,
     &           362.5648698799,
     &         -3393.2414511771,
     &          5178.1173784156,
     &         -7863.7073601031,
     &          6560.3607376437,
     &          1113.8316789772,
     &         -5017.3030972528,
     &          -440.8591101941,
     &         25489.4179184216,
     &         11173.7792402553,
     &        -35732.6711370254,
     &         30244.5599087393,
     &        -48056.3359768260,
     &         77877.3998465644,
     &        -71808.6974729424,
     &          5764.8295676300,
     &          1006.7053378493,
     &         28332.5137592468,
     &         -8092.2103234528,
     &        -79498.1299204307,
     &        -31293.4097147527,
     &        164370.3230362090,
     &       -191847.1898188904,
     &        -85926.7361277200,
     &        260275.4135666166,
     &       -133673.0490531631,
     &        263387.6756272943,
     &       -632639.9064772174,
     &        603041.5077466281,
     &       -236922.6937759927,
     &          6700.8245146343,
     &       -113425.9637954208,
     &         60455.4406802392,
     &        378520.5104870343,
     &         86984.0196323107,
     &         55343.8072617274,
     &       -503235.4170471951,
     &        285008.8945165113,
     &        516410.1142055998,
     &      -1024788.4006496501,
     &       1354991.4595010185,
     &       -235933.4600355194,
     &       -868190.4855862955,
     &        603123.2438464863,
     &      -1252224.6800955555,
     &       3314794.4823525730,
     &      -3330758.8231016411,
     &       2355595.5164276436,
     &       -207196.8778448466,
     &        101681.1108576819,
     &        257460.6576209367,
     &       -201212.0670082429,
     &       -580128.3865388532,
     &       -126806.8597396127,
     &         54164.8740416092,
     &       1010250.0667563326,
     &      -2188435.0117557826,
     &        797679.9871993139,
     &      -1793697.7770774949,
     &       2611739.3747180039,
     &       -883251.7685729633,
     &      -1349145.2436505784,
     &       3009615.9249107814,
     &      -7052322.0718244938,
     &       3874587.9290700238,
     &       1225199.8923960077,
     &      -3507124.1024165265,
     &       5796231.3624979686,
     &      -8891846.5984294061,
     &      10832850.4398638811,
     &     -12425853.1764493100,
     &        301736.6636589296,
     &       -371127.7362507661,
     &       -321645.7644776623,
     &        374599.8369838075,
     &       1973036.8288541059,
     &       -225429.4068309222,
     &        581127.4474386526,
     &       -459020.9465607482,
     &      -1293190.5723435762,
     &       1565452.0461066423,
     &      -1254345.0027045240,
     &       2769830.0552521981,
     &      -4183276.3368961909,
     &       9865817.5040407330,
     &      -5687710.6781170880,
     &       6331461.7343169162,
     &      -6199060.3498003324,
     &       2366805.4481070125,
     &       -295291.9032663772,
     &      -4626557.7481302628,
     &      21978300.2186226137,
     &     -14309379.2452538181,
     &         96083.8516814366,
     &      14169688.6588801872,
     &     -20227745.2251465656,
     &       4905125.0786394291,
     &     -17851102.8897421658,
     &      38924717.4709330797,
     &       -576834.7877907237,
     &        -46335.9803832678,
     &        504037.2236824127,
     &        202416.2680335442,
     &       -401805.4777432226,
     &      -1465800.8280503147,
     &        738207.4260055233,
     &       -986097.7081596345,
     &        840428.3375779393,
     &        991264.5669337681,
     &      -5995665.7003214313,
     &       3238025.1562320245,
     &       -198540.6010630105,
     &      -2406737.5763311395,
     &       4214995.0096759945,
     &      -2530526.3486393159,
     &        276406.6675839963,
     &      -4399887.9279605504,
     &       6257808.5376846725,
     &     -38630207.8119526505,
     &      29587282.2978549488,
     &     -17989117.3686880358,
     &      10206444.9309776053,
     &      -7956727.4088536510,
     &      11935639.2603963744,
     &        272587.5750498392,
     &     -33205003.3066417538,
     &      17814420.1847252585,
     &       1717949.1537269219,
     &     -29524642.9344762526,
     &      41664281.4065521657,
     &      30649831.9779503755,
     &       6411040.1861357680,
     &     -72601162.1412622780,
     &        301573.8474337890,
     &       -151845.9726148361,
     &       -312104.3259274737,
     &        -41842.7503088380,
     &        232165.1746365801,
     &       2329864.8484232700,
     &       -431076.9872573870,
     &        -81290.3065144457,
     &        668919.9198792472,
     &       -695227.8774063820,
     &       -391985.4736881022,
     &       -617857.7470953518,
     &        868434.0566439229,
     &      -1795080.1182449430,
     &       2099682.7952767937,
     &      -2570819.2477941490,
     &      39547802.8107521236,
     &     -26430084.7782766521,
     &      14040722.8520985711,
     &      -2630006.2582660136,
     &      -2890426.7765964856,
     &      -4285618.4932269054,
     &      -1483471.0828852160,
     &       9488176.4738190025,
     &      -7094951.3043547636,
     &      89999375.1135660559,
     &     -53602764.0800236166,
     &      16091618.2000474222,
     &      -6629563.5433165785,
     &      -1253521.0871043403,
     &     -11254990.7965670228,
     &       5098017.7156343600,
     &      34072707.4960278347,
     &      -2872237.4304915103,
     &     -11445974.2646609601,
     &      25413934.4136321582,
     &     -42223432.0668093711,
     &     -68299818.6303414106,
     &      19227087.4643655121,
     &      74483220.7264803648,
     &       -179309.5475248325,
     &         17617.4766226398,
     &         53230.0044518733,
     &         78078.6615530894,
     &         -7811.8766089654,
     &        -55943.5555708998,
     &       -267671.7329121784,
     &       -117273.5346655442,
     &        160907.6128860498,
     &       -260296.6368125519,
     &        243850.4724256268,
     &         51209.2539915440,
     &     -21857101.1706669591,
     &      18320151.9379566535,
     &      -9529829.7924672943,
     &       3690041.8095493489,
     &      -1307964.2923558105,
     &        749871.0760208279,
     &       1844403.7695033511,
     &      -6647893.7860741401,
     &       2730696.9130423600,
     &       -112201.7731048101,
     &        689349.4466145103,
     &    -124924487.2901448756,
     &      99219617.4259525836,
     &     -41801056.6200024039,
     &       6218141.7051764801,
     &       1148692.0085307148,
     &       8811748.7360878177,
     &     -19204502.2225882821,
     &       7627010.5606003199,
     &       1193908.8264275370,
     &     -96724132.1153392941,
     &      64398757.3842907548,
     &     -15602214.6525493693,
     &        859044.0599557476,
     &       8829806.1955996510,
     &      -4666527.9020765796,
     &       -379501.8213393019,
     &     -24132122.5725693628,
     &      -2759943.9894857286,
     &      10162478.5972349979,
     &      -4359214.8902296340,
     &      15410672.6411291305,
     &      43276836.9157600775,
     &     -19009630.1279357336,
     &     -32327109.3021873273,
     &             0.8545339729,
     &             0.5385533051,
     &             8.6596046408 /

      DOUBLE PRECISION BETA1,BETA2,BETA3
      DOUBLE PRECISION RO1,RO2,RO3
      DOUBLE PRECISION RHO1,RHO2,RHO3,RHO_0,RHOU
      DOUBLE PRECISION YMOD,DYDA
      DOUBLE PRECISION RA,RB,RC
      DOUBLE PRECISION RSCAT,COSANG
      DOUBLE PRECISION FACTOR,DENOM,VLONG,V3B
      DOUBLE PRECISION DIP,Q,DELALP_0,ALP_2,P1,P2
      DOUBLE PRECISION HLIDIP,HLIPOL,HLIPOL2,HLIQUAD
      INTEGER MAL,I,J,K,M,JJ,KK,IM,N1,N2,N3

      MAL=MA-MLONG-2
      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))

      RHO1=R1*DEXP(-BETA1*R1)
      RHO2=R2*DEXP(-BETA2*R2)
      RHO3=R3*DEXP(-BETA3*R3)
      
      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

      RA=MIN(R2,R3)
      RB=MAX(R2,R3)
      RC=R1

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                    H
C                   /|
C                  / |
C              RB /  | RC
C                /   |
C               /  THETA 
C            Li------H
C                RA

      RSCAT=DSQRT(0.50D0*RB**2 + 0.50D0*RC**2+
     &     -0.25D0*RA**2 )
      COSANG=(RB**2-RC**2)/RA/RSCAT/2.0D0

C     COMPUTE SWITCHING FUNCTION
      IF(RSCAT.LE.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 LIH DIPOLE, QUADRUPOLE AND POLARIZABILITY FOR R=RLIH
C     DIP IS THE DIPOLE MOMENT
C     Q   IS THE QZZ COMPONENT OF THE QUADRUPOLE TENSOR
C            WITH RESPECT TO THE GEOM CENTER OF LIH BOND
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_1(X)=X  P_2(X)=1/2*(3*X**2-1)
C
         DIP=HLIDIP(RA)
         Q=HLIQUAD(RA)
         DELALP_0=HLIPOL(RA)
         ALP_2= HLIPOL2(RA)
         P1=COSANG
         P2=0.5*(3*COSANG**2-1.0)
         VLONG= - DIP*P1/RSCAT**2                  +
     &        + 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     --------------------------------------------------------------

      EX3BODY=V3B

      RETURN
      END

   
      FUNCTION HLIDIP(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NA= 5)
      DIMENSION A(NA),FACT(2)
C     LIH FCI DIPOLE CURVE IN A.U.
C     STANDARD DEVIATION..  0.0115522485
C     NOTE THE NEED OF 'FACTOR' IN ORDER
C     TO AVOID UNPHYSICAL BEHAVIOUR
C     WITH THE SIMPLE AG ANSATZ
      DATA A /
     &         3.81107729,
     &        -3.07621218,
     &         1.20552107,
     &        -0.213831318,
     &         0.0142348131/
      DATA BETA/  0.0307255411/
      DATA FACT/  1.664664,  6.68701348/
      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
      HLIDIP=YMOD*FACTOR
      RETURN
      END

      FUNCTION HLIQUAD(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NA= 6)
      DIMENSION A(NA),FACT(2)
C     LIH QUADRUPOLE FUNCTION (IN A.U.)
C     WITH RESPECT THE GEOMETRIC CENTER
C     OF THE MOLECULE. FROM  FULLCI
C     CALCUTATIONS 
C     STANDARD DEVIATION  = 0.000592303467 A.U.
      DATA A /
     &        -4.05169074,
     &         7.010268,
     &        -4.53005171,
     &         1.53830649,
     &        -0.26118265,
     &         0.0176316218/
      DATA BETA/  0.0381326063/
      DATA FACT/  1.72807586,  6.54214083/
      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
      HLIQUAD=YMOD*FACTOR
      RETURN
      END

      FUNCTION HLIPOL(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NA= 5)
C     SPHERICAL 3-BODY POLARIZABILITY
C     STANDARD DEVIATION 0.427422896
      DIMENSION A(NA),FACT(2)
      DATA A /
     &        -137.958769,
     &         68.7249887,
     &        -14.5416019,
     &         1.24211167,
     &        -0.00608506371/
      DATA BETA/  0.00608598433/
      DATA FACT/  1.64572959,  6.74611063/
      DESP=DEXP(-BETA*R)
      RHO=R*DESP
      DESP1=DEXP(- -0.0040014362*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+ -52.7688091*DESP1
      HLIPOL=YMOD*FACTOR
      RETURN
      END
      
      FUNCTION HLIPOL2(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NA= 5)
      DIMENSION A(NA),FACT(2)
      DATA A /
     &        -3.10727044,
     &        -4.16298439,
     &         3.20124828,
     &        -0.724437507,
     &         0.0614349836/
      DATA BETA/ -0.0274272462/
      DATA FACT/  1.77848435,  6.69233524/
      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
      HLIPOL2=YMOD*FACTOR
      RETURN
      END


      FUNCTION H2P(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     Fitting of the exact H2+
C     Wind's ab-initio values.
C     We considered only data point
C     which are 0.05 a.u. above dissociation
C     (381 out of 400).
C     STANDARD DEVIATION..  3.36611192E-05
C     (i.e. 7.39 cm^-1)
C
      PARAMETER(NA= 8)
      DIMENSION A(NA),FACT(2)
      DATA A /
     &        -0.433932717,
     &        -1.46625756,
     &         26.1872609,
     &        -227.903359,
     &         1135.07033,
     &        -3281.68774,
     &         5086.2199,
     &        -3273.43527/
       DATA BETA/  0.900001396/
       DATA FACT/  10.9526825,  1./
       DESP=DEXP(-BETA*R)
       RHO=R*DESP
       FACTOR=1.0D0
       factor1=0.0d0
       rhou=  10.9526825/R
       if(rhou.lt.1.0d0)then
         factor1=dexp(-1.0d0/(1.0d0-rhou** 8))
         factor1=factor1*dexp(1.0d0)
       endif
       YMOD=0.0D0
       DO M=1,NA
             YMOD=YMOD+A(M)*RHO**M
       ENDDO
       DESP1=DEXP(-  0.899999999*R)/R
       YMOD=YMOD+  0.616429616*DESP1
       VLONG=-  4.5/2.0D0/R**4
       YMOD=YMOD*(1.0d0-FACTOR1)+FACTOR1*VLONG
       H2P=YMOD*FACTOR
       RETURN
       END

      FUNCTION HLI(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NA= 8)
      DIMENSION A(NA),FACT(2)
C
C     Standard deviation 0.000253105716 a.u.
C     (55.6 cm^-1)
C
      DATA A /
     &        -10.7239499,
     &        -117.330715,
     &        -1700.19362,
     &         26575.4401,
     &        -273771.612,
     &         1500872.94,
     &        -4367963.76,
     &         5024043.4/
       DATA BETA/  1.2/
       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+  216.039947*DESP1
       HLI=YMOD*FACTOR
       RETURN
       END
      
      FUNCTION HLIP(R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NA= 6)
      DIMENSION A(NA),FACT(2)
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..  8.97343356E-05
C     (i.e. 19.69 cm^-1)
C 
      DATA A /
     &        -0.161505456,
     &         0.151370422,
     &         1.09938176,
     &        -3.30073329,
     &         5.16501195,
     &        -2.98494189/
       DATA BETA/  0.479976824/
       DATA FACT/  8.81582567,  1./
       DESP=DEXP(-BETA*R)
       RHO=R*DESP
       FACTOR=1.0D0
       factor1=0.0d0
       rhou=  8.81582567/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(-  0.480098406*R)/R
       YMOD=YMOD+  0.976997582*DESP1
       VLONG=-  165./2.0D0/R**4
       YMOD=YMOD*(1.0d0-FACTOR1)+FACTOR1*VLONG
       HLIP=YMOD*FACTOR
       RETURN
       END


