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