c**********************************************************************
c** R.J. Le Roy's program 'RKR1' for calculating RKR turning points in
c either simple first order or first-order Kaiser approximation.
c** 'G(v)' and 'B(v)' may be generated using conventional (Dunham)
c polynomials in (v+1/2), using Near-Dissociation expansions (NDE's),
c or using Tellinguisen's MXS 'mixed' Dunha,-at-low-v and NDE-at-high-v
c functional representations.
c** If desired (see manual), this program will also extrapolate up the
c inner wall past a chosen point with an exponential fitted to the 3
c preceeding points, & adjust 'RT(outer)' accordingly.
c******************** Version of 8 January 2016 **********************
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c COPYRIGHT 2003-2016 by Robert J. Le Roy +
c Dept. of Chemistry, Univ. of Waterloo, Waterloo, Ontario, Canada +
c This software may not be sold or any other commercial use made +
c of it without the explicit written permission of the author. +
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INTEGER MXTP,MXLEV,MXDUN
PARAMETER (MXLEV=500,MXTP=1001,MXDUN=25)
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c** NOTE: Current dimensioning allows calculation of pairs of turning
c points for up to MXLEV 500 vibrational levels with a total of
c MXTP=2*MXLEV+1 turning points, and for Dunham and NDE
c polynomial expansions of defined by up to MXDUN= 25 parameters.
c These limits may may be changed by the user.
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CHARACTER*78 TITLE
CHARACTER*6 PTYPE(6)
CHARACTER*2 NAME1,NAME2,CCDC(0:1)
INTEGER IAN1,IAN2,IMN1,IMN2,GEL1,GEL2,GNS1,GNS2,CHARGE,
1 KORDR,MORSE,NGP,Kaiser,IBIS,NBIS,IDIV,NDIV,NSV,NTP,
2 IVIB,NVIB,IVB,I634,I638,I,J,K,L,NN,I1,I2,I1P1,I1P2,I640,I644,
3 IR1,IDR
c
REAL*8 CU,SQCU,MASS1,MASS2,ZMASE,ABUND1,ABUND2,TOLER,ZMU,Req,Be,
1 AE,De,WE,WEXE,BETA,V00,DV0,Sw,SwLR,VEXT,F,FF,FB,G,GG,GB,VST,VFN,
2 RANGE,TSTF,TSTG,TSTFB,TSTGB,VUP,GUP,BUP,DGUP,RMIN,RMAX,RMIN2,
3 RMIN3,VRAT,CEXT,DCEXT,ADCEXT,BDCEXT,EXP1,EXP2,EXP3,FUN,
4 DFUN,AEXT,BEXT,REXT,DR1,DEI,DEIB,DRI,DRIB,D1,D2,TT,HEL,
5 V1(9),V2(9),DV(9),V(MXLEV),RT(-4:MXTP),
4 ET(-4:MXTP),VX(MXDUN),BV(MXDUN),GV(MXDUN),DGDV(MXDUN),WW(16)
c** Common block for Dunham & MXS function parameters
INTEGER LMAXGv,LMAXBv,NDEGv,NDEBv
REAL*8 VS,DVS,Y00,YL0(0:MXDUN),YL1(0:MXDUN)
COMMON /DUNPRM/Y00,VS,DVS,YL0,YL1,LMAXGv,LMAXBv,NDEGv,NDEBv
c** Common block for NDE function parameters
INTEGER NLR,ITYPE,ITYPB,IZP0,IZQ0,IZP1,IZQ1,NP0,NQ0,NP1,NQ1
REAL*8 VD,DLIM,XCN0,XCN1,P0(MXDUN),Q0(MXDUN),P1(MXDUN),Q1(MXDUN)
COMMON /NDEPRM/VD,DLIM,XCN0,XCN1,P0,Q0,P1,Q1,NLR,ITYPE,ITYPB,
1 IZP0,IZQ0,IZP1,IZQ1,NP0,NQ0,NP1,NQ1
c** Common block for quadrature weights & points
REAL*8 XG(32),WG(32),X2(16),W2(16)
COMMON /GWGHT/XG,WG,X2,W2
c
DATA CCDC/'Gv','Bv'/
DATA PTYPE/' OUTER',' INNER','Expone',' Pade ',' Pade ','ntial '/
DATA ZMASE /5.4857990946D-04/ !! 2010 physical constants d:mohr12
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c** Quadratures performed using NGP point Gaussian integrations.
NGP= 16
CALL WGHT(NGP)
c** Quadrature convergence criterion is TOLER
TOLER= 1.d-10
c** Use up to NBIS interval bisections in the NGP-point gaussian
c integration for the "f" and "g" integrals (typically set NBIS=3-5).
NBIS= 5
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c** Begin by reading in the (integer) atomic numbers and mass numbers
c defining the effective reduced mass of the system considered.
c** IAN1 & IAM2, and IMN1 & IMN2 are, respectively, the atomic numbers
c and the mass numbers identifying the atoms forming the molecule.
c Their masses are extracted from data subroutine MASSES and used
c to generate the the reduced mass ZMU.
c** If IMN1 or IMN2 lie outside the range of mass numbers for normal
c stable isotopes of that species, subroutine MASSES returns the
c average atomic mass based on the natural isotope abundance.
c** If the read-in value of IAN1 and/or IAN2 is .LE.0, then instead of
c using the MASS table, read an actual particle mass for it/them.
c** CHARGE (integer) is the charge on the molecule (=0 for neutral). If
c (CHARGE.ne.0) generate & use Watson's charge-adjusted reduced mass.
c** NDEGv(s) & NDEBv(s) specify whether Gv & Bv are represented by:
c (a) pure Dunham expansions when NDEXv(s) = 0
c (b) pure NDE expressions when NDEXv(s) = 1
c (c) MXS mixed NDE/Dunham expressions when NDEXv(s) > 1
c!! SPECIAL CASE: if have NO v-dependent rotational data and wish to
c define inner potential wall by Morse function, set: NDEBv =-1 !!
c!! NOTE: Program does NOT allow cases with NDEBv > NDEGv
c----------------------------------------------------------------------
2 READ(5,*,END=99) IAN1, IMN1, IAN2, IMN2, CHARGE, NDEGv, NDEBv
c----------------------------------------------------------------------
Y00= 0.d0
V00= -0.5d0
MORSE= 0
IF(NDEBv.GT.NDEGv) THEN
WRITE(6,598) NDEGv, NDEBv
STOP
ENDIF
I640= 0
I644= 0
c** Subroutine MASSES returns the names of the atoms NAMEi,ground
c electronic state degeneracy GELi, nuclear spin degeneracy GNSi,
c mass MASSi, and isotopic abundance ABUNDi for a given atomic isotope.
IF((IAN1.GT.0).AND.(IAN1.LE.109)) THEN
CALL MASSES(IAN1,IMN1,NAME1,GEL1,GNS1,MASS1,ABUND1)
ELSE
c** If particle-i is not a normal atomic isotope, read a 2-character
c NAME (enclosed between ', as in 'mu') and its actual mass.
c----------------------------------------------------------------------
READ(5,*) NAME1, MASS1
c----------------------------------------------------------------------
ENDIF
IF((IAN2.GT.0).AND.(IAN2.LE.109)) THEN
CALL MASSES(IAN2,IMN2,NAME2,GEL2,GNS2,MASS2,ABUND2)
ELSE
c----------------------------------------------------------------------
READ(5,*) NAME2, MASS2
c----------------------------------------------------------------------
ENDIF
ZMU= (MASS1*MASS2)/(MASS1+MASS2 - DFLOAT(CHARGE)*ZMASE)
c** Numerical factor 16.857629206 (+/- 0.000,000,013) based on Compton
c wavelength of proton & proton mass (u) from 2011 physical constants.
CU= 16.857629206d0/ZMU
SQCU= DSQRT(CU)
c=======================================================================
c TITLE is a title or output header of up to 78 characters, read on a
c single line enclosed between single quotes: e.g. 'title of problem'
c=======================================================================
READ(5,*) TITLE
c-----------------------------------------------------------------------
WRITE(6,600) TITLE,NAME1,IMN1,NAME2,IMN2,CHARGE,ZMU,CU,MASS1,
1 MASS2
WRITE(6,602) TOLER,NBIS,NGP
c=======================================================================
IF((NDEGv.LE.0).OR.(NDEGv.GE.2)) THEN
c** For Dunham or MXS vibrational energy function, read in order LMAXGv
c and values YL0(i) (i=1,LMAXGv) of Dunham vibrational coefficients
c-----------------------------------------------------------------------
READ(5,*) LMAXGv
READ(5,*) (YL0(L),L= 1,LMAXGv)
c-----------------------------------------------------------------------
WE= YL0(1)
WEXE= -YL0(2)
c=======================================================================
c** If using Tellinghuisen-style Mixed Representations ... read
c VS the v value (floating point) where Dunham switches to NDE,
c DVS the width parameter on the switching function, and
c DLIM the well depth, or energy at the asymptote assuming an energy
c zero at v= -1/2
c=======================================================================
IF(NDEGv.GE.2) READ(5,*) VS, DVS, DLIM
c-----------------------------------------------------------------------
IF(NDEGv.GE.2) WRITE(6,604) CCDC(0),LMAXGv,VS,DVS,DLIM
WRITE(6,606) LMAXGv,CCDC(0),(YL0(I),I= 1,LMAXGv)
ENDIF
IF(NDEGv.GE.1) THEN
c
c** If use NDE or MXS expression for vibrational energies, read NDE
c control parameters and expansion constants here.
c* For an "outer" Pade expansion, ITYPE=1 ; for an "inner" Pade,
c ITYPE=2 ; if ITYPE=3, use an exponential polynomial NDE .
c* Expansion variable is (vD-v), and the leading non-zero contribution
c to the NP-term numerator polynomial is the power IZP0 of (vD-v),
c while the corresponding leading term in the NQ0-term denominator
c polynomial is (vD-v)**IZQ0
c-----------------------------------------------------------------------
READ(5,*) NLR, ITYPE, IZP0, IZQ0, NP0, NQ0, VD, XCN0
IF(NP0.GT.0) READ(5,*) (P0(I),I= 1,NP0)
IF(NQ0.GT.0) READ(5,*) (Q0(I),I= 1,NQ0)
c-----------------------------------------------------------------------
IF(NDEGv.EQ.1) DLIM= 0.d0
IZP0= IZP0- 1
IZQ0= IZQ0- 1
KORDR= 0
FF= V00 - 1.d-5
DO I= 1,3
CALL NDEDKM(FF,KORDR,GV(I),DGDV(I),NLR,XCN0,DLIM,VD,
1 IZP0,IZQ0,ITYPE,NP0,NQ0,P0,Q0)
FF= FF+ 1.d-5
ENDDO
IF(NDEGv.EQ.1) DLIM= -GV(2)
WRITE(6,608) CCDC(0),NP0,NQ0,(PTYPE(I),I= ITYPE,6,3),0,NLR,
1 XCN0,IZP0+1,IZQ0+1,VD,DLIM
IF(NP0.GT.0) WRITE(6,610) (P0(I),I= 1,NP0)
IF(NQ0.GT.0) WRITE(6,612) (Q0(I),I= 1,NQ0)
WE= DGDV(2)
WEXE= (DGDV(1)-DGDV(3))/4.D-5
IF(NDEGv.GE.2) THEN
SwLR= DEXP(-(VS+ 0.5d0)/DVS)
Sw= 1.d0/(1.d0+ SwLR)
SwLR= SwLR/(1.d0 + SwLR)
WE= SwLR*WE + Sw*YL0(1)
WEXE= SwLR*WEXE - Sw*YL0(2)
ENDIF
ENDIF
c
IF(NDEBv.LT.0) THEN
c** If have NO v-dependent rotational data and wish to define inner
c potential wall by Morse function, read Req value to define position
c of potential minimum.
c-----------------------------------------------------------------------
READ(5,*) Req
c-----------------------------------------------------------------------
c** Use vibrational constants to determine Morse parameters at minimum.
BE= CU/Req**2
AE= (SQRT(WEXE/BE) - 1.d0)*6.d0*BE**2/WE
De= 0.25d0*WE**2/WEXE
BETA= DSQRT(DABS(WEXE)/CU)
WRITE(6,614) Req,WE,WEXE,De,BETA
MORSE= 1
ENDIF
IF((NDEBv.EQ.0).OR.(NDEBv.GE.2)) THEN
c** For Dunham or MXS Bv function, read in order LMAXBv of Dunham
c expansion in (v+1/2) and coefficients YL1(i) (i=0,LMAXBv)
c-----------------------------------------------------------------------
READ(5,*) LMAXBv
IF(LMAXBv.GE.0) READ(5,*) (YL1(L),L= 0,LMAXBv)
c-----------------------------------------------------------------------
c** Use vibrational constants to determine Morse parameters at minimum.
BE= YL1(0)
AE= -YL1(1)
IF(NDEBv.GE.2) WRITE(6,604) CCDC(1),LMAXBv,VS
WRITE(6,606) LMAXBv+1,CCDC(1),(YL1(I),I= 0,LMAXBv)
ENDIF
IF(NDEBv.GE.1) THEN
c** If use NDE or MXS expression for Bv values, read NDE control
c parameters and expansion constants which are defined in exactly the
c same way as those for the vibrational NDE.
c-----------------------------------------------------------------------
READ(5,*) ITYPB, IZP1, IZQ1, NP1, NQ1, XCN1
IF(NP1.GT.0) READ(5,*) (P1(I),I= 1,NP1)
IF(NQ1.GT.0) READ(5,*) (Q1(I),I= 1,NQ1)
c-----------------------------------------------------------------------
WRITE(6,608) CCDC(1),NP1,NQ1,(PTYPE(I),I= ITYPB,6,3),1,NLR,
1 XCN1,IZP1,IZQ1
IF(NP1.GT.0) WRITE(6,610) (P1(I),I= 1,NP1)
IF(NQ1.GT.0) WRITE(6,612) (Q1(I),I= 1,NQ1)
IZP1= IZP1- 1
IZQ1= IZQ1- 1
KORDR= 1
FF= V00
CALL NDEDKM(FF,KORDR,BE,AE,NLR,XCN1,DLIM,VD,
1 IZP1,IZQ1,ITYPB,NP1,NQ1,P1,Q1)
AE= -AE
IF(NDEBv.GE.2) THEN
IF(MAX(DABS(BE/YL1(0)),DABS(AE/YL1(1))).GT.0.01d0/SwLR)
1 WRITE(6,609) FF,SwLR,YL1(0),BE,-YL1(1),AE
BE= SwLR*BE + Sw*YL1(0)
AE= SwLR*AE - Sw*YL1(1)
ENDIF
ENDIF
c** Kaiser (integer) controls sophistication of RKR calculation.
c If(Kaiser.le.0) do simple first order with Y00=0
c If(Kaiser > 0) apply "Kaiser" correction by setting lower bound of
c integration as v(min)= -1/2 - Y00/Y10 where Dunham Y00 & Y10
c calculated from energy derivatives evaluated at v=-1/2.
c** NSV is the number of different v-increments to be used in defining
c the set of v-values for which turning points are to be calculated.
c** If(VEXT.le.0) perform RKR calculation with no inner wall smoothing
c but calculate & print exponent coefficients cEXT Cext of exponential
c functions approximately fitted to the 3 preceeding inner turning pts.
c** If(VEXT.gt.0) extrapolate inner wall above v=VEXT with an
c exponential exactly fitted to last 3 inner T.P. with v.le.VEXT
c-----------------------------------------------------------------------
READ(5,*) Kaiser, NSV, VEXT
c-----------------------------------------------------------------------
IF(MORSE.EQ.1) THEN
VEXT= -1.d0
Kaiser= 0
ENDIF
IF(Kaiser.GT.0) then
c** For (Kaiser > 0) calculate Y00 and vib quantum no. at minimum V00 and
c in NDE case, adjust D value to actual well depth De= G(v=-1/2)+Y00
Y00= AE*WE/(BE*12.d0)
Y00= (BE- WEXE)/4.d0 + Y00 + Y00**2/BE
IF(NDEGv.GE.1) DLIM= DLIM + Y00
DV0= -Y00/WE
V00= V00+ DV0
c* Iterate to correct for linear Y00/WE approximation
I= 1
VX(1)= V00
CALL GVBV(I,VX,GV,DGDV,BV)
DV0= -GV(1)/WE
V00= V00+ DV0
DV0= DV0- Y00/WE
WRITE(6,618) Y00,DV0,V00,WE,WEXE,BE,AE
IF(NDEGv.GE.1) WRITE(6,620) DLIM
ENDIF
I=1
VX(I)= V00
CALL GVBV(I,VX,GV,DGDV,BV)
IF(MORSE.NE.1) BE= BV(1)
REQ= DSQRT(CU/BE)
WRITE(6,622) V00,GV(1),DGDV(1),-WEXE,BE,REQ,AE
IF(VEXT.GT.0.d0) WRITE(6,624) VEXT
DO I= 1, NSV
c** Read and generate v-s at which turning points desired
c** For each of the NSV cases, generate v values ranging from V1(j)
c to V2(j) with (positive) increment DV(j) (for j= 1 to NSV).
c* The resulting values should increase monotonically.
c-----------------------------------------------------------------------
READ(5,*,end=99) V1(I), DV(I), V2(I)
c-----------------------------------------------------------------------
ENDDO
IF((NDEGv.GE.1).AND.(V2(NSV).GT.VD)) V2(NSV)= DINT(VD)
I= 0
DO J= 1,NSV
NN= (V2(J)- V1(J)+ 1.d-4)/DV(J)
IF(J.LT.NSV) NN= (V1(J+1)-V1(J))/DV(J)- 1
I= I+1
V(I)= V1(J)
DO K= 1,NN
I= I+1
V(I)= V(I-1) + DV(J)
ENDDO
IF(DABS(V(I)).LT.1.D-12) V(I)= 0.d0
ENDDO
NVIB= I
NTP= 2*NVIB+1
WRITE(6,626) NVIB,(V(I),I= 1,NVIB)
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c++ Commence actual turning point calculation here +++++++++++++++++++++
WRITE(6,630)
I638= 0
DEI= 0.d0
DRI= 0.d0
c** Commence loop over the NVIB values of vib quantum no. v ********
DO IVB= 1,NVIB
I634= 0
VUP= V(IVB)
VX(1)= VUP
I= 1
CALL GVBV(I,VX,GV,DGDV,BV)
GUP= GV(1)
BUP= BV(1)
DGUP= DGDV(1)
IF(DGUP.LE.0.d0) THEN
WRITE(6,632) VUP,GUP,DGUP
IVIB= IVB-1
NTP= 2*IVIB+ 1
I1= NVIB- IVIB+ 1
GO TO 72
ENDIF
F= 0.d0
G= 0.d0
NDIV= 1
TSTF= 9.d99
TSTG= 9.d99
c** In attempting to achieve integral convergence, consider up to NBIS
c bisections of the interval.
DO 50 IBIS= 0,NBIS
TSTFB= TSTF
TSTGB= TSTG
VFN= V00
FF= 0.d0
GG= 0.d0
RANGE= (VUP-V00)/DFLOAT(NDIV)
c** Sum over contributions from the NDIV segments of the interval
DO IDIV= 1,NDIV
VST= VFN
VFN= VST+RANGE
IF(IDIV.NE.NDIV) THEN
c** Quadrature weights and points for range with no singularity
D1= 0.5d0*(VFN+ VST)
D2= 0.5d0*(VFN- VST)
DO I= 1,NGP
WW(I)= WG(I)
VX(I)= D1+D2*XG(I)
ENDDO
ELSE
VFN= VUP
FF= FF/2.d0
GG= GG/2.d0
c** Weights and points for range with 1/sqrt(E-V) singularity
D2= RANGE
DO I= 1,NGP
WW(I)= W2(I)
VX(I)= VST+X2(I)*D2
ENDDO
ENDIF
CALL GVBV(NGP,VX,GV,DGDV,BV)
c** Actual quadrature loop starts here
DO I= 1,NGP
TT= WW(I)/DSQRT(GUP-GV(I))
FF= FF+TT
GG= GG+TT*BV(I)
ENDDO
ENDDO
c** End of quadrature - now test for convergence
FB= F
GB= G
F= FF*RANGE*SQCU
G= GG*RANGE/SQCU
TSTF= DABS(1.d0- FB/F)
TSTG= 0.d0
IF(((NDEGv.NE.0).OR.(LMAXBv.GE.1)).AND.(VUP.le.VEXT))
1 TSTG= DABS(1.d0- GB/G)
IF(MORSE.EQ.1) TSTG= 0.d0
IF(DMAX1(TSTF,TSTG).LT.TOLER) GO TO 54
IF((TSTF.GT.TSTFB).OR.(TSTG.GT.TSTGB)) THEN
I634= I634+ 1
IF((I634.GE.1).AND.(NDIV.GT.4)) THEN
WRITE(6,634) NDIV, TSTF,TSTFB, TSTG,TSTGB
GOTO 54
ENDIF
ENDIF
NDIV= 2*NDIV
50 CONTINUE
NDIV= NDIV/2
WRITE(6,636) NDIV,TSTF,TSTG,TOLER
54 IF(MORSE.EQ.1) THEN
c** Use Morse inner turning point when no rotational data input
RMIN= Req - DLOG(1.d0 + DSQRT(GUP/De))/BETA
RMAX= RMIN + F + F
BUP= BE
ELSE
HEL= DSQRT(F/G+F*F)
RMIN= HEL-F
RMAX= HEL+F
ENDIF
I1= NVIB+1-IVB
I2= NVIB+1+IVB
IF(IVB.LE.1) GO TO 60
DEIB= DEI
DRIB= DRI
I1P1= I1+1
RMIN2= RT(I1P1)
DEI= GUP- ET(I1P1)
DRI= RMIN- RMIN2
IF(IVB.LE.2) GO TO 60
IF((VEXT.GT.0.d0).AND.(VUP.GT.VEXT)) GO TO 56
IF(RMIN.GE.RT(I1P1)) THEN
c** If inner wall turns over, print warning
IF(I638.LE.0) THEN
I638= I1P1
WRITE(6,638) VUP
ENDIF
GOTO 60
ENDIF
IF((MORSE.GT.0).OR.(I638.GT.0)) GO TO 60
c** Use differences to get rough estimate of exponent coefficient CEXT
c for local exponential fit to inner wall.
I1P2= I1+2
RMIN3= RT(I1P2)
VRAT= DEI/DEIB
CEXT= 2.d0*(DRI-VRAT*DRIB)/(DRI**2 + VRAT*DRIB**2)
c** Iteratively converge on exact (to machine precision) value of CEXT
IF((DABS(CEXT*RMIN3).GT.70.D0).AND.(VUP.GT.0.d0)) THEN
IF(I640.LE.0) WRITE(6,640)
I640= 1
GO TO 55
ENDIF
VRAT= (GUP-ET(I1P1))/(GUP-ET(I1P2))
ADCEXT= 1.d99
DO I= 1,15
BDCEXT= ADCEXT
EXP1= 1.d0
EXP2= DEXP(-CEXT*(RT(I1P1)-RMIN))
EXP3= DEXP(-CEXT*(RT(I1P2)-RMIN))
FUN= (EXP1-EXP2)/(EXP1-EXP3)-VRAT
DFUN= ((RMIN*EXP1-RT(I1P1)*EXP2) - (RMIN*EXP1-
1 RT(I1P2)*EXP3)*(EXP1-EXP2)/(EXP1-EXP3))/(EXP1-EXP3)
DCEXT= FUN/DFUN
ADCEXT= DABS(DCEXT)
IF((ADCEXT.GE.BDCEXT).AND.(ADCEXT.LT.1.d-10)) GO TO 230
IF(ADCEXT.LE.0.d0) GO TO 230
CEXT= CEXT+DCEXT
ENDDO
WRITE(6,642) DCEXT
230 CONTINUE
IF(CEXT.LE.0.d0) THEN
I644= I644+ 1
IF(I644.EQ.1) WRITE(6,644) VUP
ENDIF
EXP1= DEXP(-CEXT*RMIN)
EXP2= EXP2*EXP1
BEXT= (GUP-ET(I1P1))/(EXP1-EXP2)
AEXT= GUP-BEXT*EXP1
c** Parameters for possible inner wall extrapolation now determined.
c
55 WRITE(6,646) VUP,GUP,DGUP,BUP,RMIN,RMAX,NDIV,TSTF,TSTG,CEXT
GO TO 64
c** Apply option to smoothly extrapolate inner wall above v=VEXT
c using a simple exponential
56 REXT= DLOG(BEXT/(GUP-AEXT))/CEXT
DR1= REXT-RMIN
RMIN= REXT
RMAX= RMAX+DR1
WRITE(6,646)VUP,GUP,DGUP,BUP,RMIN,RMAX,NDIV,TSTF,TSTG,CEXT,DR1
GO TO 64
60 WRITE(6,646) VUP,GUP,DGUP,BUP,RMIN,RMAX,NDIV,TSTF,TSTG
64 RT(I1)= RMIN
RT(I2)= RMAX
ET(I1)= GUP
ET(I2)= GUP
ENDDO
I1= 1
I2= NTP
72 IF(VEXT.GT.0.d0) WRITE(6,648) VEXT,AEXT,BEXT,CEXT
c** Write ordered turning points compactly on channel 7
c** To facilitate use of resulting potential array, add 5 extrapolated
c inner-wall points in the output file
IF((CEXT.GT.0.d0).AND.(CEXT.LT.20.d0)) THEN
IDR= -INT(1.0d3*ET(I1)*(RT(I1+1)-RT(I1))/(ET(I1+1)-ET(I1)))
IR1= INT(1.d4*RT(I1))
DO I= 1,5
RT(I1-I)= 1.D-4*(IR1- I*IDR)
ET(I1-I)= AEXT + BEXT*DEXP(-CEXT*RT(I1-I))
ENDDO
I1= I1 - 5
I2= NTP + I1 + 4
ENDIF
RT(NVIB+1)= REQ
ET(NVIB+1)= 0.d0
IF(I1.LT.0) WRITE(7,700) TITLE,NTP+5,ZMU
IF(I1.GT.0) WRITE(7,700) TITLE,NTP,ZMU
WRITE(7,702) (RT(I),ET(I),I= I1,I2)
WRITE(6,650)
GO TO 2
99 STOP
598 FORMAT(/' *** Input ERROR *** NDEBv=',i3,' .GT. NDEGv=',i3,
1 ' is NOT allowed!')
600 FORMAT(1x,A78/1x,35('**')/' RKR potential for ',A2,'(',I3,')-',
1 A2,'(',I3,') with Charge=',I2/' Reduced mass ZMU=',
2 F15.11,' and constant C_u/ZMU =',F16.12/5x,'from atomic masse
3s:',f15.10,' & ',F15.10,'(u)')
602 FORMAT(/' Seek relative quadrature convergence',1PD8.1,
1 '. Bisect interval up to',i3,' times.'/5x,'performing',i3,
2 '-point Gaussian quadrature in each segment')
604 FORMAT(/' Represent ',A2,"'s by Tellinghuisen-type MXS mixed rep
1resentation:"/ 1x,16('==')/I5,"'th order Dunham for v .le. VS &
2 NDE for v > VS, with VS=",F8.4:/4x,'with switching functio
3n F_s = 1/[1 + exp{(v-VS)/DVS}] with DVS=',f7.4/5x,'and a
4sympotote energy (dissociation limit) DLIM=',F11.4,' [cm-1]')
606 FORMAT(/' The',i3,' Dunham ',A2,' expansion coefficients are'/
1 (6X,1P4D18.10:))
608 FORMAT(/' NDE for ',A2,' is an (NP=',I2,'/NQ=',I2,') ',2A6,
1 ' expansion in (vD-v) with'/8x,'X',I1,'(n=',I1,')=',1PD14.7:
2 ' and leading num. and denom. powers',I3,' & ',I3:/8X,'vD=',
3 0PF12.6:' D-G(v=-1/2)=',F14.6)
609 FORMAT(/' *** CAUTION *** DVS may be too large. At v=',f7.3,
1 ' where Sw(LR)=',1PD9.2/6x,'Be(Dun)=',d8.1,' Be(NDE)=',
2 d8.1,' ae(Dun)=',d8.1,' ae(NDE)=',d8.1)
610 FORMAT(5X,'Numerator coefficients are:',2(1PD20.12:)/
1 (12X,3D20.12:))
612 FORMAT(5X,'Denominator coefficients : ',2(1PD20.12:)/
1 (12X,3D20.12:))
614 FORMAT(/' NO rotational constants input, so inner wall of potentia
1l is Morse function.'/' Input Req=',f9.6,'(Angst) plus we=
2',f9.3,' & wexe=',f9.5,' [cm-1]'/' yields Morse with De=',
3 f10.3,' [cm-1] and beta=',f9.6,' [1/Angst.]')
618 FORMAT(/' Calculate Y00=',F13.9,4X,'v(cor)=',F14.10,4x,'v(min)='
1 ,F14.10/5x,'using we=',f9.4,' wexe=',f9.6,' Be=',f10.6,
2 ' ae=',f11.8)
620 FORMAT(5x,'and corrected effective De=',F13.6,' (after adding Y0
10)')
622 FORMAT(/' At v00=',F9.5,' Gv=',F12.8,' dG/dv=',F9.4,
1 ' (1/2)d2G/dv2=',F10.6/21x,'Bv=',F12.8,' { ==> Req=',F12.9,
2 '(A) }'/21x,'alpha_e =',F12.9)
624 FORMAT(/' Above v =',f8.3,' extrapolate inner wall with expone
1ntial'/9x,'fitted to last 3 points ( & shift RMAX accordingly)')
626 FORMAT(/' Calculate turning points at the',i4,' v-values'/
1 (1X,11F7.2:))
630 FORMAT(/' Resulting Turning Points:'/
1 5X,'v',8X,'E(v)',4X,'dE(v)/dv',7X,'B(v)',10X,'Rmin(v)',8X,
2 'Rmax(v)',3X,'NDIV tst(f) tst(g)',4X,'C(exp)',7X,'d(RMIN)'/
3 1X,61('**'))
632 format(/' STOP calculation at v=',f6.2,' where E(v)=',f8.2,
1 ' & dE(v)/dv =',G14.7)
634 FORMAT(' *** STOP ITERATION: At NDIV=',i3,
1 ' tst(f)/(previous)=',1Pd8.1,'/',d7.1,' tst(g)/(previous)=',
2 d8.1,'/',d7.1)
636 FORMAT(' *** CAUTION:',i3,' interval incomplete convergence: tst(
1f) & tst(g)=',1P2d8.1,' while TOLER=',d8.1)
638 FORMAT(' *** WARNING *** inner wall becomes unstable at v =',
1 F6.2,' where RMIN turns outward!')
640 FORMAT(' *** CAUTION *** inner wall exponent parameter becomes ver
1y large so skip converging it.')
642 FORMAT(' *** CAUTION *** FAIL to converge C(exp) after 15 tries.',
1 ' Last step =',1PD10.2)
644 FORMAT(' *** CAUTION *** Inner potential wall has negative curvat
1ure and requires smoothing for VEXT .ge.',f7.2)
646 FORMAT(F7.3,F12.4,F11.4,F15.10,2F15.10,I4,1P2D9.1,0PF11.6,F15.10)
648 FORMAT(1X,61('**')//' For v .GE.',f6.2,' inner wall extrapolate
1d as: V(R) =',F13.4,' +',D15.8,'*exp(-',f12.8,'*R)')
650 FORMAT(1x,61('**')/////)
700 FORMAT(/1x,A78/' NTP=',I5,' RKR turning points for mu=',f14.10)
702 FORMAT((F20.14,F19.11))
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c***********************************************************************
SUBROUTINE GVBV(NL,VX,GV,DG,BV)
c** At the NL values of the vibrational quantum number VX, generate
c values of Bv, Gv and their first derivatives w.r.t. v from
c (i) pure Dunham (NDEXv= 0), (ii) pure NDE (NDEXv= 1) or
c (iii) Tellinghuisen 'mixed' MXS function NDEXv .ge. 2 where
c require NDEBv .le. NDEGv
c
INTEGER MXDUN !! max polynomial order for vib expansions
PARAMETER (MXDUN=25)
c*** define types for local variables
INTEGER I,L,KORDR,NL
REAL*8 VX(NL),GV(NL),DG(NL),BV(NL),KM,DKM,Sw,SwLR,VPH,YY,DY
c** Common block for Dunham & MXS function parameters
INTEGER LMAXGv,LMAXBv,NDEGv,NDEBv
REAL*8 VS,DVS,Y00,YL0(0:MXDUN),YL1(0:MXDUN)
COMMON /DUNPRM/Y00,VS,DVS,YL0,YL1,LMAXGv,LMAXBv,NDEGv,NDEBv
c** Common block for NDE function parameters
INTEGER NLR,ITYPE,ITYPB,IZP0,IZQ0,IZP1,IZQ1,NP0,NQ0,NP1,NQ1
REAL*8 VD,DLIM,XCN0,XCN1,P0(MXDUN),Q0(MXDUN),P1(MXDUN),Q1(MXDUN)
COMMON /NDEPRM/VD,DLIM,XCN0,XCN1,P0,Q0,P1,Q1,NLR,ITYPE,ITYPB,
1 IZP0,IZQ0,IZP1,IZQ1,NP0,NQ0,NP1,NQ1
c-----------------------------------------------------------------------
DO I= 1, NL
c** Loop over all levels, calculating energies, Bv's and their first
c derivatives w.r.t. v
VPH= VX(I)+ 0.5d0
c** First - for Dunham vib energy (or Dunham part of MXS energy)
IF(NDEGv.NE.1) THEN
YY= YL0(LMAXGv)
DY= LMAXGv*YY
DO L= LMAXGv-1, 1, -1
YY= YY*VPH + YL0(L)
dY= DY*VPH + L*YL0(L)
ENDDO
YY= YY*VPH + Y00
GV(I)= YY
DG(I)= DY
ENDIF
IF(NDEGv.GE.1) THEN
c** For NDE or MXS vibrational energy .....
KORDR= 0
CALL NDEDKM(VX(I),KORDR,KM,DKM,NLR,XCN0,DLIM,VD,
1 IZP0,IZQ0,ITYPE,NP0,NQ0,P0,Q0)
IF(NDEGv.EQ.1) THEN
GV(I)= KM
DG(I)= DKM
ELSE
SwLR= DEXP((VX(i)- VS)/DVS)
Sw= 1.d0/(1.d0+ SwLR)
SwLR= SwLR * Sw
GV(I)= Sw*YY + SwLR*KM
DG(I)= Sw*DY + SwLR*DKM - (YY - KM)*Sw*SwLR/DVS
ENDIF
ENDIF
IF(NDEBv.NE.1) THEN
c
c** First - for Dunham Bv value (or Dunham part of MXS function for Bv)
DY= 1.d0
YY= 0.d0
IF(LMAXBv.GE.0) THEN
YY= YL1(LMAXBv)
DY= LMAXBv*YY
DO L= LMAXBv-1, 0, -1
YY= YY*VPH + YL1(L)
DY= DY*VPH + L*YL1(L)
ENDDO
BV(I)= YY
ENDIF
ENDIF
IF(NDEBv.GE.1) THEN
c** For NDE or MXS function for Bv value ...
KORDR= 1
CALL NDEDKm(VX(I),KORDR,KM,DKM,NLR,XCN1,DLIM,VD,
1 IZP1,IZQ1,ITYPB,NP1,NQ1,P1,Q1)
IF(NDEBv.EQ.1) THEN
BV(I)= KM
ELSE
BV(I)= Sw*YY + SwLR*KM
ENDIF
ENDIF
ENDDO
RETURN
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c**********************************************************************
SUBROUTINE NDEDKM(VV,KORDR,KM,DKM,NLR,XCN,DLIM,VD,IZP,IZQ,ITYP,
1 NP,NQ,P,Q)
c** Subroutine which uses Near-Dissociation Expansions to generate
c predicted value KM and vibrational first derivative DKM of band
c constant of rotational order KORDR (KORDR= 0 for Gv, =1 for Bv,
c =2 for -Dv, =3 for Hv, ... etc., for vibrational level (integer) IV
c-----------------------------------------------------------------------
INTEGER KORDR,IZP,IZQ,ITYP,IPW,J,NLR,NP,NQ
REAL*8 KM,DKM,XCN,DLIM,VD,P(NP),Q(NQ), PW,PAB,DV,DVP,SMN,SMD,DSMN,
1 DSMD,VV
IPW= 2*NLR/(NLR-2)
PW= 2.d0*NLR/(DFLOAT(NLR)- 2.d0)
IF(KORDR.GE.1) THEN
IPW= IPW - 2*KORDR
PW= PW - 2.d0*KORDR
ENDIF
PAB= 1.d0
IF(ITYP.EQ.2) PAB= PW
DV= VD- VV
c** Evaluate the NDE numerator & denominator polynomials
SMN= 1.d0
SMD= 1.d0
DSMN= 0.d0
DSMD= 0.d0
IF(NP.GT.0) THEN
c ... numerator polynomial ...
DVP= 1.d0
IF(IZP.GT.0) DVP= DV**IZP
DO J= 1,NP
DSMN= DSMN+(J+IZP)*P(J)*DVP
DVP= DVP*DV
SMN= SMN+P(J)*DVP
ENDDO
ENDIF
IF(NQ.GT.0) THEN
c ... denominator polynomial ...
DVP= 1.d0
IF(IZQ.GT.0) DVP= DV**IZQ
DO J= 1,NQ
DSMD= DSMD+(J+IZQ)*Q(J)*DVP
DVP= DVP*DV
SMD= SMD+Q(J)*DVP
ENDDO
ENDIF
IF(ITYP.EQ.3) THEN
KM= XCN*DV**PW *DEXP(SMN- 1.d0)
DKM= -KM*(PW/DV + DSMN)
IF(KORDR.EQ.0) THEN
KM= DLIM - KM
DKM= -DKM
ENDIF
ELSE
IF(ITYP.EQ.1) THEN
KM= XCN*SMN/SMD
DKM= DV
ENDIF
IF(ITYP.EQ.2) THEN
KM= XCN
DKM= DV*SMN/SMD
ENDIF
IF(NLR.EQ.5) KM= KM*DABS(DKM)**PW
IF(NLR.NE.5) KM= KM*DKM**IPW
DKM= - KM*(PW/DV + PAB*(DSMN/SMN- DSMD/SMD))
IF(KORDR.EQ.0) THEN
KM= DLIM - KM
DKM= -DKM
ENDIF
ENDIF
RETURN
END
c=======================================================================
c***********************************************************************
SUBROUTINE MASSES(IAN,IMN,NAME,GELGS,DGNS,MASS,ABUND)
c***********************************************************************
c** For isotope with (input) atomic number IAN and mass number IMN,
c return (output): (i) as the right-adjusted 2-character variable NAME
c the alphabetic symbol for that element, (ii) the ground state
c electronic degeneracy GELGS, (iii) the nuclear spin degeneracy DGNS,
c (iv) the atomic mass MASS [amu], and (v) the natural isotopic
c abundance ABUND [in percent]. GELGS values based on atomic states
c in Moore's "Atomic Energy Level" tables, the isotope masses are taken
c from the 2012 mass table [Wang, Audi, Wapstra, Kondev, MacCormick, Xu
c & Pfeiffer, Chin.Phys.C 36, 1603-2014 (2012)] ,the proton, deuteron,
c and triton masses are taken from the 2010 fundamental constants table
c [Mohr, Taylor, & Newell, Rev. Mod. Phys. 84, 1587-1591 (2012)] and other
c quantities from Tables 6.2 and 6.3 of "Quantities, Units and Symbols in
c Physical Chemistry", by Mills et al.(Blackwell,2'nd Edition, Oxford,1993).
c** If the input value of IMN does not equal one of the tabulated values
c for atomic species IAN, return the abundance-averaged standard atomic
c weight of that atom and set DGNS=-1 and ABUND=-1.
c** For Atomic number IAN=0 and isotope mass numbers IMN=1-3, return the
c masses of the proton, deuteron, and triton, p,d & t, respectively
c Masses and properties of selected Halo nuclei an unstable nuclei included
c COPYRIGHT 2005-2015 : last updated 10 January 2016
c** By R.J. Le Roy, with assistance from
c G.T. Kraemer, J.Y. Seto and K.V. Slaughter.
c***********************************************************************
REAL*8 zm(0:123,0:15),mass,ab(0:123,15),abund
INTEGER i,ian,imn,gel(0:123),nmn(0:123),mn(0:123,15),
1 gns(0:123,15),DGNS,gelgs
CHARACTER*2 NAME,AT(0:123)
cc
DATA at(0),gel(0),nmn(0),(mn(0,i),i=1,3)/' p',1,3,1,2,3/
DATA (zm(0,i),i=0,3)/1.008d0,1.007276466812d0,2.013553212712d0,
2 3.0155007134d0/
DATA (gns(0,i),i=1,3)/2,3,2/
DATA (ab(0,i),i=1,3)/0.d0, 0.d0, 0.d0/
c
DATA at(1),gel(1),nmn(1),(mn(1,i),i=1,3)/' H',2,3,1,2,3/
DATA (zm(1,i),i=0,3)/1.00794d0, 1.00782503223d0, 2.01410177812d0,
1 3.0160492779d0/
DATA (gns(1,i),i=1,3)/2,3,2/
DATA (ab(1,i),i=1,3)/99.985d0,0.015d0,0.d0/
c
DATA at(2),gel(2),nmn(2),(mn(2,i),i=1,4)/'He',1,4,3,4,6,8/
DATA (zm(2,i),i=0,4)/4.002602d0, 3.0160293201d0, 4.00260325413d0,
1 6.0188891d0, 8.033922d0/
DATA (gns(2,i),i=1,4)/2,1,1,1/
DATA (ab(2,i),i=1,4)/0.000137d0,99.999863d0, 2*0.d0/
c
DATA at(3),gel(3),nmn(3),(mn(3,i),i=1,6)/'Li',2,6,6,7,8,9,11,12/
DATA (zm(3,i),i=0,6)/6.941d0, 6.0151228874d0, 7.016003437d0,
1 8.02248736d0,9.0267895d0,11.043798d0,12.05378d0/
DATA (gns(3,i),i=1,6)/3,4,5,4,4,1/
DATA (ab(3,i),i=1,6)/7.5d0, 92.5d0, 4*0.d0/
c
DATA at(4),gel(4),nmn(4),(mn(4,i),i=1,8)/'Be',1,8,7,9,10,11,12,
1 14,15,16/
DATA (zm(4,i),i=0,8)/9.012182d0, 7.01692983d0, 9.01218307d0,
1 10.0135338d0, 11.021658d0, 12.026921d0, 14.04289d0, 15.05346d0,
2 16.06192d0/
DATA (gns(4,i),i=1,8)/4,4,3,2,1,1,2,1/
DATA (ab(4,i),i=1,8)/0.d0, 100.d0, 6*0.d0/
c
DATA at(5),gel(5),nmn(5),(mn(5,i),i=1,10)/' B',2,10,8,10,11,12,
1 13,14,15,17,18,19/
DATA (zm(5,i),i=0,10)/10.811d0, 8.0246072d0, 10.0129369d0,
1 11.0093054d0, 12.0143521d0, 13.0177802d0, 14.025404d0,
2 15.031103d0, 17.04699d0, 18.05617d0,19.06373d0/
DATA (gns(5,i),i=1,10)/5,7,4,3,4,5,4,4,1,4/
DATA (ab(5,i),i=1,10)/0.d0, 19.9d0,80.1d0, 7*0.d0/
c
DATA at(6),gel(6),nmn(6),(mn(6,i),i=1,14)/' C',1,14,9,10,11,12,13,
1 14,15,16,17,18,19,20,21,22/
DATA (zm(6,i),i=0,14)/12.011d0, 9.0310367d0, 10.0168532d0,
1 11.0114336d0, 12.d0, 13.00335483507d0, 14.003241989d0,
1 15.0105993d0, 16.014701d0, 17.022586d0, 18.02676d0, 19.03481d0,
2 20.04032d0, 21.04934d0, 22.05720d0/
DATA (gns(6,i),i=1,14)/4,1,4,1,2,1,2,1,4,1,2,1,2,1/
DATA (ab(6,i),i=1,14)/3*0.d0, 98.90d0,1.10d0, 9*0.d0/
c
DATA at(7),gel(7),nmn(7),(mn(7,i),i=1,2)/' N',4,2,14,15/
DATA (zm(7,i),i=0,2)/14.00674d0, 14.00307400443d0,15.0001088989d0/
DATA (gns(7,i),i=1,2)/3,2/
DATA (ab(7,i),i=1,2)/99.634d0,0.366d0/
c
DATA at(8),gel(8),nmn(8),(mn(8,i),i=1,3)/' O',5,3,16,17,18/
DATA (zm(8,i),i=0,3)/15.9994d0, 15.99491461957d0, 16.9991317565d0,
1 17.9991596129d0/
DATA (gns(8,i),i=1,3)/1,6,1/
DATA (ab(8,i),i=1,3)/99.762d0, 0.038d0, 0.200d0/
c
DATA at(9),gel(9),nmn(9),(mn(9,i),i=1,1)/' F',4,1,19/
DATA (zm(9,i),i=0,1)/18.9984032d0, 18.9984031627d0/
DATA (gns(9,i),i=1,1)/2/
DATA (ab(9,i),i=1,1)/100.d0/
c
DATA at(10),gel(10),nmn(10),(mn(10,i),i=1,4)/'Ne',1,4,17,20,21,22/
DATA (zm(10,i),i=0,4)/20.1797d0, 17.017672d0, 19.9924401762d0,
1 20.99384669d0,21.991385115d0/
DATA (gns(10,i),i=1,4)/2,1,4,1/
DATA (ab(10,i),i=1,4)/0.d0, 90.48d0, 0.27d0, 9.25d0/
c
DATA at(11),gel(11),nmn(11),(mn(11,i),i=1,1)/'Na',2,1,23/
DATA (zm(11,i),i=0,1)/22.989768d0, 22.9897692820d0/
DATA (gns(11,i),i=1,1)/4/
DATA (ab(11,i),i=1,1)/100.d0/
c
DATA at(12),gel(12),nmn(12),(mn(12,i),i=1,3)/'Mg',1,3,24,25,26/
DATA (zm(12,i),i=0,3)/24.3050d0, 23.985041698d0, 24.98583698d0,
1 25.98259297d0/
DATA (gns(12,i),i=1,3)/1,6,1/
DATA (ab(12,i),i=1,3)/78.99d0, 10.00d0, 11.01d0/
c
DATA at(13),gel(13),nmn(13),(mn(13,i),i=1,1)/'Al',2,1,27/
DATA (zm(13,i),i=0,1)/26.981539d0, 26.98153853d0/
DATA (gns(13,i),i=1,1)/6/
DATA (ab(13,i),i=1,1)/100.d0/
c
DATA at(14),gel(14),nmn(14),(mn(14,i),i=1,3)/'Si',1,3,28,29,30/
DATA (zm(14,i),i=0,3)/28.0855d0, 27.9769265346d0, 28.9764946649d0,
1 29.973770136d0/
DATA (gns(14,i),i=1,3)/1,2,1/
DATA (ab(14,i),i=1,3)/92.23d0, 4.67d0, 3.10d0/
DATA at(15),gel(15),nmn(15),(mn(15,i),i=1,2)/' P',4,2,26,31/
DATA (zm(15,i),i=0,2)/30.973762d0, 26.01178d0, 30.9737619984d0/
DATA (gns(15,i),i=1,2)/15,2/
DATA (ab(15,i),i=1,2)/0.d0, 100.d0/
c
DATA at(16),gel(16),nmn(16),(mn(16,i),i=1,5)/' S',5,5,27,32,33,
1 34,36/
DATA (zm(16,i),i=0,5)/32.066d0, 27.01883d0, 31.9720711744d0,
1 32.9714589098d0,33.96786700d0, 35.96708071d0/
DATA (gns(16,i),i=1,5)/6,1,4,1,1/
DATA (ab(16,i),i=1,5)/0.d0, 95.02d0, 0.75d0, 4.21d0, 0.02d0/
c
DATA at(17),gel(17),nmn(17),(mn(17,i),i=1,2)/'Cl',4,2,35,37/
DATA (zm(17,i),i=0,2)/35.4527d0, 34.96885268d0, 36.96590260d0/
DATA (gns(17,i),i=1,2)/4,4/
DATA (ab(17,i),i=1,2)/75.77d0, 24.23d0/
c
DATA at(18),gel(18),nmn(18),(mn(18,i),i=1,3)/'Ar',1,3,36,38,40/
DATA (zm(18,i),i=0,3)/39.948d0, 35.967545105d0, 37.96273211d0,
1 39.9623831237d0/
DATA (gns(18,i),i=1,3)/1,1,1/
DATA (ab(18,i),i=1,3)/0.337d0, 0.063d0, 99.600d0/
c
DATA at(19),gel(19),nmn(19),(mn(19,i),i=1,3)/' K',2,3,39,40,41/
DATA (zm(19,i),i=0,3)/39.0983d0, 38.963706486d0, 39.96399817d0,
1 40.961825258d0/
DATA (gns(19,i),i=1,3)/4,9,4/
DATA (ab(19,i),i=1,3)/93.2581d0, 0.0117d0, 6.7302d0/
DATA at(20),gel(20),nmn(20),(mn(20,i),i=1,6)/'Ca',1,6,40,42,43,44,
1 46,48/
DATA (zm(20,i),i=0,6)/40.078d0, 39.962590864d0, 41.95861783d0,
1 42.95876644d0, 43.9554816d0, 45.9536890d0, 47.95252277d0/
DATA (gns(20,i),i=1,6)/1,1,8,1,1,1/
DATA (ab(20,i),i=1,6)/96.941d0, 0.647d0, 0.135d0, 2.086d0,
1 0.004d0, 0.187d0/
c
DATA at(21),gel(21),nmn(21),(mn(21,i),i=1,1)/'Sc',4,1,45/
DATA (zm(21,i),i=0,1)/44.955910d0, 44.9559083d0/
DATA (gns(21,i),i=1,1)/8/
DATA (ab(21,i),i=1,1)/100.d0/
c
DATA at(22),gel(22),nmn(22),(mn(22,i),i=1,5)/'Ti',5,5,46,47,48,49,
1 50/
DATA (zm(22,i),i=0,5)/47.88d0, 45.9526277d0, 46.9517588d0,
1 47.9479420d0, 48.9478657d0, 49.9447869d0/
DATA (gns(22,i),i=1,5)/1,6,1,8,1/
DATA (ab(22,i),i=1,5)/8.0d0, 7.3d0, 73.8d0, 5.5d0, 5.4d0/
c
DATA at(23),gel(23),nmn(23),(mn(23,i),i=1,2)/' V',4,2,50,51/
DATA (zm(23,i),i=0,2)/50.9415d0, 49.9471560d0, 50.9439570d0/
DATA (gns(23,i),i=1,2)/13,8/
DATA (ab(23,i),i=1,2)/0.250d0, 99.750d0/
c
DATA at(24),gel(24),nmn(24),(mn(24,i),i=1,4)/'Cr',7,4,50,52,53,54/
DATA (zm(24,i),i=0,4)/51.9961d0, 49.9460418d0, 51.9405062d0,
1 52.9406481d0, 53.9388792d0/
DATA (gns(24,i),i=1,4)/1,1,4,1/
DATA (ab(24,i),i=1,4)/4.345d0, 83.789d0, 9.501d0, 2.365d0/
c
DATA at(25),gel(25),nmn(25),(mn(25,i),i=1,1)/'Mn',6,1,55/
DATA (zm(25,i),i=0,1)/54.93805d0, 54.938049d0/
DATA (gns(25,i),i=1,1)/6/
DATA (ab(25,i),i=1,1)/100.d0/
c
DATA at(26),gel(26),nmn(26),(mn(26,i),i=1,4)/'Fe',9,4,54,56,57,58/
DATA (zm(26,i),i=0,4)/55.847d0, 53.9396090d0, 55.9349363d0,
1 56.9353928d0, 57.9332744d0/
DATA (gns(26,i),i=1,4)/1,1,2,1/
DATA (ab(26,i),i=1,4)/5.8d0, 91.72d0, 2.2d0, 0.28d0/
c
DATA at(27),gel(27),nmn(27),(mn(27,i),i=1,1)/'Co',10,1,59/
DATA (zm(27,i),i=0,1)/58.93320d0, 58.9331943d0/
DATA (gns(27,i),i=1,1)/8/
DATA (ab(27,i),i=1,1)/100.d0/
c
DATA at(28),gel(28),nmn(28),(mn(28,i),i=1,5)/'Ni',9,5,58,60,61,62,
1 64/
DATA (zm(28,i),i=0,5)/58.69d0, 57.9353424d0, 59.9307859d0,
1 60.9310556d0, 61.9283454d0, 63.9279668d0/
DATA (gns(28,i),i=1,5)/1,1,4,1,1/
DATA (ab(28,i),i=1,5)/68.077d0,26.223d0,1.140d0,3.634d0,0.926d0/
c
DATA at(29),gel(29),nmn(29),(mn(29,i),i=1,2)/'Cu',2,2,63,65/
DATA (zm(29,i),i=0,2)/63.546d0, 62.9295977d0,64.9277897d0/
DATA (gns(29,i),i=1,2)/4,4/
DATA (ab(29,i),i=1,2)/69.17d0, 30.83d0/
c
DATA at(30),gel(30),nmn(30),(mn(30,i),i=1,5)/'Zn',1,5,64,66,67,68,
1 70/
DATA (zm(30,i),i=0,5)/65.40d0, 63.9291420d0, 65.9260338d0,
1 66.9271277d0, 67.9248446d0, 69.9253192d0/
DATA (gns(30,i),i=1,5)/1,1,6,1,1/
DATA (ab(30,i),i=1,5)/48.6d0, 27.9d0, 4.1d0, 18.8d0, 0.6d0/
c
DATA at(31),gel(31),nmn(31),(mn(31,i),i=1,2)/'Ga',2,2,69,71/
DATA (zm(31,i),i=0,2)/69.723d0, 68.9255735d0, 70.9247026d0/
DATA (gns(31,i),i=1,2)/4,4/
DATA (ab(31,i),i=1,2)/60.108d0, 39.892d0/
c
DATA at(32),gel(32),nmn(32),(mn(32,i),i=1,5)/'Ge',1,5,70,72,73,74,
1 76/
DATA (zm(32,i),i=0,5)/72.61d0, 69.9242488d0, 71.92207583d0,
1 72.92345896d0, 73.921177762d0, 75.921402726d0/
DATA (gns(32,i),i=1,5)/1,1,10,1,1/
DATA (ab(32,i),i=1,5)/21.23d0, 27.66d0, 7.73d0, 35.94d0, 7.44d0/
c
DATA at(33),gel(33),nmn(33),(mn(33,i),i=1,1)/'As',4,1,75/
DATA (zm(33,i),i=0,1)/74.92159d0, 74.9215946d0/
DATA (gns(33,i),i=1,1)/4/
DATA (ab(33,i),i=1,1)/100.d0/
c
DATA at(34),gel(34),nmn(34),(mn(34,i),i=1,6)/'Se',5,6,74,76,77,78,
1 80,82/
DATA (zm(34,i),i=0,6)/78.96d0, 73.922475935d0, 75.919213704d0,
1 76.91991415d0, 77.91730928d0, 79.9165218d0, 81.9166995d0/
DATA (gns(34,i),i=1,6)/1,1,2,1,1,1/
DATA (ab(34,i),i=1,6)/0.89d0, 9.36d0, 7.63d0, 23.78d0, 49.61d0,
1 8.73d0/
c
DATA at(35),gel(35),nmn(35),(mn(35,i),i=1,2)/'Br',4,2,79,81/
DATA (zm(35,i),i=0,2)/79.904d0, 78.9183376d0, 80.9162897d0/
DATA (gns(35,i),i=1,2)/4,4/
DATA (ab(35,i),i=1,2)/50.69d0, 49.31d0/
c
DATA at(36),gel(36),nmn(36),(mn(36,i),i=1,6)/'Kr',1,6,78,80,82,83,
1 84,86/
DATA (zm(36,i),i=0,6)/83.80d0, 77.9203649d0, 79.9163781d0,
1 81.9134827d0, 82.9141272d0, 83.911497728d0, 85.910610627d0/
DATA (gns(36,i),i=1,6)/1,1,1,10,1,1/
DATA (ab(36,i),i=1,6)/0.35d0, 2.25d0, 11.6d0, 11.5d0, 57.0d0,
1 17.3d0/
c
DATA at(37),gel(37),nmn(37),(mn(37,i),i=1,2)/'Rb',2,2,85,87/
DATA (zm(37,i),i=0,2)/85.4678d0, 84.911789738d0, 86.909180532d0/
DATA (gns(37,i),i=1,2)/6,4/
DATA (ab(37,i),i=1,2)/72.165d0, 27.835d0/
c
DATA at(38),gel(38),nmn(38),(mn(38,i),i=1,4)/'Sr',1,4,84,86,87,88/
DATA (zm(38,i),i=0,4)/87.62d0, 83.9134191d0, 85.9092606d0,
1 86.9088775d0, 87.9056125d0/
DATA (gns(38,i),i=1,4)/1,1,10,1/
DATA (ab(38,i),i=1,4)/0.56d0, 9.86d0, 7.00d0, 82.58d0/
c
DATA at(39),gel(39),nmn(39),(mn(39,i),i=1,1)/' Y',4,1,89/
DATA (zm(39,i),i=0,1)/88.90585d0, 88.9058403d0/
DATA (gns(39,i),i=1,1)/2/
DATA (ab(39,i),i=1,1)/100.d0/
c
DATA at(40),gel(40),nmn(40),(mn(40,i),i=1,5)/'Zr',5,5,90,91,92,94,
1 96/
DATA (zm(40,i),i=0,5)/91.224d0, 89.9046977d0, 90.9056396d0,
1 91.9050347d0, 93.9063108d0, 95.9082714d0/
DATA (gns(40,i),i=1,5)/1,6,1,1,1/
DATA (ab(40,i),i=1,5)/51.45d0, 11.22d0, 17.15d0, 17.38d0, 2.80d0/
c
DATA at(41),gel(41),nmn(41),(mn(41,i),i=1,1)/'Nb',2,1,93/
DATA (zm(41,i),i=0,1)/92.90638d0, 92.9063730d0/
DATA (gns(41,i),i=1,1)/10/
DATA (ab(41,i),i=1,1)/100.d0/
c
DATA at(42),gel(42),nmn(42),(mn(42,i),i=1,7)/'Mo',7,7,92,94,95,96,
1 97,98,100/
DATA (zm(42,i),i=0,7)/95.94d0, 91.9068080d0, 93.9050849d0,
1 94.9058388d0, 95.9046761d0, 96.9060181d0, 97.9054048d0,
2 99.9074718d0/
DATA (gns(42,i),i=1,7)/1,1,6,1,6,1,1/
DATA (ab(42,i),i=1,7)/14.84d0, 9.25d0, 15.92d0, 16.68d0, 9.55d0,
1 24.13d0, 9.63d0/
c
DATA at(43),gel(43),nmn(43),(mn(43,i),i=1,1)/'Tc',6,1,98/
DATA (zm(43,i),i=0,1)/97.907215d0, 97.907212d0/
DATA (gns(43,i),i=1,1)/13/
DATA (ab(43,i),i=1,1)/100.d0/
c
DATA at(44),gel(44),nmn(44),(mn(44,i),i=1,7)/'Ru',11,7,96,98,99,
1 100,101,102,104/
DATA (zm(44,i),i=0,7)/101.07d0, 95.9075903d0, 97.905287d0,
1 98.9059341d0, 99.9042143d0, 100.9055769d0, 101.9043441d0,
2 103.9054275d0/
DATA (gns(44,i),i=1,7)/1,1,6,1,6,1,1/
DATA (ab(44,i),i=1,7)/5.52d0, 1.88d0, 12.7d0, 12.6d0, 17.0d0,
1 31.6d0, 18.7d0/
c
DATA at(45),gel(45),nmn(45),(mn(45,i),i=1,1)/'Rh',10,1,103/
DATA (zm(45,i),i=0,1)/102.90550d0, 102.9054980d0/
DATA (gns(45,i),i=1,1)/2/
DATA (ab(45,i),i=1,1)/100.d0/
c
DATA at(46),gel(46),nmn(46),(mn(46,i),i=1,6)/'Pd',1,6,102,104,105,
1 106,108,110/
DATA (zm(46,i),i=0,6)/106.42d0, 101.9056022d0, 103.9040305d0,
1 104.9050796d0, 105.9034804d0, 107.9038916d0, 109.9051722d0/
DATA (gns(46,i),i=1,6)/1,1,6,1,1,1/
DATA (ab(46,i),i=1,6)/1.02d0, 11.14d0, 22.33d0, 27.33d0, 26.46d0,
1 11.72d0/
c
DATA at(47),gel(47),nmn(47),(mn(47,i),i=1,2)/'Ag',2,2,107,109/
DATA (zm(47,i),i=0,2)/107.8682d0, 106.9050916d0, 108.9047553d0/
DATA (gns(47,i),i=1,2)/2,2/
DATA (ab(47,i),i=1,2)/51.839d0, 48.161d0/
c
DATA at(48),gel(48),nmn(48),(mn(48,i),i=1,8)/'Cd',1,8,106,108,110,
1 111,112,113,114,116/
DATA (zm(48,i),i=0,8)/112.411d0, 105.9064599d0, 107.9041834d0,
1 109.9030066d0, 110.9041829d0, 111.9027629d0, 112.9044081d0,
2 113.9033651d0, 115.90476315d0/
DATA (gns(48,i),i=1,8)/1,1,1,2,1,2,1,1/
DATA (ab(48,i),i=1,8)/1.25d0, 0.89d0, 12.49d0, 12.80d0, 24.13d0,
1 12.22d0, 28.73d0, 7.49d0/
c
DATA at(49),gel(49),nmn(49),(mn(49,i),i=1,2)/'In',2,2,113,115/
DATA (zm(49,i),i=0,2)/114.818d0, 112.9040618d0, 114.903878776d0/
DATA (gns(49,i),i=1,2)/10,10/
DATA (ab(49,i),i=1,2)/4.3d0, 95.7d0/
c
DATA at(50),gel(50),nmn(50),(mn(50,i),i=1,10)/'Sn',1,10,112,114,
1 115,116,117,118,119,120,122,124/
DATA (zm(50,i),i=0,10)/118.710d0, 111.9048239d0, 113.9027827d0,
1 114.903344699d0, 115.90174280d0, 116.9029540d0, 117.9016066d0,
2 118.9033112d0, 119.9022016d0, 121.9034438d0, 123.9052766d0/
DATA (gns(50,i),i=1,10)/1,1,2,1,2,1,2,1,1,1/
DATA (ab(50,i),i=1,10)/0.97d0, 0.65d0, 0.34d0, 14.53d0, 7.68d0,
1 24.23d0, 8.59d0, 32.59d0, 4.63d0, 5.79d0/
c
DATA at(51),gel(51),nmn(51),(mn(51,i),i=1,2)/'Sb',4,2,121,123/
DATA (zm(51,i),i=0,2)/121.757d0, 120.903812d0, 122.9042132d0/
DATA (gns(51,i),i=1,2)/6,8/
DATA (ab(51,i),i=1,2)/57.36d0, 42.64d0/
c
DATA at(52),gel(52),nmn(52),(mn(52,i),i=1,8)/'Te',5,8,120,122,123,
1 124,125,126,128,130/
DATA (zm(52,i),i=0,8)/127.60d0, 119.904059d0, 121.9030435d0,
1 122.9042698d0, 123.9028171d0, 124.9044299d0, 125.9033109d0,
2 127.9044613d0, 129.906222749d0/
DATA (gns(52,i),i=1,8)/1,1,2,1,2,1,1,1/
DATA (ab(52,i),i=1,8)/0.096d0, 2.603d0, 0.908d0, 4.816d0,
1 7.139d0, 18.95d0, 31.69d0, 33.80d0/
c
DATA at(53),gel(53),nmn(53),(mn(53,i),i=1,2)/' I',4,2,127,129/
DATA (zm(53,i),i=0,2)/126.90447d0, 126.904472d0, 128.904984d0/
DATA (gns(53,i),i=1,2)/6,8/
DATA (ab(53,i),i=1,2)/100.d0,0.d0/
c
DATA at(54),gel(54),nmn(54),(mn(54,i),i=1,9)/'Xe',1,9,124,126,128,
1 129,130,131,132,134,136/
DATA (zm(54,i),i=0,9)/131.29d0, 123.9058920d0, 125.904298d0,
1 127.9035310d0, 128.904780861d0,129.903509350d0,130.90508406d0,
2 131.904155086d0, 133.9053947d0, 135.907214484d0/
DATA (gns(54,i),i=1,9)/1,1,1,2,1,4,1,1,1/
DATA (ab(54,i),i=1,9)/0.10d0, 0.09d0, 1.91d0, 26.4d0, 4.1d0,
1 21.2d0, 26.9d0, 10.4d0, 8.9d0/
c
DATA at(55),gel(55),nmn(55),(mn(55,i),i=1,1)/'Cs',2,1,133/
DATA (zm(55,i),i=0,1)/132.90543d0, 132.905451961d0/
DATA (gns(55,i),i=1,1)/8/
DATA (ab(55,i),i=1,1)/100.d0/
c
DATA at(56),gel(56),nmn(56),(mn(56,i),i=1,7)/'Ba',1,7,130,132,134,
1 135,136,137,138/
DATA (zm(56,i),i=0,7)/137.327d0, 129.9063207d0, 131.9050611d0,
1 133.90450818d0, 134.90568838d0, 135.90457573d0, 136.9058271d0,
2 137.9052470d0/
DATA (gns(56,i),i=1,7)/1,1,1,4,1,4,1/
DATA (ab(56,i),i=1,7)/0.106d0, 0.101d0, 2.417d0, 6.592d0,
1 7.854d0, 11.23d0, 71.70d0/
c
DATA at(57),gel(57),nmn(57),(mn(57,i),i=1,2)/'La',4,2,138,139/
DATA (zm(57,i),i=0,2)/138.9055d0, 137.907115d0, 138.9063563d0/
DATA (gns(57,i),i=1,2)/11,8/
DATA (ab(57,i),i=1,2)/0.0902d0, 99.9098d0/
c
DATA at(58),gel(58),nmn(58),(mn(58,i),i=1,4)/'Ce',9,4,136,138,140,
1 142/
DATA (zm(58,i),i=0,4)/140.115d0, 135.9071292d0, 137.905991d0,
1 139.9054431d0, 141.9092504d0/
DATA (gns(58,i),i=1,4)/1,1,1,1/
DATA (ab(58,i),i=1,4)/0.19d0, 0.25d0, 88.48d0, 11.08d0/
c
DATA at(59),gel(59),nmn(59),(mn(59,i),i=1,1)/'Pr',10,1,141/
DATA (zm(59,i),i=0,1)/140.90765d0, 140.9076576d0/
DATA (gns(59,i),i=1,1)/6/
DATA (ab(59,i),i=1,1)/100.d0/
c
DATA at(60),gel(60),nmn(60),(mn(60,i),i=1,7)/'Nd',9,7,142,143,144,
1 145,146,148,150/
DATA (zm(60,i),i=0,7)/144.24d0, 141.9077290d0, 142.9098200d0,
1 143.9100930d0, 144.9125793d0, 145.9131226d0, 147.9168993d0,
2 149.9209022d0/
DATA (gns(60,i),i=1,7)/1,8,1,8,1,1,1/
DATA (ab(60,i),i=1,7)/27.13d0, 12.18d0, 23.80d0, 8.30d0, 17.19d0,
1 5.76d0, 5.64d0/
c
DATA at(61),gel(61),nmn(61),(mn(61,i),i=1,1)/'Pm',6,1,145/
DATA (zm(61,i),i=0,1)/144.912743d0, 144.912756d0/
DATA (gns(61,i),i=1,1)/6/
DATA (ab(61,i),i=1,1)/100.d0/
c
DATA at(62),gel(62),nmn(62),(mn(62,i),i=1,7)/'Sm',1,7,144,147,148,
1 149,150,152,154/
DATA (zm(62,i),i=0,7)/150.36d0, 143.9120065d0, 146.9149044d0,
1 147.9148292d0, 148.9171921d0, 149.9172829d0, 151.9197397d0,
2 153.9222169d0/
DATA (gns(62,i),i=1,7)/1,8,1,8,1,1,1/
DATA (ab(62,i),i=1,7)/3.1d0, 15.0d0, 11.3d0, 13.8d0, 7.4d0,
1 26.7d0, 22.7d0/
c
DATA at(63),gel(63),nmn(63),(mn(63,i),i=1,2)/'Eu',8,2,151,153/
DATA (zm(63,i),i=0,2)/151.965d0, 150.9198578d0, 152.9212380d0/
DATA (gns(63,i),i=1,2)/6,6/
DATA (ab(63,i),i=1,2)/47.8d0, 52.2d0/
c
DATA at(64),gel(64),nmn(64),(mn(64,i),i=1,7)/'Gd',5,7,152,154,155,
1 156,157,158,160/
DATA (zm(64,i),i=0,7)/157.25d0, 151.9197995d0, 153.9208741d0,
1 154.9226305d0, 155.9221312d0, 156.9239686d0, 157.9241123d0,
2 159.9270624d0/
DATA (gns(64,i),i=1,7)/1,1,4,1,4,1,1/
DATA (ab(64,i),i=1,7)/0.20d0, 2.18d0, 14.80d0, 20.47d0, 15.65d0,
1 24.84d0, 21.86d0/
c
DATA at(65),gel(65),nmn(65),(mn(65,i),i=1,1)/'Tb',16,1,159/
DATA (zm(65,i),i=0,1)/158.92534d0, 158.9253547d0/
DATA (gns(65,i),i=1,1)/4/
DATA (ab(65,i),i=1,1)/100.d0/
c
DATA at(66),gel(66),nmn(66),(mn(66,i),i=1,7)/'Dy',17,7,156,158,
1 160,161,162,163,164/
DATA (zm(66,i),i=0,7)/162.50d0, 155.9242847d0, 157.924416d0,
1 159.9252046d0, 160.9269405d0, 161.9268056d0, 162.9287383d0,
2 163.9291819d0/
DATA (gns(66,i),i=1,7)/1,1,1,6,1,6,1/
DATA (ab(66,i),i=1,7)/0.06d0, 0.10d0, 2.34d0, 18.9d0, 25.5d0,
1 24.9d0, 28.2d0/
c
DATA at(67),gel(67),nmn(67),(mn(67,i),i=1,1)/'Ho',16,1,165/
DATA (zm(67,i),i=0,1)/164.93032d0, 164.9303288d0/
DATA (gns(67,i),i=1,1)/8/
DATA (ab(67,i),i=1,1)/100.d0/
DATA at(68),gel(68),nmn(68),(mn(68,i),i=1,6)/'Er',13,6,162,164,
1 166,167,168,170/
DATA (zm(68,i),i=0,6)/167.26d0, 161.9287884d0, 163.9292088d0,
1 165.9302995d0, 166.9320546d0, 167.9323767d0, 169.9354702d0/
DATA (gns(68,i),i=1,6)/1,1,1,8,1,1/
DATA (ab(68,i),i=1,6)/0.14d0, 1.61d0, 33.6d0, 22.95d0, 26.8d0,
1 14.9d0/
c
DATA at(69),gel(69),nmn(69),(mn(69,i),i=1,1)/'Tm',8,1,169/
DATA (zm(69,i),i=0,1)/168.93421d0, 168.9342179d0/
DATA (gns(69,i),i=1,1)/2/
DATA (ab(69,i),i=1,1)/100.d0/
c
DATA at(70),gel(70),nmn(70),(mn(70,i),i=1,7)/'Yb',1,7,168,170,171,
1 172,173,174,176/
DATA (zm(70,i),i=0,7)/173.04d0, 167.9338896d0, 169.9347664d0,
1 170.9363302d0, 171.9363859d0, 172.9382151d0, 173.9388664d0,
2 175.9425764d0/
DATA (gns(70,i),i=1,7)/1,1,2,1,6,1,1/
DATA (ab(70,i),i=1,7)/0.13d0, 3.05d0, 14.3d0, 21.9d0, 16.12d0,
1 31.8d0, 12.7d0/
c
DATA at(71),gel(71),nmn(71),(mn(71,i),i=1,2)/'Lu',4,2,175,176/
DATA (zm(71,i),i=0,2)/174.967d0, 174.9407752d0, 175.9426897d0/
DATA (gns(71,i),i=1,2)/6,15/
DATA (ab(71,i),i=1,2)/97.41d0, 2.59d0/
c
DATA at(72),gel(72),nmn(72),(mn(72,i),i=1,6)/'Hf',5,6,174,176,177,
1 178,179,180/
DATA (zm(72,i),i=0,6)/178.49d0, 173.9400461d0, 175.9414076d0,
1 176.9432277d0, 177.9437058d0, 178.9458232d0, 179.9465570d0/
DATA (gns(72,i),i=1,6)/1,1,8,1,10,1/
DATA (ab(72,i),i=1,6)/0.162d0, 5.206d0, 18.606d0, 27.297d0,
1 13.629d0, 35.100d0/
c
DATA at(73),gel(73),nmn(73),(mn(73,i),i=1,2)/'Ta',4,2,180,181/
DATA (zm(73,i),i=0,2)/180.9479d0, 179.9474648d0, 180.9479958d0/
DATA (gns(73,i),i=1,2)/17,8/
DATA (ab(73,i),i=1,2)/0.012d0, 99.988d0/
c
DATA at(74),gel(74),nmn(74),(mn(74,i),i=1,5)/' W',1,5,180,182,183,
1 184,186/
DATA (zm(74,i),i=0,5)/183.84d0, 179.9467108d0, 181.9482039d0,
1 182.9502227d0, 183.9509309d0, 185.9543628d0/
DATA (gns(74,i),i=1,5)/1,1,2,1,1/
DATA (ab(74,i),i=1,5)/0.13d0, 26.3d0, 14.3d0, 30.67d0, 28.6d0/
c
DATA at(75),gel(75),nmn(75),(mn(75,i),i=1,2)/'Re',6,2,185,187/
DATA (zm(75,i),i=0,2)/186.207d0, 184.9529545d0, 186.9557501d0/
DATA (gns(75,i),i=1,2)/6,6/
DATA (ab(75,i),i=1,2)/37.40d0, 62.60d0/
c
DATA at(76),gel(76),nmn(76),(mn(76,i),i=1,7)/'Os',9,7,184,186,187,
1 188,189,190,192/
DATA (zm(76,i),i=0,7)/190.23d0, 183.9524885d0, 185.9538350d0,
1 186.9557474d0, 187.9558352d0, 188.9581442d0, 189.9584437d0,
2 191.9614770d0/
DATA (gns(76,i),i=1,7)/1,1,2,1,4,1,1/
DATA (ab(76,i),i=1,7)/0.02d0, 1.58d0, 1.6d0, 13.3d0, 16.1d0,
1 26.4d0, 41.0d0/
c
DATA at(77),gel(77),nmn(77),(mn(77,i),i=1,2)/'Ir',10,2,191,193/
DATA (zm(77,i),i=0,2)/192.22d0, 190.9605893d0, 192.9629216d0/
DATA (gns(77,i),i=1,2)/4,4/
DATA (ab(77,i),i=1,2)/37.3d0, 62.7d0/
c
c
DATA at(78),gel(78),nmn(78),(mn(78,i),i=1,6)/'Pt',7,6,190,192,194,
1 195,196,198/
DATA (zm(78,i),i=0,6)/195.08d0, 189.959930d0, 191.961039d0,
1 193.9626809d0, 194.9647917d0, 195.9649521d0, 197.9678949d0/
DATA (gns(78,i),i=1,6)/1,1,1,2,1,1/
DATA (ab(78,i),i=1,6)/0.01d0,0.79d0,32.9d0,33.8d0,25.3d0,7.2d0/
c
DATA at(79),gel(79),nmn(79),(mn(79,i),i=1,1)/'Au',2,1,197/
DATA (zm(79,i),i=0,1)/196.96654d0, 196.9665688d0/
DATA (gns(79,i),i=1,1)/4/
DATA (ab(79,i),i=1,1)/100.d0/
c
DATA at(80),gel(80),nmn(80),(mn(80,i),i=1,7)/'Hg',1,7,196,198,199,
1 200,201,202,204/
DATA (zm(80,i),i=0,7)/200.59d0, 195.965833d0, 197.9667686d0,
1 198.9682806d0, 199.9683266d0, 200.9703028d0, 201.9706434d0,
2 203.9734940d0/
DATA (gns(80,i),i=1,7)/1,1,2,1,4,1,1/
DATA (ab(80,i),i=1,7)/0.15d0, 9.97d0, 16.87d0, 23.10d0, 13.18d0,
1 29.86d0, 6.87d0/
c
DATA at(81),gel(81),nmn(81),(mn(81,i),i=1,2)/'Tl',2,2,203,205/
DATA (zm(81,i),i=0,2)/204.3833d0, 202.9723446d0, 204.9744278d0/
DATA (gns(81,i),i=1,2)/2,2/
DATA (ab(81,i),i=1,2)/29.524d0, 70.476d0/
c
DATA at(82),gel(82),nmn(82),(mn(82,i),i=1,4)/'Pb',1,4,204,206,207,
1 208/
DATA (zm(82,i),i=0,4)/207.2d0, 203.9730440d0, 205.9744657d0,
1 206.9758973d0, 207.9766525d0/
DATA (gns(82,i),i=1,4)/1,1,2,1/
DATA (ab(82,i),i=1,4)/1.4d0, 24.1d0, 22.1d0, 52.4d0/
c
DATA at(83),gel(83),nmn(83),(mn(83,i),i=1,1)/'Bi',4,1,209/
DATA (zm(83,i),i=0,1)/208.98037d0, 208.9803991d0/
DATA (gns(83,i),i=1,1)/10/
DATA (ab(83,i),i=1,1)/100.d0/
c
DATA at(84),gel(84),nmn(84),(mn(84,i),i=1,1)/'Po',5,1,209/
DATA (zm(84,i),i=0,1)/208.982404d0, 208.9824308d0/
DATA (gns(84,i),i=1,1)/2/
DATA (ab(84,i),i=1,1)/100.d0/
c
DATA at(85),gel(85),nmn(85),(mn(85,i),i=1,1)/'At',-1,1,210/
DATA (zm(85,i),i=0,1)/209.987126d0, 209.987148d0/
DATA (gns(85,i),i=1,1)/11/
DATA (ab(85,i),i=1,1)/100.d0/
c
DATA at(86),gel(86),nmn(86),(mn(86,i),i=1,1)/'Rn',1,1,222/
DATA (zm(86,i),i=0,1)/222.017571d0, 222.0175782d0/
DATA (gns(86,i),i=1,1)/1/
DATA (ab(86,i),i=1,1)/100.d0/
c
DATA at(87),gel(87),nmn(87),(mn(87,i),i=1,1)/'Fr',-1,1,223/
DATA (zm(87,i),i=0,1)/223.019733d0, 223.0197360d0/
DATA (gns(87,i),i=1,1)/4/
DATA (ab(87,i),i=1,1)/100.d0/
c
DATA at(88),gel(88),nmn(88),(mn(88,i),i=1,1)/'Ra',1,1,226/
DATA (zm(88,i),i=0,1)/226.025403d0, 226.0254103d0/
DATA (gns(88,i),i=1,1)/1/
DATA (ab(88,i),i=1,1)/100.d0/
c
DATA at(89),gel(89),nmn(89),(mn(89,i),i=1,1)/'Ac',4,1,227/
DATA (zm(89,i),i=0,1)/227.027750d0, 227.0277523d0/
DATA (gns(89,i),i=1,1)/4/
DATA (ab(89,i),i=1,1)/100.d0/
c
DATA at(90),gel(90),nmn(90),(mn(90,i),i=1,1)/'Th',-1,1,232/
DATA (zm(90,i),i=0,1)/232.038d0, 232.0380558d0/
DATA (gns(90,i),i=1,1)/1/
DATA (ab(90,i),i=1,1)/100.d0/
c
DATA at(91),gel(91),nmn(91),(mn(91,i),i=1,1)/'Pa',-1,1,231/
DATA (zm(91,i),i=0,1)/231.03588d0, 231.0358842d0/
DATA (gns(91,i),i=1,1)/4/
DATA (ab(91,i),i=1,1)/100.d0/
c
DATA at(92),gel(92),nmn(92),(mn(92,i),i=1,4)/' U',-1,4,233,234,
1 235,238/
DATA (zm(92,i),i=0,4)/238.0289d0, 233.0396355d0, 234.0409523d0,
1 235.0439301d0, 238.0507884d0/
DATA (gns(92,i),i=1,4)/6,1,8,1/
DATA (ab(92,i),i=1,4)/0.d0, 0.0055d0, 0.7200d0, 99.2745d0/
c
DATA at(93),gel(93),nmn(93),(mn(93,i),i=1,1)/'Np',-1,1,237/
DATA (zm(93,i),i=0,1)/237.0481678d0, 237.0481736d0/
DATA (gns(93,i),i=1,1)/6/
DATA (ab(93,i),i=1,1)/100.d0/
c
DATA at(94),gel(94),nmn(94),(mn(94,i),i=1,1)/'Pu',-1,1,244/
DATA (zm(94,i),i=0,1)/244.064199d0, 244.064205d0/
DATA (gns(94,i),i=1,1)/1/
DATA (ab(94,i),i=1,1)/100.d0/
c
DATA at(95),gel(95),nmn(95),(mn(95,i),i=1,1)/'Am',-1,1,243/
DATA (zm(95,i),i=0,1)/243.061375d0, 243.0613815d0/
DATA (gns(95,i),i=1,1)/6/
DATA (ab(95,i),i=1,1)/100.d0/
c
DATA at(96),gel(96),nmn(96),(mn(96,i),i=1,1)/'Cm',-1,1,247/
DATA (zm(96,i),i=0,1)/247.070347d0, 247.070354d0/
DATA (gns(96,i),i=1,1)/10/
DATA (ab(96,i),i=1,1)/100.d0/
c
DATA at(97),gel(97),nmn(97),(mn(97,i),i=1,1)/'Bk',-1,1,247/
DATA (zm(97,i),i=0,1)/247.070300d0, 247.070307d0/
DATA (gns(97,i),i=1,1)/4/
DATA (ab(97,i),i=1,1)/100.d0/
c
DATA at(98),gel(98),nmn(98),(mn(98,i),i=1,1)/'Cf',-1,1,251/
DATA (zm(98,i),i=0,1)/251.079580d0, 251.079589d0/
DATA (gns(98,i),i=1,1)/2/
DATA (ab(98,i),i=1,1)/100.d0/
c
DATA at(99),gel(99),nmn(99),(mn(99,i),i=1,1)/'Es',-1,1,252/
DATA (zm(99,i),i=0,1)/252.082944d0, 252.082980d0/
DATA (gns(99,i),i=1,1)/11/
DATA (ab(99,i),i=1,1)/100.d0/
c
DATA at(100),gel(100),nmn(100),(mn(100,i),i=1,1)/'Fm',-1,1,257/
DATA (zm(100,i),i=0,1)/257.095099d0, 257.095106d0/
DATA (gns(100,i),i=1,1)/10/
DATA (ab(100,i),i=1,1)/100.d0/
c
DATA at(101),gel(101),nmn(101),(mn(101,i),i=1,1)/'Md',-1,1,258/
DATA (zm(101,i),i=0,1)/258.09857d0, 258.098431d0/
DATA (gns(101,i),i=1,1)/17/
DATA (ab(101,i),i=1,1)/100.d0/
c
DATA at(102),gel(102),nmn(102),(mn(102,i),i=1,1)/'No',-1,1,259/
DATA (zm(102,i),i=0,1)/259.100931d0, 259.101030d0/
DATA (gns(102,i),i=1,1)/10/
DATA (ab(102,i),i=1,1)/100.d0/
c
DATA at(103),gel(103),nmn(103),(mn(103,i),i=1,1)/'Lr',-1,1,260/
DATA (zm(103,i),i=0,1)/260.105320d0, 260.105510d0/
DATA (gns(103,i),i=1,1)/-1/
DATA (ab(103,i),i=1,1)/100.d0/
c
DATA at(104),gel(104),nmn(104),(mn(104,i),i=1,1)/'Rf',-1,1,261/
DATA (zm(104,i),i=0,1)/261.10869d0, 261.108770d0/
DATA (gns(104,i),i=1,1)/-1/
DATA (ab(104,i),i=1,1)/100.d0/
c
DATA at(105),gel(105),nmn(105),(mn(105,i),i=1,1)/'Db',-1,1,262/
DATA (zm(105,i),i=0,1)/262.11376d0, 262.114070d0/
DATA (gns(105,i),i=1,1)/-1/
DATA (ab(105,i),i=1,1)/100.d0/
c
DATA at(106),gel(106),nmn(106),(mn(106,i),i=1,1)/'Sg',-1,1,263/
DATA (zm(106,i),i=0,1)/263.11822d0, 263.118290d0/
DATA (gns(106,i),i=1,1)/-1/
DATA (ab(106,i),i=1,1)/100.d0/
c
DATA at(107),gel(107),nmn(107),(mn(107,i),i=1,1)/'Bh',-1,1,262/
DATA (zm(107,i),i=0,1)/262.12293d0, 262.122970d0/
DATA (gns(107,i),i=1,1)/-1/
DATA (ab(107,i),i=1,1)/100.d0/
c
DATA at(108),gel(108),nmn(108),(mn(108,i),i=1,1)/'Hs',-1,1,265/
DATA (zm(108,i),i=0,1)/265.13016d0, 265.129793d0/
DATA (gns(108,i),i=1,1)/-1/
DATA (ab(108,i),i=1,1)/100.d0/
c
DATA at(109),gel(109),nmn(109),(mn(109,i),i=1,1)/'Mt',-1,1,266/
DATA (zm(109,i),i=0,1)/266.13764d0, 266.137370d0/
DATA (gns(109,i),i=1,1)/-1/
DATA (ab(109,i),i=1,1)/100.d0/
c
IF((IAN.LT.0).OR.(IAN.GT.109)) THEN
MASS= 0.d0
NAME= 'XX'
IMN= 0
WRITE(6,601) IAN
RETURN
ELSE
NAME= AT(IAN)
ENDIF
IF((IAN.EQ.1).AND.(IMN.GT.1)) THEN
c** Special case: insert common name for deuterium or tritium
IF(IMN.EQ.2) NAME=' D'
IF(IMN.EQ.3) NAME=' T'
ENDIF
IF((IAN.EQ.0).AND.(IMN.GT.1)) THEN
IF(IMN.EQ.2) NAME=' d'
IF(IMN.EQ.3) NAME=' t'
ENDIF
GELGS= GEL(IAN)
MASS= -1.d0
DGNS= -1
ABUND = -1.d0
DO I= 1,NMN(IAN)
if(i.gt.15) write(6,606) ian,imn,nmn(ian)
IF(IMN.EQ.MN(IAN,I)) THEN
MASS= ZM(IAN,I)
DGNS= gns(IAN,I)
ABUND = AB(IAN,I)
ENDIF
ENDDO
IF(MASS.LT.0.d0) THEN
MASS= ZM(IAN,0)
IF(IMN.NE.0) WRITE(6,602) AT(IAN),IMN
IMN= 0
ENDIF
RETURN
601 FORMAT(' *** MASSES Data base does not include Atomic Number=',i4)
602 FORMAT(' *** MASSES Data base does not include ',A2,'(',i3,
1 '), so use average atomic mass.')
606 format(/' *** ERROR *** called MASSES for atom with AN=',I4,
1 ' MN=',I4,'n(MN)=',I4)
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c***********************************************************************
SUBROUTINE WGHT(NGP)
c** Subroutine to generate points (XG & X2) and weights (WG & W2) for
c both regular (XG & WG) and singular 1/sqrt(1-X) (X2 & W2)
c Gaussian integration for NGP=8 and 16 (for both) and for NGP=7 & 32
c (for regular Gaussian only)
c-----------------------------------------------------------------------
INTEGER I,J,NGP,NGPH
REAL*8 x7(4),w7(4),aa(4),bb(4),A(8),B(8),XX(16),WX(16)
c** Common block for quadrature weights & points
REAL*8 XG(32),WG(32),X2(16),W2(16)
COMMON /GWGHT/XG,WG,X2,W2
c
data x7/0.949107912342759d0, 0.741531185599394d0,
1 0.405845151377397d0, 0.d0/,
2 w7/0.129484966168870d0, 0.279705391489277d0,
3 0.381830550505119d0, 0.417959183673469d0/
data aa/0.960289856497536d0, 0.796666477413627d0,
1 0.525532409916329d0, 0.183434642495650d0/,
2 bb/0.101228536290376d0, 0.222381034453374d0,
3 0.313706645877887d0, 0.362683783378362d0/
DATA A/0.989400934991649932596154D0,0.944575023073232576077988D0,
1 0.865631202387831743880468D0,0.755404408355003033895101D0,
2 0.617876244402643748446672D0,0.458016777657227386342420D0,
3 0.281603550779258913230461D0,0.095012509837637440185319D0/,
4 B/0.02715245941175409485178D0,0.06225352393864789286284D0,
5 0.09515851168249278480993D0,0.12462897125553387205248D0,
6 0.14959598881657673208150D0,0.16915651939500253818931D0,
7 0.18260341504492358886676D0,0.18945061045506849628540D0/
DATA XX/0.997263861849481563544981D0,0.985611511545268335400175D0,
1 0.964762255587506430773812D0,0.934906075937739689170919D0,
2 0.896321155766052123965307D0,0.849367613732569970133693D0,
3 0.794483795967942406963097D0,0.732182118740289680387427D0,
4 0.663044266930215200975115D0,0.587715757240762329040746D0,
5 0.506899908932229390023747D0,0.421351276130635345364120D0,
6 0.331868602282127649779917D0,0.239287362252137074544603D0,
7 0.144471961582796493485186D0,0.048307665687738316234813D0/
8 ,WX/0.00701861000947009660041D0,0.01627439473090567060517D0,
9 0.02539206530926205945575D0,0.03427386291302143310269D0,
1 0.04283589802222668065688D0,0.05099805926237617619616D0,
2 0.05868409347853554714528D0,0.06582222277636184683765D0,
3 0.07234579410884850622540D0,0.07819389578707030647174D0,
4 0.08331192422694675522220D0,0.08765209300440381114277D0,
5 0.09117387869576388471287D0,0.09384439908080456563918D0,
6 0.09563872007927485941908D0,0.09654008851472780056676D0/
c** For simple Gaussian case
NGPH= (NGP+1)/2
J= NGP+1
if(ngp.eq.7) then
DO I= 1,4
J= J-1
XG(I)= -x7(I)
WG(I)= w7(I)
XG(J)= x7(I)
WG(J)= w7(I)
ENDDO
c%% WRITE(6,601) NGP,(XG(I),WG(I),I= 1,NGP)
c*** For 1/SQRT(E-V) case
c%% WRITE(6,603) (XX(I),WX(I),I= 1,NGP)
write(6,606) ngp
go to 99
endif
if(ngp.eq.8) then
DO I= 1,4
J= J-1
XG(I)= -aa(I)
WG(I)= bb(I)
XG(J)= aa(I)
WG(J)= bb(I)
ENDDO
c%% WRITE(6,601) NGP,(XG(I),WG(I),I= 1,NGP)
c*** For 1/SQRT(E-V) case
c%% WRITE(6,603) (XX(I),WX(I),I= 1,NGP)
DO I= 1,8
X2(I)= 1.d0-A(I)**2
c** Include SQRT(1-X) factor in weight, not integrand.
W2(I)= 2.d0*B(I)*A(I)
ENDDO
c%% WRITE(6,604) NGP,(X2(I),W2(I),I=1,NGP)
go to 99
ENDIF
if(ngp.eq.16) then
DO I=1,8
J=J-1
XG(I)=-A(I)
WG(I)=B(I)
XG(J)=A(I)
WG(J)=B(I)
ENDDO
c%% WRITE(6,601) NGP,(XG(I),WG(I),I=1,NGP)
c*** For 1/SQRT(E-V) case
c%% WRITE(6,603) (XX(I),WX(I),I=1,NGP)
DO I=1,16
X2(I)=1.d0-XX(I)**2
c** Include SQRT(1-X) factor in weight, not integrand.
W2(I)=2.d0*WX(I)*XX(I)
ENDDO
c%% WRITE(6,604) NGP,(X2(I),W2(I),I=1,NGP)
go to 99
endif
if(ngp.eq.32) then
DO I=1,NGPH
J=J-1
XG(I)=-XX(I)
WG(I)=WX(I)
XG(J)=XX(I)
WG(J)=WX(I)
ENDDO
c%% WRITE(6,601) NGP,(XG(I),WG(I),I=1,NGP)
write(6,606) ngp
go to 99
endif
write(6,608) ngp
99 RETURN
c 601 FORMAT(//2X,'Points and weights for simple Gaussian',i3,
c 1 '-point quadrature'/(1x,2(F25.20,F23.20)))
c 603 FORMAT(/' Simple Gaussian points and weights used for generating G
c 1auss-Mehler points and weights:'/(5X,2(F25.20,F23.20)))
c 604 FORMAT(/' Points and weights for singular integrand',i3,
c 1 '-point quadratures'/(1x,2(F25.20,F23.20)))
606 format(/' NOTE: cannot generate singular integrand Gaussian point
1s and weights for NGP =',i3)
608 format(/' NOTE: cannot generate Gaussian points & weights for',
1 ' NGP =',i3)
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12