C RELHF.FOR D. Sober 12-Apr-89; 25-Nov-92 C Relativistic Hartree-Fock atomic form factors from C J.H. Hubbell and I. Overbo, J.Phys.Chem.Ref.Data 8, (1979) 69. C Table gives Z*F(x) vs. x = 20.60744 * q/mc. C Entries: C RHFSET(Z,FLAG,DUM) initializes log-log spline interpolation C for element Z and returns FLAG=1. if table for Z exists. C RELHF(Z,Q,F) returns F(Q), where C Q = momentum transfer in units of mc C F --> 1 as Q --> 0. C Subroutines called: C SPLINE, SPLINT (from Num. Rec.) -- included in this file C 25-Nov-92: Allow change of Z by returning F=-100. from RELHF C (previously stopped if Z changed) SUBROUTINE RHFSET(Z,FLAG,DUM) PARAMETER (NZ=3, XDQ=20.60744) c SAVE REAL X(100) REAL FTAB(100,NZ) REAL ZTAB(NZ) REAL QX(100) REAL Y(100) REAL Y2(100) REAL XL(100) REAL YL(100) REAL YL2(100) common /savedx/X,FTAB,ZTAB,QX,Y,Y2,XL,YL,YL2 vector zuse(1) DATA ZUSE/0./ DATA ZTAB/6.,50.,79./ DATA X/ 1 .0001,.01,.02,.03,.04,.05,.06,.07,.08,.09, .10,.11,.12,.13, 2 .14,.15,.16,.17,.18,.19, .20,.22,.24,.25,.26,.28,.30,.32,.34, 3 .35, .36,.38,.40,.42,.44,.45,.46,.48,.50,.55, .60,.65,.70,.80, 4 .90,1.,1.1,1.2,1.3,1.4, 1.5,1.6,1.7,1.8,1.9,2.,2.2,2.4,2.5,2.6, 5 2.8,3.0,3.3,3.5,3.6,3.9,4.0,4.2,4.6,5.0, 5.4,5.5,5.8,6.0,6.2, 6 6.6,7.0,7.4,8.0,9.0, 10.,11.,12.,14.,16.,18.,20.,22.,25.,28., 7 31.,35.,40.,45.,50.,70.,100.,1000.,1.0E6,1.0E9/ C Z=6 (C) C Z=50 (Sn) C Z=79 (Au) DATA FTAB/ 1 6.000,5.990,5.958,5.907,5.837,5.749,5.645,5.526,5.396,5.255, 2 5.107,4.952,4.794,4.633,4.472,4.311,4.153,3.998,3.847,3.701, 3 3.560,3.297,3.058,2.949,2.846,2.658,2.494,2.351,2.227,2.171, 4 2.120,2.028,1.948,1.880,1.821,1.794,1.770,1.725,1.685,1.603, 5 1.537,1.479,1.426,1.322,1.219,1.114,1.012,0.914,0.822,0.736, 6 .659,.588,.525,.468,.418,.373,.2986,.2402,.21587,.19430, 7 .1582,.1297,.0974,.080994,.07406,.05729,.05274,.04486,.03295, 8 .02465, .01877,.017577,.01451,.012822,.01137,.009031,.007256, 9 .005892,.004389,.002802, .001868,.001292,9.209E-4,5.014E-4, A 2.984E-4,1.877E-4,1.238E-4,8.498E-5,5.126E-5,3.273E-5, B 2.188E-5,1.353E-5,7.976E-6,5.005E-6,3.300E-6,8.745E-7, C 2.152E-7,3.665E-11,1.692E-20,1.713E-29, 1 50.,49.955,49.821,49.601,49.303,48.934,48.504,48.022,47.498, 2 46.942, 46.361,45.764,45.155,44.541,43.924,43.309,42.696,42.088, 3 41.486,40.891, 40.302,39.145,38.016,37.462,36.915,35.841,34.794, 4 33.775,32.786,32.303, 31.828,30.902,30.011,29.154,28.334,27.938, 5 27.551,26.805,26.096,24.482, 23.081,21.868,20.815,19.073,17.646, 6 16.384,15.201,14.062,12.962,11.913, 10.933,10.034,9.227,8.516, 7 7.897,7.367,6.5403,5.9401,5.7020,5.4959, 5.1491,4.8530,4.4497, 8 4.1920,4.0646,3.6880,3.5650,3.3246,2.8832,2.5090, 2.2011,2.1329, 9 1.9463,1.835,1.742,1.585,1.461,1.364,1.251,1.116, 1.011,.9179, A .8313,.6732,.5385,.4284,.3411,.2726,.1971,.1449, .1084,.07561, B .05012,.03453,.02458,8.1178E-3,2.472E-3,1.769E-6,4.07E-15, C 1.055E-23, 1 79.,78.957,78.826,78.609,78.311,77.936,77.491,76.981,76.414, 2 75.797, 75.135,74.437,73.706,72.950,72.173,71.380,70.575,69.761, 3 68.941,68.119, 67.296,65.657,64.039,63.241,62.452,60.902, 4 59.395,57.935,56.523,55.835, 55.160,53.846,52.581,51.363,50.191, 5 49.622,49.063,47.976,46.929,44.469, 42.207,40.110,38.153,34.581, 6 31.387,28.530,25.998,23.789,21.892,20.287, 18.943,17.821,16.880, 7 16.081,15.388,14.770,13.671,12.660,12.168,11.680, 10.729,9.826, 8 8.594,7.878,7.5538,6.7207,6.4890,6.0876,5.4687,5.0100, 9 4.6520,4.5745,4.3659,4.244,4.110,3.848,3.591,3.340,2.984,2.470, A 2.070,1.777,1.567,1.303,1.144,1.023,.9169,.8195,.6878,.5742, B .4788,.3769,.2842,.2150,.1662,.06907,.02576,4.408E-5,5.210E-13, C 6.519E-21/ FLAG=DUM DO 20 I=1,NZ IF(Z.EQ.ZTAB(I)) THEN IZ=I ZUSE(1)=Z WRITE(*,*) 'Intializing RELHF for Z = ',Z GO TO 30 END IF 20 CONTINUE FLAG=-1. WRITE(*,*) 'Element Z =',Z,' not found in RHFSET.' RETURN 30 FLAG=0. C Set up X and Y arrays for both linear and log-log interpolations: DO 40 J=1,100 QX(J)=X(J)/XDQ Y(J) =FTAB(J,IZ)/ZUSE(1) XL(J)=LOG(X(J)/XDQ) YL(J)=LOG(FTAB(J,IZ)/ZUSE(1)) 40 CONTINUE C Set up for cubic spline interpolations (linear and log-log): CALL SPLINE(QX,Y ,100,1.E30,1.E30,Y2) CALL SPLINE(XL,YL,100,1.E30,1.E30,YL2) C Value of Q above which log-log interp. is needed: C QXLIM=0. QXLIM=X(95)/XDQ END SUBROUTINE RELHF(Z,Q,F) REAL Z,Q,F PARAMETER (NZ=3, XDQ=20.60744) c SAVE REAL X(100) REAL FTAB(100,NZ) REAL ZTAB(NZ) REAL QX(100) REAL Y(100) REAL Y2(100) REAL XL(100) REAL YL(100) REAL YL2(100) common /savedx/X,FTAB,ZTAB,QX,Y,Y2,XL,YL,YL2 vector zuse(1) IF(Z.NE.ZUSE(1)) THEN IF(ZUSE(1).EQ.0.) THEN C On first call w/o initialization, return F=-99. as flag. WRITE(*,*) 'RELHF called before RHFSET' F=-99. RETURN ELSE WRITE(*,*) 'Change of Z in RELHF, Z =',Z F=-100. RETURN END IF END IF IF(Q.EQ.0.) THEN F=1. RETURN END IF C Spline interpolation using previously initialized derivatives: C Direct interp. (for max. speed), except at high end of table where C log-log interp. needed because of crude step size: IF(Q.GT.QXLIM) THEN QL=LOG(Q) CALL SPLINT(XL,YL,YL2,100,QL,FF) F=EXP(FF) ELSE CALL SPLINT(QX,Y,Y2,100,Q,F) END IF RETURN END C SPLINE.FOR (Numerical Recipes) SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2) PARAMETER (NMAX=100) REAL X(N),Y(N),Y2(N),U(NMAX) IF (YP1.GT..99E30) THEN Y2(1)=0. U(1)=0. ELSE Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE IF (YPN.GT..99E30) THEN QN=0. UN=0. ELSE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END C SPLINT.FOR (Numerical Recipes) SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y) REAL XA(N),YA(N),Y2A(N) KLO=1 KHI=N 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) IF (H.EQ.0.) PRINT *,'Bad XA input.' A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. RETURN END real function frelhf(Z,Q) real Q,Z real F real dum integer iflag vector zuse(1) if (Z.ne.zuse(1)) then call RHFSET(Z,iflag,dum) zuse(1)=Z endif call RELHF(Z,Q,F) if (F.lt.1) then frelhf = F else frelhf = 1 endif end real function ffdipole(Z,Q) real Q,Z real betaFF betaFF = 111*Z**(-1/3.) ffdipole=1/(1+(Q*betaFF)**2) end