Ideal Gas Thermodynamic Functions For Water

The calculation of ideal gas thermodynamic properties for steam to 10,000 K is examined. Centrifugal effects are included using spectroscopic data for the lowest vibrational levels, with extension to higher bending levels based on estimates from a bending model. Modifications are examined for rotational and vibrational cut-off effects. Uncertainties in obtaining a suitably regularized representation of energy versus bond stretching vibration in approaching the dissociation energy region appear relevant to the reliability of the extrapolation.


Introduction
A calculation made in 1979 of the ideal gas thermodynamic properties for water [1]' is the background for this study on improving the extrapolation to higher temperatures. The procedure under consideration is to calculate the internal partition function and its first two temperature derivatives, using direct summation over vibrational levels.
Direct summation was also used for the 1979 work. It was used not only for vibrational levels, but for rotational ones as well for temperatures up to 230 K. In the earlier work of Friedman and Haar [2], sums over vibrational levels were computed according to a procedure similar to that of Mayer and Mayer [3] for diatomic molecules. Such a formula-About the Author: Harold W, Woolley, a physicist, retired from NBS some years ago but continues his affiliation with the Bureau as a guest worker in its Thermophysics Division, Center for Chemical Engineering. ' tion is in principle a low temperature form, analogous to a power series in increasing powers of temperature, and as such encounters questions of convergence in extension to very high temperature.
While it is not too difficult in such an approach to include effects of lowest order anharmonicities in their first power contribution, the inclusion of higher power contributions of even the lowest order of anharmonicities including cross product effects and rotational dependence is more complex. The extension to cover still higher order terms such as might satisfactorily describe the approach to the energy of dissociation would involve too many terms of cross product power type for easy confidence in the adequacy of their enumeration. There is in fact no formal end to the number of possible orders of the expansion, and each succeeding order would be enormously more complex than the one before it.
In the presentation that follows, the problem of the myriads of correction terms is circumvented by using direct summation over the three vibrational Journal of Research of the National Bureau of Standards quantum numbers for this triatomic molecule. The problems which remain are the physical ones of obtaining reliable energies for the vibrational levels and reliable rotational parameters, including the way in which the rotational energy behaves in the extension to high rotational quantum numbers.
1) As a first step, this last item is discussed with a study of an empirical bending model, including its implications as to bending energy.
2) The next topic to be covered is the implied centrifugal thermodynamic effect based on empirical constants fitting the spectroscopically-determined rotational levels as represented by a Watson Hamiltonian [4]. This applies to eight low lying vibrational levels, extending up to 5331 cm '.
3) Following this, the problem of extending the representation of rotational details to high energy and high temperature is considered, including the termination of rotational levels due to lack of bond stability. 4) Finally, considerations are introduced pertaining to a suitable extrapolation for energies of vibrational levels in approaching the dissociation region, a problem that is not regarded as solved.
The calculation range of new estimates for ideal gas thermodynamic functions for water is from 200 K to 10,000 K. It uses Q=I,Q, eKp(^Gyhc/kT), (1) with V indicating all sets of values for v,, Vj and Vj giving stable vibrational states. Gv is the vibrational energy in cm^' units. Q, is the rotational partition function for the given vibrational state. This is taken as 0.=G,V.4dQ, with the usual semi-classical relation where the symmetry number cr is equal to 2 for the H2O molecule. A,,B, and C" are the principal rotational constants and dQ, is a "low temperature" quantum correction such as that of Stripp and Kirkwood [5] as used in reference [1]. The centrifugal distortion and stretching effects are here represented by the factor/,, related to the "Wilson [6] centrifugal effect constant p," by /v = exp(/iO, where h,=p,T. The treatment of centrifugal effects for eight of the lowest vibrational levels is based on spectroscopic data interpreted with a Watson type rotational Hamiltonian. The extension to high bending quantum number is made on the basis of model estimates.
An examination is also made as to plausible magnitudes for effects of rotational cut-off in the dissociation region. Ad hoc adjustments in the approach to the dissociation energy region have been made to preserve approximate symmetry versus v, and V3 quantum numbers.
For temperatures upward of 200 K, the data and empirical representations of Camy-Peyret et al. [7] could be used for seven vibrational levels above the ground state for which detailed parameters for rotational Hamiltonians are available. These include values for the principal rotational constants Ay, By and C" and for G,, the energy for the vibrational level at zero rotation. Their reported values for G, were used also for four other vibrational levels. Their G, data for the (1,1,1) state were not used, as its partner in resonance, (0,3,1), had not apparently been similarly covered. A placement estimate for a (0,4,0) level based on a resonance shift from Benedict [8] was used to complete G, values for the resonating triad including the (2,0,0) and (1,2,0) levels. For the ground vibrational state, however, the slightly differing results of the more recent analysis by Kyro [9] were accepted in the later calculations. The overall course of vibrational energies versus vibrational quantum number was taken to follow an empirical data fit by Benedict [8], but with some adjustment in the higher energy regions to be consistent with other data, such as dissociation and heats of reaction.

The Bending Model
On the basis of spectroscopic data it can be inferred that vibrational bending produces large effects on f, and special effects on the total energy. The lack of extensive data and the impossibility of making reliable long range extrapolations of directly fitted polynomial representations of data have led to the present numerical exploration based on a simple bending model. For this, somewhat crude evaluations of WKBJ integrals have been used, based on an approximate bending potential.
The potential U is taken as the product of an empirical basic function Uo and an empirical correction function Ur,, as U=U,-U^.

(4)
The basic function Ua here involves two Lorentz type terms,  The calculations were made with bonds of fixed length "d," taking the atoms as point masses. The quantum phase integral (7) may be put in the form (8) for which conventional wave number energy units are convenient. The quantities mo and WH represent masses of the respective atoms.
In the absence of rotation, U{i) is identical with U of eq (4). For rotation purely about the principal axes there are three cases: A value of 75.73 degrees was used for <^, in adjusting between A, and Bg indications of a preliminary data fit of Benedict [8].
Energy values on effective potential curves including rotation were approximated by where 0.0001 has been added in two of the denominators to avoid accidental division by zero.
Values of V2+I/2 were computed at 17 to 40 suitably spaced values for energy "£"' or "G" up into the 100,000 cm~' region for these three cases of rotation about the principal axes, for several values of each /,^ ranging from zero to 400 in the "a " and "c" cases and somewhat further in the "b" instance. Four point Lagrangian interpolation was then used to obtain energies corresponding to integer Vi values at each of the chosen /,^ values. Effective values for the principal rotational constants were estimated according to  Estimates for the centrifugal effect contribution associated with rotation about separate principal axes were obtained in the following way. In the case of the "J?" rotation, for example, and empirical representation for B versus Jt^ was used for each of various v^ values. Numerical quadratures were performed for a partition function contribution as for various temperatures where fi=hc/kT. This determined an effective B according to 5e,r=77/(4^(e)/).

(20)
A similar procedure was used for the "A " rotation.
No appreciable contribution came from the "C" case.
The combined contributions to the Wilson centrifugal effect constant based on the rigid bender model are shown in figure 5 for values of vj ranging where t~T/1000. For these, the numerical parameters that follow are based on combined estimates as from eq (21), fitted approximately by inspection and graphical processes. With v representing vt, the constants *, and rf, were taken as 6o=2.2616(l +0.242V +O.O37v'+0.0OO83f') as an approximate representation of the calculated values. Values for ^" 5" and Q as interpreted via data of figures 2, 3, and 4 are shown in figures 6, 7, and 8 by solid circles. The solid curves are from a fit of spectroscopic data by Benedict [8] running up to V2=4, with Vi and Vj also extending up to 4. The X 's in figure 6 show individual estimates at the quantum number ^" = 1, indicating that for A, extrapolation above V2 = 7 encounters some imperfection in the traditional representation.
The dashed curves in figures 6, 7, and 8 are given respectively by with 7^=0.1473v/9.2806, where v represents V2. The expressions for A, and B, are roughly representative of the directly indicated model results as shown by the solid circles. For C" the equation here is basically a rearrangement of Benedict's equation, although the model results suggest that a different curve might be better.
The open circles show values that have been adjusted from the solid circles according to an allowance for bond stretching based on OH bond data. The effects may be summarized in part as due This causes the adjusted value for the "true" v corresponding to each model estimate to occur for a larger value than the v (efQ. The change in bond length also reduces the magnitudes of principal ro-tational constants from the model according to the square of the same ratio. An indication via the bending model as to the dependence of vibrational energy on the Vj quantum number is also given directly by the Vj+l/l values versus Eo values from the WKBJ integrals with Ji^ = 0. The results are shown graphically in figure 9 by the large open circles, obtained directly using the bending model without any bond stretching allowance. Ad hoc adjustments for bond   Stretching were then taken from vj +1/2 differences between the small open circles for free rotation without bond stretch and the small solid circles for free rotation with bond stretch obtained with bond data transferred from the OH bond of the OH diatomic molecule.
As to the curves of figure 9, the one labelled "a" is for three terms in v, when v, = V3 = 0, from a fit by Benedict with 2v,'s up to 4. Curve "c," as G = 1608.034V2-11.748v/-^ 1.643vj' + 0.0937v2^ is obtained versus V2 alone from basically the same data with a fourth term included in the fitting. Curve "b," used in the ideal gas calculations of 1979 [1], was obtained from "a" by adding the two terms Curve "d" is represented by a rational function with coefficients chosen to fit the large open circles V2 from the potential model without bond stretch.
Curve "e" is also by a rational function, but fitted only to the large solid circles at low quantum number. The last curve, "f," involves a combined locus asymptotic to a rational function curve at low Vi and to a straight line at large Vj. It is given by where fl=2300(v ^7), c = 16(X), cf = 157.4661132, and b = 1601.337v(l -0.4105368v + 0.0706926v' where v represents V2. It is this curve that represents values used in the present calculation of thermodynamic functions.
In comment on the many digits used for these constants, this is to be attributed not to any extreme accuracy, obviously, but at least in part to a regard for correlation between coefficients and to a desire to retain significance in difference type effects.
It is perhaps well to admit at this point that great accuracy is not claimed for the bending potential used. The application of the bending model is seen as quite successful, however, in providing a clear indication of the rather moderate magnitude of change in distant extrapolation as compared with the results for free rotation. Still, the uncertainty in these extrapolations must be very appreciable.

Centrifugal Data
As shown in an earlier publication [1], the effects on the rotational partition function for a given vibrational state due to centrifugal distortion and bond stretching may be obtained in semi-classical approximation from the integral g = 87r' /i -' I j Lxp{-H/kT)dP, dP,. dPj (28) with H=hciWo+n{), where fFo=2B,/',' and where H, represents the remaining part of the Watson-type rotational Hamiltonian. This uses Si = C" B2=B, and B}~A" the principal rotational constants, with P,-Px, Pi=Py and P}=P,.
The present application of the method has been carried to the evaluation of five coefficients versus temperature. In the 1979 application of the method, the evaluation was made to three coefficients, only. The hope was that the added detail would provide a better overall representation of thermal effects implied by the spectroscopic data.
There may typically be about 20 to 30 terms in current realizations of i/i, involving coefficients of various powers or products of powers of P^==P^^+P/ + P^\ '<p^i'-^p^i_pi^ and P,\ The factor e\\>{-hcHx/kT) may be expanded as a Taylor series in powers of //,. With 20 terms in H^, running generally up to the 10th power in P/s, but with a 12th power in P,, there are found to be 98 terms in Hy and 35 in i/,' in the range through the 12th power in P,. When the terms mp^ and P^y" are expanded in powers of P;,', Py'^ and /'^^ the number of separate terms to use in evaluating Gaussian integrals becomes quite large. For each, the integrand is a product over i = \, 2, and 3 of F,^"' exp{-hcB,Pi^/kT). As a result, each separate term is of the form e. = 0v" ,^n ^ F(«,) {kT/hcB,r where A computer program has been arranged for carrying out the preparation of the correction factor as a series in powers of temperature using the empirical constants of the Watson-type Hamiltonian. With terms up to the 5th power of Hy covered, the corresponding coefficients are computed to give the logarithm of the correction factor as a series h,=p,T{\+a, T^a^ T'+a, T'+a, T).

(30)
In practical application, this has not appeared to provide a well-behaved form when used for moderately high temperatures. Accordingly, the program next computes the corresponding coefficients in a Fade approximant or rational function form, The Fade form appears to be much better adapted to computation in ordinary circumstances. For some higher vibrational states, however, there can still be comphcations such as the occurrence of negative values for the coefficients Cj and di which are for the highest powers of T. Happenings of this type appear to be somewhat dependent on the source of the empirical Hamiltonian constants used.

K=p,n\+c,T+c^Ty{\+d,T+d^T\
The program is interactive in asking for values for the principal rotational constants, A, B, and C, and then for the highest power of temperature to be covered (up to 6 as set up). It then asks progressively, in a selected order for values for 30 coefficients in the Watson Hamiltonian, identifying each by a coefficient name in "string-variable" form, e. g., DELJ, etc. These happen to be in the order shown schematically by:  where J**2 is JX**2+JY**2 + JZ**2 and JXY**2 is JX**2-JY**2.
Values for the Pade constants as obtained from the available Watson Hamiltonian constants for the eight observed vibrational levels, based largely on the work of Camy-Peyret and Flaud [7] are given in table 1.
A listing of the program is included in the appendix. Further discussion of the results will be reserved for a later section dealing with table comparisons.

Rotation at High Temperature
While the Pade form for centrifugal effect seems better adapted for calculation than the simple power series form which encounters convergence problems of an erratically varying sign type, there are other consideration if extrapolation to very high temperature is required.
The rotational quantum numbers can increase up to some limiting large values as "/"," "/*," or "Z^," beyond which centrifugal force would cause the molecule to break apart by bond rupture. The limiting rotational energy for rotation about any principal axis would be of a general magnitude indicated by where 5, refers to ^v, B, or C", according to the axis involved. However, the affected moment of inertia at bond rupture would be appreciably increased by bond stretching over its ordinary value. A semi-classical OH bond model study on the principal rotational constants at zero vibration suggests that the ratio "r" between principal constants at maximum versus at low rotational quantum number should be about 0.195 as r^ for A, 0.260 as r^ for 5, and 0.310 as re for C.   where R, the factor of centrifugal increase, is estimated as

R=[l+iRo^l)iD,/Don
with Ro=i.O and possibly 5 = 1 may be a useful speculation as to plausible behavior. Here D, is the additional energy Do-G,, to reach dissociation for the vibrational state without rotation. As to acceptable values for DR for higher vibrational states, it appears useful, with DM as dissociation energy including rotation, to note that by a logarithmic plot of y={DM-D)/D versus x={D-G,)/D for the OH diatomic potential, an approximate representation is y={l/9)x'^'*. Estimates of a similar magnitude can also be found from the expression;; ={\/A)x, which is a form easier to use. The latter choice provides an approximate relation Dn^D(x+y), or D=D-G, + (D-G,)/4=1.25(£>-G,).
As in the discussion leading to eq (34), an estimate for the rotational partition function may be based on an integral using an ellipsoidal shell with semi-axes n{E/Ay'\ M,E/BY^ and fi{E/CY\ The "volume" of the shell between energies hcE and hc(E+dE) is In fi' {ABCY^ E"^ dE. In the evaluation of density of states as measured by II (/»"' dpidqi), there is a factor 47r for orientation of the total momentum vector and Iv for position of the rotator in making one revolution. Thus the number of states available within the energy shell needs the factor STT^A"' to be included, giving dN = l{ABC)''^E'"-dE.
The last factor is /, or exp(pr) for centrifugal effects according to eq (2).
If the integration is extended only to a rotational energy E=D, or for x to x=hcD/kT, the result for each term involves an incomplete gamma function. QS=l{kT/hcy'\ABCy\y{yi^,) The recurrence relation y (a +1 ,x)=a y (a ,x) -x''exp(-x) is used to relate all later terms to the first one. Values can be found for 7 ( with Q°==^n'^KkT/hcy'\A^,Cy)-"^ cr~\ as for a rigid rotator, tr being the symmetry number, with 0, = 27r-'^^y (3/2,x,) and with/, representing exp(p7') as in eq (2)  As a variant study based on the derivation leading to eq (38), one may remove the Boltzmann factor exp{-hcE/kT) and consider the integration up to an energy E.
This was tried on the ground state and some others. Conversion of the polynomial to a Fade-Wilson exponential form quieted a term-wise sign fluctuation effect. However, in extension to very large E, a condition of excessively large computed Q"^ was encountered. This was due to the exponential factor becoming grossly over-sized. For the ground state, as E rises from \.E + 5 to l.£+6 cm"', the computed centrifugal factor rises from 6.5 to over l.E + 1. This result is contrary to the previous estimate of a limit for the centrifugal factor of the order of /? = 8 or less, as in eq (36). The catastrophe can obviously be avoided by enlarging the denominator by including a term fi?3 E\ with the ratio bl/di near to 2 or to In/?. Logically, the parameters should be chosen again in such a way that the expansion into a power series would remain unchanged through the first five coefficients. If s represents b2/d1, and with B\, B2,D\ and D2 to represent original values of b 1, b2, d\ and d2, respectively, the revised coefficients can be obtained from s 1 = * 1 * * 3 2*bQ*b\''b2 \ b0**2*b2*dl -bO*bl**2*d\+bO**2*bl*d2 S2=b0*b\*dl*d2-bl**2*b2-b0*b2*dl"2 -bl*b2*d\+2*hQ*b2*d2-b2**2 -b0**2*d2**2
The same type of control adjustment should apparently be applicable to the Fade form in terms of temperature in a normal computation. However, in actual application to a multitude of levels, there could seem to be a possibility that the final Fade constants might not always be positive, due to numerical accident. A requirement that </3 and ^2 be positive can be met by using the absolute value of This might preserve only four instead of five coefficients of the series leading to a Fade development and no absolute guarantee is known to exist against occurrence of a zero denominator.
Another simple scheme for keeping the rotational Q below the limiting QM value for a given vibrational level has been patterned after the familiar relation of a hyperbola to its asymptotes represented as a combination of loci. With QM as an excessive rotational partition function without cutoff, an estimate with cut-off included might be

(46)
The quantity q is to be taken in a convenient form showing an acceptable temperature dependence. Results of a graphical study for rotations about principal axes, using bond stretching of the OH molecule, lead to a provisional suggestion that a usable form might be q=q\Z exp(-z) where z-q2hcD/kT, with q\=OA and qi=\.%. Other representation schemes may reasonably be more suitable, however.
It is conceded that direct rotational cut-off effects are fairly small even for temperatures at the top of the range of the present tabulation. However, an indirect effect in the extrapolation is not quite so negligible. In making evaluations based on the empirical constants of Benedict, it was found that gross differences in behavior between v 1 and v3 dependences were produced with vl or v3 large, particularly as v2 was increased so as to be more than a small integer. This characteristic is attributed to the effect of the long range of the extrapolation with equations fitted to data at low quantum numbers only. A more uniform behavior has been obtained by a revised procedure for treating the empirical vibrational energy.

Energies for High Vibrational Levels
In the last several years the method for estimation of high vibrational levels appears to be changing, involving such new developments as are referred to as localized bond excitation and local mode description [10]. Special potential forms can be used for such calculations with constants converted [11] from empirical values found with a conventional valence bond system formulation and normal coordinate analysis [12]. Potential improvement based on direct comparison between computed and "observed" levels could be an ultimate objective. A hazard at the outset in this approach may be a sensitivity to the correctness of identification or assignment of spectroscopic data on which at least the original numerical constants are based. As to direct a priori quantum mechanical calculation of levels for the molecule as a collection of nuclei and electrons, based on general physical constants, it appears that significant advances have been made on this intrinsically difficult endeavor.
Somewhat approximate agreement with vibrational fundamentals has been obtained [13,14] but whether a similar quality of prediction could be achieved for higher vibrational levels may be in an area of pure speculation. A priori calculation appears informative in regard to excited electronic states [15], in an energy domain beyond the range of the present treatment.
Even if a local mode description will prove ultimately more reliable than the conventional approach, it has appeared expedient to continue for the time being with the older formulation, for which the necessary parameters are at hand. It appears plausible that newly and correctly calculated levels should on the average agree tolerably well with the old values of corresponding description. This is thought to be the usual situation for a group of "interacting levels" in a so-called resonance situation.
The vibrational constants used here are based on a formulation by the late Prof W. S. Benedict Journal of Research of the National Bureau of Standards for levels with resonance shifts removed as indicated earlier.
The present proposed innovation in regard to vibrational energy is to suppress the long range effects of Benedict's fitting on the basis that the fine details of fit while relevant in the region of fit in the low quantum number range (v/<5) may still not be numerically reliable when extrapolated to large vi. For the various "small" quadratic and cubic terms, involving "v products" (=p), an extrapolation by replacement of "p" by p/[\+(p/a)'^ has been used with k=6 and with the parameter "a" chosen differently for terms quadratic and cubic in the v's. (65 versus 460) This causes these terms to become small in the approach to the dissociation region.
For the estimation of vibrational levels in the region of large v 1 and v3, the procedure adopted was to take the energy as given primarily by a quadratic jointly in vl and v3, much as in the case with a Morse potential in a diatomic molecule. Thus, in the case with v2=0, the form for this main part of the vibrational energy becomes GL(vl,0,v3)=i>l vl + »)3v3 -;cllvl(vl-l) ^;c33v3(v3-l)-xl3vlv3. (49) The anharmonicities for this were chosen so as to agree with energies of dissociative reactions based on thermochemical data. For vl or v3 increasing singly with the other at zero, there is dissociation according to H20=0+2H at about 76721 cm"'=£>. For vl and v3 equal and advancing together, dissociation is taken to be according to H20 = H-KOH at about 41280 cm~'=Dm.
In the cases of vl and v3 advancing singly, the energy is given as in give X 13 = 248.95779 kaysers but the partially suppressed residual contributions of other constants at this dissociation energy (near vl = v3 -11) raise xl3 to 251.3489 cm-'. A multiplier factor [1^0.0028 v2-0.00013 v2 (v2-1)] has been introduced for i), vl+i>3 v3 to allow for a diminishing energy increment to dissociation as v2 advances upward above v2=0. All such adjustments are compensated for in the low quantum number range so as to preserve the behavior there according to the empirical data fit of Benedict [8], Thermodynamic Tables for H2O Two sets of tabular values have been included as prospective thermodynamic quantities for the ideal gas state of the light isotopic water molecule. These are here designated by their dates of computation, which were 1982 and 1984.
For the 1982 table, appearing here as table 2, parameter values used were influenced by results of computations for a rigid bender model, adjusted further for bond length increase by centrifugal stretching due to a rotational character of motion in the bending vibration. These included indications as to the v2 dependence of the principal rotational constants, the extrapolation of vibrational energy to high v2 values, and the course of the ordinary centrifugal effects to high v2 and elevated temperatures, using a five parameter Fade formulation.
The 1984 table, shown here as table 3, includes  the innovations of the 1982 table, and a few others, also. In the approach to dissociation at high v 1 and v3, the behavior of Gv was taken as essentially quadratic in vl and v3, in resemblance to the known diatomic behavior with a Morse potential. Special functions were used to fade out the detailed higher order terms arising out of Benedict's Gv fit at low vibrational quantum numbers. For eight low-lying vibrational states, numerical values were inserted via the computer program for observed  The final column gives the number of vibrational levels involved in the final state sum. The next to the last column gives the fractional reduction in the state sum due to rotational cut-off, according to the "locus-asymptote" estimate used. values for vibrational energy, principal rotational constants, and their five member Pade centrifugal parameters, based on reported spectroscopic data analyses using the Watson Hamiltonian formulation. A rotational cut-off approximation of a "locus-asymptote" type was also introduced, but with little apparent effect up to 10000 K.
It is natural to see the difference in values between the two tables as relevant to their uncertainty. It is presumed that the disagreement in values should be attributed to effects in changes in level distribution, which may reflect the ad hoc modification of level description for the later table.
It had appeared reasonable to maintain a favorable view of progress in raising the number of constants based on the Watson Hamiltonian data from three to five. However, it is now recognized in retrospect that some basis for reserve exists. As used, the program for finding Pade constants was able to produce the five parameters as desired even when the Hamiltonian parameters were not complete to a corresponding extent. This might be termed a "spill-over" effect akin to the forming of product terms in a series development. The highest power of r directly included as a contribution to "J times the Wilson constant" may be obtained by taking the highest net power of J's in the Hamiltonian, dividing by 2, and subtracting 1. On this basis, the ground state and first excited vibrational state, (000) and (010), may be "complete" through the 5th power. The states (020) and (030) show fitting in the 4th power, and the states (100), (001), (110) and (Oil) include only into the 3rd power. One may hope that a moving of the Pade process into the Hamiltonian will lead to a more uniform treatment [16].
As comment on our present use of a "preliminary" 1972 data formulation received from Professor Benedict [8], we accepted his view that his was better than that of Khachkuruzov [17], of 1959. We note that a more recent vibrational energy formulation presumably of comparable quality was published in 1983 by Bykov, Makushkin and Ulenikov [18], and could in all probability provide a similar basis for a table of thermodynamic quantities.
It appears that greater consideration should be given to recent work such as that by Child and Lawton [19] on local mode representations of vibrational states. However, at this time it is not clear how energies for the entire manifold of vibrational states would be reliably and conveniently given for the calculation of thermodynamic functions on such a basis.

Conclusion
The objective in this study has been to obtain an improved extrapolation of the ideal gas table to higher temperatures. The procedure has made use of direct data, augmented with numerical estimates based on simple physical models. It is hoped that this might provide a realistic approach to better sum of state estimation.
Although the models have involved some numerical choices that were not at all rigorous, the results may allow such comparisons as may lead to an informed appreciation of the problems remaining for the reduction of uncertainties.