Crib-Sheet to SPFIT and SPCAT of Herbert Pickett
  

        The SPFIT/SPCAT program suite of Pickett (available from JPL Molecular Spectroscopy) is currently the most widely used tool for analysis of high-resolution spectra and for their prediction, most notably for spectroscopy databases for atmospheric and astrophysical applications.  The appeal of the package lies in flexible construction of the molecular Hamiltonian and in very efficient internal matrix factorization.  This allows a single package to deal with widely varying molecular problems and for very fast operation.  The package is well known for its steep learning curve so that various help documents (such as this one) are available.

        This Crib-Sheet contains loose notes based only on applications of Pickett's SPFIT/SPCAT (or CALFIT/CALCAT) program package tried by kisiel@ifpan.edu.pl and dating back to the 1990s. These notes should only be used once you are familiar with the documentation of the package:

        Several versions of the programs have been published by HMP since these notes started. Some comments may therefore have been superseded by developments, and may not apply to the latest program versions. The doublet PI subroutine set has not yet been used and statistical weights on intensities have generally been disregarded.  Many years have passed since these note were started and I currently find that I most often tend to jump directly to the Parameter indices and names section.  I also recommend the SPFIT Notes and Hints and SPCAT Notes and Hints, since these cover points that either pose difficulties or are not obvious to the users.

updated: 26.I.2022

        Notes:

     

      

    

SPFIT = the fitting program

             
             
input:   .par   -   control flags and declarations + values of constants
    .lin   -   transition listing
output:          .bak   -   backup of the original .par file created at the beginning of each fitting run
    .par   -   modified .par file with constants resulting from the fit
    .var   -   as .par but with errors and Cholesky decomposition of correlation matrix for use in SPCAT
    .fit   -   main output file
    .bin   -   binary output file containing parameters and variances for use by the predictive program (obsolete)

       SPFIT in the 1994 version worked only for J up to 100, since the March1997 version the J limits are (-270,369) for SPCAT, (-99,999) for SPFIT

     

       Linear top

        Third line in the .par  file (single quadrupolar nucleus, I=3/2):

       s -4 1 0 ,0 , , , , , , ,

        (three quadrupolar nuclei, I1=3/2, I2=1, I3=1, and F1,F2,F coupling):

       s -334 1 0 ,0 , , , , , , ,

        (three quadrupolar nuclei, I1=3/2, I2=1, I3=1, and I,I1,F coupling):

       s -334 1 0 ,0 , , -1 , , , , ,

       Remember to set KNMAX explicitly to zero. Also in case of symmetry equivalent nuclei use WTPL and WTML. For example for A1 and A2 states of 14N2...H2O the third lines of the  .par  file are:

       s -33 1   0 , 0 , , -1 , 1. , 0. ,,,,,,,,,,,    <---- A1
s -33 1 0 , 0 , , -1 , 0. , 1. ,,,,,,,,,,, <---- A2

       Important constants:

        100   B   or
      20000   B
-30000 C (set to follow the value of B)

       Note that other similar fitted constants (i.e. those that have identical values by symmetry) all have to be declared in this way. If this is not done the resulting constant will be equal to the conventional value multiplied by its degeneracy. For spin-rotation, for example, use

   10020000  -M.bb
-10030000 -M.cc

       In case of coupled constants the deviation of the second constant gives the overall deviation for both constants in the fit i.e. the deviation and value of constant 100 is the same as that of -30000 if the pair (20000, -30000) is used; since the Dec1996 version only the combined deviation is printed for both constants.

        Line in the .lin  file (single quadrupolar nucleus):

     J' F' J'' F'' 0 0 0 0 0 0 0 0  freq  error  weight
      

       Symmetric top, with fitting of h3 and effects of near-sphericity

       Third line in the .par  file(-1 defines symmetric top, 1 is for prolate symmetry, change to -1 for oblate top):

       s -1 1 0   ,,,,,,,,,,,,,,,,,,

       Important constants (prolate):

      10000   A
20000 B
-30000 C (set to follow the value of B)
60000 h3

       Alternatively:

       1000  A-B      :prolate
3000 B-C :oblate
100 B

        Line in the .lin  file:

     J' K' J'' K'' 0 0 0 0 0 0 0 0  freq  error  weight

        Note that K=3 (K=6 etc.) doublets are to be specified as

     J+1 -3 J  3   0 0 0 0 0 0 0 0
J+1 3 J -3 0 0 0 0 0 0 0 0 etc.

       The sign of K alternates between odd and even values of J. For a positive value of h3 and even J it is the K=3 ← 3 component in SPFIT that is the higher frequency component, and 3 ← 3 is the higher frequency component for odd J. Also note that due to this cross type of declaration it is necessary to explicitly declare both components of unresolved K=3, K=6 doublets earlier than in other programs.

     

      Doubly degenerate v=1 excited state in a C3v symmetric top

       Treatment of l-doubled v=1 states became possible in the March1997 version of the program. The declaration for this problem has evolved with SPFIT and is currently (2008) somewhat different than originally used.  The key change is that many different tops are recognised and a C3v top is now declared by using IAX=6.
       A doubly degenerate state is specified as two states by using the flag EWT1 in thethird line of the .par file: EWT1=1 is l=+1, EWT1=2 is l=1. Singly degenerate states are specified with EWT1=0 for l=0. The EWT1=1,2 states should be declared on adjacent cards. It is also useful to include the ground state in the declaration and the example below is for such a case. The third line in the .par  file is to have two continuation lines (the settings of WTPL, WTMN and EWT0 ensure correct relative intensities for a case with off-axis protons):

       Prolate top:

  s -1, 3, , , , 6,  2,  2,   -1,   1, 1,      ; this is the g.s.
, 1, , , , 6, 2, 2, -1, 101, , ; this is l=+1
, 2, , , , 6, 2, 2, , 201, , ; this is l=-1

       Oblate top:

  s -1,-3, , , , 6,  2,  2,   -1,   1, 1,
, 1, , , , 6, 2, 2, -1, 101, ,
, 2, , , , 6, 2, 2, , 200, ,

       Note that the diagonalization option DIAG (twelfth parameter on the first of these lines) is to be set to 1 i.e. for a full projection assignment. Of the pair of doubled states the diagonal constants are to be declared only for the state EWT1=1 and the relevant constants are then:

           prolate                oblate
            1011  A-B               3011  B-C
111 B 100 B
or:
10011 A 10011 A
20011 B -20011 B
-30011 C 30011 C
          200011 -2 A zeta.EE     600011 -2 C zeta.EE
200111 eta.J 600111 eta.J
201011 eta.K 601011 eta.K
200211 eta.JJ 600211 eta.JJ
201111 eta.JK 601111 eta.JK
202011 eta.KK 602011 eta.KK
40012 q/2 40012 q/2
40112 -q.j/2 40112 -q.j/2
or:
210012 q 610012 q
210112 -q.J 610112 -q.J

       for off-diagonal elements in the form 0.25[q-qJ J(J+1)..]. The second definition for q involves much more diagonalization and is therefore not recommended.

       Transitions which are normally degenerate are best specified by using positive K. Those for l=+1 are declared as state EWT1=1, while those for l=1 as EWT1=2. The K=0 (Kl1=1) is EWT1=2. The two K=l1 doublets are both EWT1=1 and their K assignment alternates with parity of J:

        Prolate (positive q and ) and EVEN J'':

 J' -1  1  J''  1  1 .......     A+ i.e. upper freq. component
J' 1 1 J'' -1 1 ....... A- i.e. lower freq. component

        Prolate (positive q and ) and ODD J'':

 J'  1  1  J'' -1  1 .......     A+ i.e. upper freq. component
J' -1 1 J'' 1 1 ....... A- i.e. lower freq. component

        Oblate (positive q, negative ):

 J' -1  1  J'' -1  1 .......     A+
J' 1 1 J'' 1 1 ....... A-

        Oblate (positive q, positive ):

 J' -1  1  J'' -1  1 .......     A-
J' 1 1 J'' 1 1 ....... A+

       If the E-symmetry states declared as v=1 and v=2 undergo Coriolis coupling with the A-symmetry state declared as v=0 then connecting constants are:

           prolate                oblate
          600001                  200001  w = 2 root2 Omega B zeta.EA
600101 200101 w.J
601001 201001 w.K
600201 200201 w.JJ
601101 201101 w.JK
602001 202001 w.KK etc.

where Ω = [sqrt(ωa/ωe) + sqrt(ωe/ωa)]/2 and ζEA is the Coriolis constant connecting the E and A states.  The oblate case has been extensively tested in Kisiel et al.  J.Mol.Spectrosc. 251(2008)235-240, which contains a more extensive table of parameter identifiers and their translation into Hamiltonian terms.

     

       Straightforward asymmetric top fit

       Third line in the .par  file:

       Reduction-A, prolate (representation Ir):

       a 1  1 0   ,,,,,,,,,,,,,,,,,,

       Reduction-A, oblate (representation IIIl):

       a 1 -1 0   ,,,,,,,,,,,,,,,,,,

       Reduction-S, prolate (representation Ir):

       s 1  1 0   ,,,,,,,,,,,,,,,,,,

       Reduction-S, oblate (representation IIIl):

       s 1 -1 0   ,,,,,,,,,,,,,,,,,,

       Standard constants are as given in APPENDIX I. The difference in the third .par line between the use of A- and S-reductions is only in declaring different .nam files for constant names (and these files are currently optional). It is also possible to fit linear combinations of constants such as prolate (B+C), (BC) declared using:

  20000  value   /(B+C)/2
-30000 value
40000 value1 /(B-C)/4

or, in an alternative way, using:

  20000  value    /(B+C)/2
-30000 value
20000 value1 /(B-C)/2
-30000 -value1

which is analogous to the way nuclear quadrupole coupling constants are declared in the 'simple examples' at the end of SPINV.PDF.  This type of declaration is possible since values from successive instances of a parameter with the same index are additive.
In a similar way it is possible to deal with values of rotational constants in different vibrational states, for example:

  10000  value    /(A0+A1)/2
-10011 value
10000 value1 /(A0-A1)/2
-10011 -value1

       Line in the .lin  file:

  J' K-1' K+1' J'' K-1'' K+1'' 0 0 0 0 0 0   freq  error  weight
      

       Coriolis coupling between two asymmetric rotor vibrational levels

       The lower of the two levels is to be given the identifier v=0, the upper v=1. The constants are then:

    10000 A rotational constant for state v=0 etc.
10011 A rotational constant for state v=1 etc.
11 Delta E i.e. difference in vibrational energies of
the two states
200001 i G.a P.a
210001 F.a (P.b P.c + P.c P.b)
400001 i G.b P.b
410001 F.b (P.c P.a + P.a P.c)
600001 i G.c P.c
610001 F.c (P.a P.b + P.b P.a)

        Third line in the .par  file (A-type, prolate):

 a 1 2 0   ,,,,,,,,,,,,,,,,,,

        Line in the .lin  file:

 J' K-1' K+1' v' J'' K-1'' K+1'' v'' 0 0 0 0   freq  error  weight
      

       Vibration-rotation transition between two vibrational states

       Specification is identical to the case above with the exception that only v''=0 to v'=1 transitions are specified (for absorption) and of the coupling terms only the energy difference term:

       11 Delta E

is used. Transition frequency and its error can be specified in cm-1 if the error is made negative.

     

       Coriolis coupling between two excited vibrational states with use of vibration-rotation and rotational spectra

        Three states need to be specified, ground state v=0, and the two excited states v=1 and v=2. The third line in the .par file is then

     a 1 3 0   ,,,,,,,,,,,,,,,,,,

        and some pertinent constants are:

      20000   B for the ground state etc.
20011 B for first excited state etc.
11 vibrational energy E for first excited state
20022 B for the second excited state etc.
22 vibrational energy E for the second excited state
600012 G.c coupling constant
610012 F.c coupling constant
200012 G.a coupling constant
      

       Asymmetric rotor with a single quadrupolar nucleus

       Values of the spin quantum number have to be rounded up to the nearest integer, so that the odd half integral spins are declared as: 1/2→1, 3/2→2, etc. The independent quadrupolar constants are specified as:

 110010000  1.5 Chi.aa             :Prolate
110040000 0.25 (Chi.bb-Chi.cc)
 110030000  1.5 Chi.cc             :Oblate
110040000 0.25 (Chi.bb-Chi.aa)
 110610000  Chi.ab                 :Prolate and oblate
110210000 Chi.bc
110410000 Chi.ac

        Third line in the .par  file (A-type, prolate, spin 3/2):

 a 4  1 0   ,,,,,,,,,,,,,,,,,,

        Third line in the .par  file (S-type, oblate, spin 3/2):

 s 4 -1 0   ,,,,,,,,,,,,,,,,,,

        Line in the .lin  file:

 J' K-1' K+1' F' J'' K-1'' K+1'' F'' 0 0 0 0   freq  error  weight
      

       Asymmetric rotor with two quadrupolar nuclei

       Notes from the single nucleus case above are applicable. Prolate rotor constants:

 110010000  1.5 Chi.aa                 \
110040000 0.25 (Chi.bb-Chi.cc) |
110610000 Chi.ab | nucleus 1
110210000 Chi.bc |
110410000 Chi.ac /
220010000 1.5 Chi.aa \
220040000 0.25 (Chi.bb-Chi.cc) |
220610000 Chi.ab | nucleus 2
220210000 Chi.bc |
220410000 Chi.ac /

        Third line in the .par  file (A-type, prolate, spin1=3/2, spin2=1, and I,F coupling):

 a 34 1 0   , , , -1 ,,,,,,,,,,,,,,,

        Line in the .lin  file (for I,F coupling):

 J' K-1' K+1' I' F' J'' K-1'' K+1'' I'' F'' 0 0   freq  error  weight
      

       Spin-spin

        The D coefficients are for Hss = I1.D.I2 and are (prolate):

 120010000  1.5 D.aa
120040000 0.25 (D.bb-D.cc)
120610000 D.ab etc.

        The first two digits define interacting nuclei. Note that the constant c3 for diatomics in Landolt-Bornstein vol.II/6,table 2.9.1.2 is Dbb=(1/2) Daa, as only Dbb is spectroscopically determinable for linear and diatomic molecules.

     

       Spin-rotation

        The M coefficients are for Hsr = + I.M.J and thus:

  10010000  M.aa for nucleus 1
10020000 M.bb for nucleus 1
10030000 M.cc for nucleus 1
20010000 M.aa for nucleus 2 etc.

        For linear molecules the use of:

  10000000  M for nucleus 1

        is equivalent to the declaration:

  10020000  M.bb for nucleus 1
-10030000 M.cc for nucleus 1 set equal to that of M.bb

       Note that, following Flygare, a reverse Hamiltonian definition is often used i.e. Hsr = −I.M.J and the resulting M's have reversed signs relative to those above.

  

        SPFIT parameter errors

        There are two relatively subtle points associated with understanding the errors in fitted parameters that are printed by SPFIT when you use the default FRAC=1 parameter (the seventh parameter in the second line of the .par file). These points are primarily of relevance for small datasets, and diminish in importance for large datasets fitted to close to RMS deviation of 1. Nevertheless if one wants to obtain exact correspondence between the results of SPFIT and those from other fitting programs then some corrections may be necessary.

        Since weighted least-squares fitting programs all use the same basic statistical treatment it is useful to refer to the seminal paper on the subject (referred to below as ASZ):

D.L.Albritton, A.L.Schmeltekopf, R.N.Zare, "An introduction to the Least-Squares Fitting of Spectroscopic Data", in Molecular Spectroscopy: Modern Research, Ed. K.N.Rao, Academic Press, 1976, Chapter.1, p.1-67.

POINT I:

        The two RMS quantities in the .fit output of SPFIT are:

       (MICROWAVE RMS) = sqrt ( Sum[(obs-calc)**2]/Nlines )
(RMS ERROR) = sqrt ( Sum[((obs-calc)/error)**2]/Nlines )

        Thus the deviations are worked out on the basis of a straightforward number of lines in the dataset, Nlines, and not the degrees of freedom Ndegf=Nlines-Nconst, where Nconst is the number of fitted constants. In ASZ, say Eq.31, the degrees of freedom, (n-m), are used. As a consequence, the deviations evaluated for small datasets by using Nlines instead of Ndegf may differ appreciably. Of course, as HMP points out, one can readily construct data sets with pathologies that make the assumption Ndegf=Nlines-Nconst inadequate, but it is nonetheless an assumption that is generally made. Even more suspicious is the statistics of relatively small numbers when Ndegf/Nlines is significantly less than unity, but this has never stopped spectroscopists from working in that regime!

POINT II:

        In the standard application of weighted least-squares (section E5 of ASZ) only relative values of assumed experimental errors are of importance, and the estimated errors on parameters of fit (square roots of diagonal elements of the variance-covariance matrix Eq.33 ASZ) are insensitive to multiplication of assumed measurement errors by a constant factor c. If statistics are sufficient then the errors are standard errors (i.e. for 67% confidence level). For a fit with equal weights (same assumed error for all frequencies) you should get the same result irrespective of the value of the assumed error. This is not the case for SPFIT in which:

       "The printed error in fitted constants is at the level of confidence assumed in declaring the uncertainties in experimental measurements in the .lin file".

        In this case the absolute, and not only relative values of such measurement errors are of importance, and the printed errors cannot be automatically assumed to be standard errors. The choice made in SPFIT has been to forego the σ2 term in Eq.32 of ASZ. This has the consequence that even for an equally weighted fit changes in the values of the assumed measurement error will lead to changes in evaluated parameter errors.

        The two points above are illustrated by fits made for D2S, where the results from the original publication (as reproduced in Table 8.24 of the 3rd Edition of the Gordy&Cook textbook) are used as the basis for the comparison. In the original work the cited errors are at the 95% confidence level, and were obtained by straightforward doubling of the derived standard errors. Those numbers can be seen to be in good correspondence with doubled errors from the ASFIT fit, which returns standard errors. The comparison is not exactly 100% only due to different rounding at the last digit and the possibility of typos (such as in the published error on C: 71 instead of possibly 61). On the other hand you can see that the SPFIT errors are directly dependent on the assumed measurement error even for (as in this case) equally weighted lines. The primary data files are collected in this archive. A similar exercise with an even more exact correspondence with the published data is possible for T2O, for which Ndegf/Nlines is even more pathological.

        What you can see is that with SPFIT the values of assumed measurement errors are more important than with many other fitting programs, in which only relative magnitudes of such errors are relevant. In addition, once RMS ERROR deviates significantly from unity, you have to think about what is the confidence interval that the errors correspond to. Fortunately, if you were to care about actual standard errors then it is very easy to apply a suitable correction, see below.

What you can do:

1/  If you want to cite standard (1*sigma = 67% confidence level) errors on parameters fitted as a result of running SPFIT then you need to multiply the printed errors in parameters by

          (RMS ERROR)*SQRT( N.lines/ (N.lines-N.fittedconsts) )

This can be done automatically with program PIFORM, which reformats the .fit output and, among other things, converts SPFIT errors to standard errors.

2/  An equivalent procedure is to multiply the assumed errors in frequencies by the same factor as above and to refit, whence you will obtain RMS ERROR of

          SQRT( (N.lines-N.fittedconsts)/N.lines )

which will be equivalent to RMS error of unity if Ndegf was used in the statistics. The printed errors on fitted parameters will then be standard (67%) errors. If the dataset is large and the fitting model adequate then this procedure is in fact the way to establish the experimental error of the spectrometer used to make the measurements, and the RMS error will approach unity. For the D2S example with 0.2 MHz errors on frequencies, this procedure will correspond to resetting the measurement errors to:

          0.2 * 0.46794 sqrt(66/44) = 1.14602

resulting in SPFIT RMS ERROR of

          sqrt(44/66) = 0.81650

3/  For large datasets you may just set the FRAC parameter (the seventh parameter in the second line of the .par file) to the RMS ERROR from a converged fit and then refit. The resulting parameter errors will only differ from standard errors by the factor SQRT(Ndegf/Nlines). You can of course also account for this factor in FRAC and for the D2S example with 0.1 MHz errors on frequencies you will obtain actual standard errors by setting

          FRAC = 0.93588 / sqrt(44/66) = 1.14602.

This procedure, as well as point (2), allow you to obtain standard errors on transition frequencies calculated by SPCAT, as printed in the second column of the .cat file. On the other hand, you have to take care, since it might be rather easy to forget that FRAC was set to some nonstandard value that was only valid for a specific fit.

4/  As of July 2007 you can use negative values of the FRAC parameter in SPFIT, in which case parameter errors will be equal to the multiple of standard error that corresponds to the magnitude of FRAC. Hence FRAC=-1 will rescale SPFIT errors internally using the relationship in point (1) above, so that you will now get standard errors in the .fit output. Similarly FRAC=-2 will produce 2sigma (95% errors) and FRAC=-3 will lead to 3sigma (99%) errors. Furthermore, the use of FRAC=-1, for example, ensures that the errors in frequencies calculated by SPCAT will also be standard errors etc.

WARNING: It has been discovered (March2016) that the negative FRAC mechanism, such as FRAC=-1, is not without some issues that are relevant to small datasets.  In fact, the values of NDFREE in the multiplier
    -FRAC*RMS*SQRT(NLINE/NDFREE)
used to correct the parameter errors (see the documentation) are somewhat specific:

  1. when ERPAR>1.E+015 on all fitted parameters then NDFREE=NLINE-(the number of free parameters)+1
  2. when ERPAR<=1.E+015 on all fitted parameters then NDFREE=NLINE
  3. intermediate values of NDFREE between those arising in the two cases above will occur if some, but not all, ERPAR are greater than 1.E+015

- - -

The bottom line is that whatever you do: "Do not automatically assume that the parameter errors printed by SPFIT are standard errors, especially when your RMS > 1, since you are then in danger of seriously underestimating your actual errors".

      

        SPFIT Notes and Hints

      
        Be warned: SPFIT poses multiple traps for unwary users!  I have seen some problems crop up regularly over the years.  A determined scientific detective would find their unintended consequences in a surprisingly large number of published scientific papers. There is now a new danger resulting from the use of SPFIT in the "under the hood" mode in spectral analysis packages.  Some issues may, therefore, escape human inspection altogether.  You are  strongly advised to check against the various known traps:

- - -

TRAP 1 results from the specific meaning of the parameter errors obtained from the fit and the associated subtleties are discussed in the section immediately above

- - -

TRAP 2 is due to the premature assumption that a fit has converged.  Note that a given fit can be considered to be properly converged only if two conditions are satisfied:

1/ The current and the predicted values of the RMS ERROR are identical and you get screen output terminating in something like: 
END OF ITERATION  1 OLD, NEW RMS ERROR=  1.19353   1.19353

2/ The THRESH parameter (the fifth parameter in the second line of the .par file) is zero.  When the fit is properly converged SPFIT will normally reset THRESH to 0.0000E+000.  If this is not the case then you can try zeroing this parameter by hand.
   
THRESH  > 1 is especially dangerous since it indicates that the data set does not contain sufficient information to unambiguously fit the desired parameters and SPFIT has kicked in the singular value decomposition procedure in order to arrive at some plausible, but not necessarily unique solution.
   
THRESH  < 1 indicates that the search range for new parameters has been artificially constrained. The parameter errors obtained with any non-zero value of THRESH should not be used.

        A confident sign of final convergence is that the fit terminates in one iteration and THRESH=0.0000E+000, irrespective of the number of iterations declared in the parameter NITR in the second line of the .par file.

- - -

TRAP 3 is posed by the use of a priori parameter uncertainties, ERPAR, as declared next to the parameter values in the .par file.  In some situations, like cases with strong intercorrelations between several constants, it is useful to reduce the search range for selected parameters with suitably chosen values of  ERPAR. Nevertheless, in the final fit such errors should be reset to be at least an order of magnitude greater than the resulting errors in the fitted constants.  Otherwise the solution may be artificially constrained and appear more precise than it really is.

- - -

TRAP 4 results from the fact that almost all quartic centrifugal distortion constants in the .par file are negative relative to the usual convention, as you can check in this table. The only exception are the two off-diagonal S-reduction constants, d1 and d2. The recommended solution is to write the parameter descriptions explicitly in the .par file:

     200 -5.007436030889545E-004 1.00000000E+000 /-DJ
    1100 -3.854253318510433E-003 1.00000000E+000 /-DJK
    2000 -9.473858803517874E-003 1.00000000E+000 /-DK
   40100 -9.570171396103077E-005 1.00000000E+000 /d1
   50000 -1.774207632318250E-005 1.00000000E+000 /d2

You can then optionally use the program PIFORM to automatically recover the values of fitted parameters and of their errors in the form normally used for publication.

- - -

TRAP 5 is associated with impoper use of the DIAG parameter (twelfth parameter in third line of .par). The use of a suitable value of DIAG is crucial when there is a high density of states and complex quantisation.  There is, regrettably, no universal solution and for this reason there are many alternative settings of DIAG.  Conversely, the use of an inapropriate DIAG , especially in simple cases, may lead to spurious quantisation and spurious (sometimes very strong) transition series.   It is usually best to commence a new problem with DIAG=0.




        If the fit goes haywire make sure to recover the original .par by copying from the .bak version automatically saved by SPFIT at the beginning of this fitting run, using the command: copy xxx.bak xxx.par. It is also good practice to keep good .par and .fit files under different names such as xxx.good.par


        In order to fix the value PAR of a parameter set its a priori error ERPAR to well below the least significant digit of PAR.  The smallest value that can be used is 1.E-37.

          It is also possible to actually zero the value of a parameter by using the following for the value and the error (PAR and ERPAR):
        0.000000000000000E+000 1.00000000E-037
        This value of PAR stays fixed in straightforward fits, but in complex fits may become incremented to some value in the 10-60 range, which is still pretty low, and will usually be reset to an actual zero on final convergence. 
        An older alternative has been to use something like:
        1.000000000000000E-025 1.00000000E-037
providing that this was not too large for your application.


        In order to omit a line from the fit or to predict a line carrying zero frequency:

  • Use the ERRTST parameter (sixth parameter in second line of .par) but make sure that the (obs-calc)/error for such lines will not fall below ERRTST
  • Assign a very large error to this line e.g. 1.E+20 (but remember to remove such lines from the final fit since their presence artificially lowers the deviation of fit)

        It is possible to assign a zero error to a line and such a line will most of the time be rejected from the fit. This mechanism is not safe, however, since if accidentally the obs-calc deviation for this line falls below a certain threshold (of the order of 0.001) then this line can be captured into the fit and in consequence will distort it.


        In order to run a prediction from known constants but without any lines set up a dummy .lin file with just one line carrying legal quantum numbers and fix all constants in the .par file. Then run SPFIT, which will generatethe .var file, allowing SPCAT to be run.


        The .lin file is read either to the number of lines declared in the .par file or to the first blank line, whichever comes first. There is thus no need to count the number of lines manually: just set an excessively large number in the .par file.

        Uncertain lines/additional information can be kept at the end of the .lin file, after at least one blank line, and will not interfere with running of the program. Furthermore, it is possible to place comments at the end of each line, and this is used for example in the annotation scheme used by the PIFORM program.


        The .par file is acted on up to the number of declared parameters, or to the first blank line, whichever comes first.  However, if there are non-blank lines after the declared parameters, then those will be echoed to the .par file resulting from the fit. This mechanism provides a very useful way of keeping an annotated history of fits in the .par file, including previous versions of this file.   There is a limitation in that only the first 80 columns of such lines are copied over.

WARNINGS:

  • For the 'official' 32-bit SPFIT executable the optional text after the end of the parameter block should not extend beyond 82 columns (even spaces) as this will generally cause SPFIT to crash.
  • There should be no comma before SPIND in the third line of  .par and in its possible continuation lines, and the consequence is a difficult to understand crash.  It is good form, however, to separate all further fields explicitly with commas (especially if those are left blank)
  • An accidental empty line in the .par file will cause a loss of all of the information that was stored below it in the file.

        SPFIT no longer requires access to descriptive names of constants that were taken from the .nam files. It is sufficient to place an appropriate alphanumeric comment at the end of the parameter declaration line and this will be carried over to resulting .par files. Note that this comment has to be separated by a space and a slash, and cannot be longer than 10 characters:

        1200 2.095827797711319E-007 1.00000000E+000 /HJK


        The use of a negative parameter identifier in the .par file has, as the primary function, the declaration of a parameter for which the value follows that of an immediately preceding parameter by a fixed ratio.

        A second feature is that it is also possible to construct a dependent parameter in such a way that it is a linear combination of independent parameters. This is because the values corresponding to the same negative parameter identifier following several different positive identifiers are additive. An example of this is given by HMP in his documentation, as the last of the 'simple examples'.


        In fitting of problems with strong intercorrelations between several constants it is useful to reduce their search range with suitably chosen a priori errors. WARNING: in the final fit such errors should be at least an order of magnitude greater than the errors in the fitted constants.  Otherwise the solution may be artificially constrained and appear more precise than it really is.

     

        Compilation (only for the FORTRAN legacy versions)

        Compilation with NDP-F-486 3.1:

   mf77 -exp -phar2 -OLM -onrc -nomap calfit.for subfit.for
ulib.for blas.for slib.for spinv.for

        The dimensioning in the top PARAMETER statement in CALFIT may be changed as required, and the original versions of routines ULIB, SLIB, SPINV which contain SAVE calls by name of COMMON block are necessary.

        Compilation with SG Fortran:

   f77 -O2 -static -o spfit calfit.for subfit.for ulib.for blas.for
slib.for spinv.for

        Compilation with MS Power-Fortran:

   fl32 -G4 -Ox calfit.for subfit.for ulib.for blas.for slib.for
spinv.for
      

SPCAT = the predictive program


             It is best to run SPFIT before running SPCAT, even if only on dummy data, as this automatically generates the .var file required for an SPCAT run.
             
input:   .var   -   constants, errors + correlation matrix from SPFIT
    .int   -   information for intensity calculation
    .bin   -   binary version of parameters and variances generated by a previous run of SPFIT (this is now obsolete and no longer used)
output:         .out   -   main output file with block by block listing of energies and transition frequencies
    .cat   -   catalogue file, calculated transitions only, 80-column format. The file serves as input for the merging program and is now sorted in frequency.
    .str   -   output of transition dipoles
    .egy   -   output of energy levels, derivatives, and eigenvectors
             

        The .int file requires the rotational partition function Qr which should by default be calculated for 300K. Intensities for other temperatures, such as those of supersonic expansion, could earlier only be recalculated externally, which can be done with PISORT.

       However, it has been possible since at least 2007 to directly calculate intensities for a specified temperature (as this is the ninth parameter in the second line of the .int file).  Remember that it is necessary to not only specify the temperature explicitly but also the Qr value for this temperature. A preliminary run of SPCAT with the desired temperature will precalculate the needed Qr value.  Some commonly used approximations for Qr are given below.

        For a linear top (Gordy & Cook, ed.3 p.57, Eq.3.64):

   Qr = kT/(hB) = 2.0837x104 T/B             ; B/MHz, T/K

        For a (prolate) symmetric top (G&C, Eq.3.68):

   Qr = 5.3311x106 (T 3/(B 2 A))1/2 /sigma      ; A,B/MHz, T/K

        For an asymmetric top (G&C, Eq.3.69):

   Qr = 5.3311x106 (T 3/(A B C))1/2 /sigma     ; A,B,C/MHz, T/K

        The line intensity (i.e. LGSTR.) tabulated in the SPCAT output is the logarithm of intensity Iba(T) as defined in Eq.(1,2) of the catdoc.pdf documentation file for the jpl catalogue, and is in units of nm2MHz. Conversion to peak absorption coefficient αmax in cm-1 is made with Eq.(5) of the documentation, i.e.

   αmax = 102.458 xIba(T)/Δν,

where Δν is the HWHH in MHz at 1 Torr pressure. If the symmetry factor σ and the nuclear spin weight gI are both unity then Iba(T) is a factor of (102.458/20)x0.95=4.87 smaller than αmax calculated by ASROT, which uses Eq.7.73, p.263 G&C, and assumes Δν=20MHz, Fv=0.95. For prediction of intensities for a vibrational state with Fv lower than 1, use Qr/Fv for Qr.

        Dipole moments are to be specified as:

   1  x.xxx     mua/D for state v=0
2 x.xxx mub/D for state v=0
3 x.xxx muc/D for state v=0
111 x.xxx mua/D for state v=1
112 x.xxx mub/D for state v=1 etc.

  

        SPCAT Notes and Hints


       
If parameters in the .par file are in MHz and if you want the energies in the .egy file to also be in MHz (and not in the default cm-1) then in the .int file set FLAGS=1101, for example. The trick is to set the most significant digit, IRFLG=1. Even though this involves pretending that the parameters are in cm-1, the result of associated conversions is that numerical values for the energies come out in MHz.


        Since it is the .var file and not the .par file that is used by SPCAT it is possible to change predictions by editing some control parameters in the .var file without the need for a new fit.



        It is possible to use more than six quantum numbers to describe energy levels.  The examples are asymmetric top cases Q=54 and Q=43 (see the "Format of Quantum Numbers" section in SPINV.PDF).  In such cases all spin quantum numbers preceding F are lumped into an aggregate quantum number nn.  All non-aggregated quantum numbers are only printed in the .egy file, while nn,F are used in the .lin and .cat files, so that the situation is compatible with ASCP_L. Conversion between all spin quantum numbers and nn,F has not been defined.



        Partition function, hyperfine structure and state multiplicity:

        The SPFIT/SPCAT package provides several ways of dealing with predictions for higher-J transitions for a molecule with a quadrupolar atom.  The lowest-J rotational transitions will show nuclear quadrupole hyperfine structure, but in higher-J transitions this will be completely collapsed and is usually ignored.  Databases set up for astrophysicists use several possible approaches so it is useful to discuss the differences.  The HCCCN molecule is used as an example, and it contains a single 14N nucleus, nuclear spin I=1.  The various possibilities for generating line lists are:

CASE 1: Hyperfine structure is ignored.  The molecule is treated as a standard linear rotor, as in the 051501 species entry in the  CDMS database.  Spin is assumed to be zero, hence the spin degeneracy parameter SPIND in the .var file is set to 1 (its value is 2I+1).  Each rotational level is also declared to have multiplicity of 1 (WTPL=1, WTMN=1 in the .var file).  The partition function value QROT in the .int file is then compatible with the linear top formula discussed above and is Q(300K)=1374.9098.

CASE 2: Hyperfine structure is implicit, as treated in the 51001 species entry in the jpl database.  Explicit nuclear spin is ignored as in CASE 1, so that spin degeneracy SPIND is also declared as 1.  But each rotational level is declared to be triply degenerate, using WTPL=3, WTMN=3.  In consequence, there are three times more levels to be counted in evaluation of the partition function and its room temperature value required in the .int file is therefore Q(300K)=4124.495.  Nevertheless, the  calculated transition intensities as given by the LGINT values in the .cat file are identical to those from CASE 1.

CASE 3: Hyperfine structure is explicit.  Nuclear spin degeneracy parameter is set to
SPIND=3 as required for nuclear spin I=1. Each rotational level is an explicit hyperfine level so that there is no hidden degeneracy and WTPL=1, WTMN=1. There are still three times the number of rotational levels to be counted than in CASE 1, so that the partition function is three times larger, and is Q(300K)=4124.7294 (as specified at the bottom of this partition functions listing for species 051501 in CDMS).  The predicted transitions will now contain the hyperfine components explicitly.  For each higher-J rotational transition there will usually be three overlapped components, each of intensity roughly a third of those in CASE 1 and CASE 2, and summing up to those from the first two cases.

     

        Compilation (only for the FORTRAN legacy versions)

        Compilation with NDP-F-486 3.1:

   mf77 -exp -phar2 -OLM -onrc -nomap calcat.for ulib.for
blas.for slib.for spinv.for

        Compilation with MS-Power Fortran:

   fl32 -G4 -Ox calcat.for ulib.for blas.for slib.for spinv.for

        The dimensioning in the top PARAMETER statement in CALCAT may be changed as required and other notes for compilation of SPFIT are also applicable.

   
   

APPENDIX I Parameter indices and names
      
 projection type TYP_  _KSQ = power of K^2
|| _NSQ1 = power of N(N+1)
(I2>=I1) I2 ||| _V2
|I1 |||| _V1 >= V2
|| _|||||
XXxXxXXxXxX
-- |
| NS = power of N.S
FF = Fourier flag
  • For NVIB>9 V1 and V2 become two digit numbers
  • Setting V1=V2=9 for NVIB<10, or V1=V2=99 otherwise, enforces the same value of this parameter for all vibrational states. If the same parameter is then declared with a different specification of V1,V2, the declared value is the difference relative to the global specification e.g. BvB
  • Negative parameter identifier specifies that this parameter is to follow the value for the immediately preceding parameter with a +ve identifier, at the ratio defined by their initial values. This allows fitting of simple two component linear combinations, as well as of blocks of linked values
----------------------------------------------------------------------
Watson-A Watson-S Hyperfine Coriolis (0<->1)
----------------------------------------------------------------------
10000 A         10000 A       spin 1:                    11 E1
20000 B 20000 B 110010000chi.aa*3/2 200001 Ga
30000 C 30000 C 110020000chi.bb*3/2 210001 Fbc
110030000chi.cc*3/2 400001 Gb
200-DeltaJ 200-DJ 110040000chi(b-c)/4 410001 Fca
1100-DeltaJK 1100-DJK 110610000chi.ab 600001 Gc
2000-DeltaK 2000-DK 110210000chi.bc 610001 Fab
40100-deltaJ 40100 d1 110410000chi.ac
41000-deltaK 50000 d2 110011000chi.k*3/2 -----------------
110010100chi.J*3/2 Fermi (0<->1)
300 PhiJ 300 HJ -----------------
1200 PhiJK 1200 HJK spin 2: 01 F
2100 PhiKJ 2100 HKJ 220010000chi.aa*3/2
3000 PhiK 3000 HK 220020000chi.bb*3/2 -----------------
40200 phiJ 40200 h1 220030000chi.cc*3/2 Internal rotation
41100 phiJK 50100 h2 etc. -----------------
42000 phiK 60000 h3 200000 D_a
spin-spin: 200100 D_aJ
400 LJ 400 LJ 120010000 D.aa*3/2 201000 D_aK etc.
1300 LJJK 1300 LJJK 120040000 D(b-c)/4
2200 LJK 2200 LJK 120610000 Dab 400000 D_b
3100 LKKJ 3100 LKKJ 600000 D_c
4000 LK 4000 LK spin-rotation (nucl.1):
40300 lJ 40300 l1 10010000M.aa {+IMJ
41200 lJK 50200 l2 10020000M.bb {+IMJ
42100 lKJ 60100 l3 10030000M.cc {+IMJ
43000 lK 70000 l4
  500 PJ          500 PJ       
1400 PJJK 1400 PJJK
2300 PJK 2300 PJK
3200 PKJ 3200 PKJ
4100 PKKJ 4100 PKKJ
5000 PK 5000 PK
40400 pJ 40400 p1
41300 pJJK 50300 p2
42200 pJK 60200 p3 <- A: up to sextics G&C p.331,
43100 pKKJ 70100 p4 octics+decadics p.333
44000 pK 80000 p5 <- S: up to sextics p.334, octics p.335
----------------------------------------------------------------------

 

   

APPENDIX II Various postprocessing programs for manipulation of data and output for SPFIT and SPCAT available from the PROSPE web site in its asymmetric top programs section at:
http://info.ifpan.edu.pl/~kisiel/asym/asym.htm#pickett





        Most often used:
PIFORM   -   Formatting of the .fit output into a form closer to that customarily used in publications
PISLIN   -   1/ Checks for line duplication and split blends in a .lin file.  2/ Sorts the .lin file according to any quantum number, frequency and magnitude of error (also possible to subsort according to a second criterion within a previous sorting criterion).
LINFIL
-
Filtering and global modifications of a .lin file according to many different criteria.
CATFIL
-
Filtering of a .cat file according to a combination of user specified options.







-
Less frequenctly used:
PICHAM   -   Reformatting of the .egy file to make it easier to check the actual Hamiltonian matrix elements set up for a given parameter in the .par file.
PMIX
PMIXC

-
Graphical identification of resonances between energy levels in interacting vibrational states using mixing coefficients from the .egy file.
PISORT
      Sorting of predictions in .out file according to frequency and optional conversion into .ASR format for the legacy graphics viewer ASCP
LINASF


Converter of a .lin file to frequency block in .asf standard of the ASFIT program.
CATBAND


Generation of line sequence/band markers on diagrams drawn with gle.









Indispensible:
ASCP_L


  -   Graphics viewer of SPCAT predictions as stick diagrams: multiple data sets, various types of highlighting, scrolling etc., reads .cat files directly
   

APPENDIX II SPFIT/SPCAT and Windows Vista (as of July 2008)

         Windows Vista introduced a new protective software layer called User Account Control (UAC). Under default settings the launching of some programs triggers UAC and you are asked to explicitly confirm that they are bona fide, even if you are an administrator. Large (but not small) runs of SPFIT/SPCAT trigger UAC. The problem is that you do not get the chance to do anything and the programs either hang up or crash. Both SPFIT and SPCAT are subject to this, and the objectionable action seems to be the generation of a scratch memory area for larger sorting operations:

SPFIT: crashes while reading in the .LIN file with the console message "scratch file write error"

SPCAT: hangs after the final "STARTING QUANTUM xx" message at the stage when it would commence sorting the .CAT file in frequency

         Switching off UAC altogether cures this problem but it is not recommended. There are various ways of dealing with the problem at the system configuration level, but those also depend on which version of Vista you might have. The simplest and noninvasive way of circumventing this problem is to:

  • place the "Command Prompt" icon on the desktop
  • right-click on the icon and open "Properties"
  • in the "Shortcut" tab click on "Advanced" and tick the box "Run as administrator"

         Launching of this "Command Prompt" window still requires one confirmatory click from UAC, but any programs run from within this window will not be subject to further UAC scrutiny. This special window is identified by the word "Administrator:" at the beginning of the text on the top bar. The same technique should work also for programs that use SPFIT/SPCAT in some batch mode, or you can just run such programs from the prepared "Command Prompt" window.

         If you run your programs from a file manager, such as Total Commander, then an even better solution is to launch the file manager itself in the "Run as administrator" mode. You can read more about UAC here and here.



  

APPENDIX IV The story of how the FORTRAN programs were extended to deal with CH2I2(and NOBr) - this was done on a purely phenomenological basis, some guidelines and explanation of key parameterized variables are now given by HMP at the bottom of SPINV.TXT

              (none of this should be necessary with C-versions of the programs although the largest cases dealt with using the manually extended versions did not want to run with those)

        With standard parameter values (i.e. program with dimensioning as obtained in 1994 from H.M.Pickett) there was a diagnostic

     "spin states truncated in SETSP 18 36"

and for I,F double quadrupole formulation the energy levels for F<I were not calculated.

        --> this diagnostic was traced to limit on parameter NX which was set to NX=17 in several parameter statements in SPINV.FOR. When this was changed to NX=36 there was a diagnostic

     "total number of sub-blocks per F is too large"

        --> this was identified to be associated with MAXQ=100 in several parameter statements in SPINV. On change to MAXQ=200 there was a diagnostic

     "GETQQ: too many subblocks 36"

        --> this was associated with the limit MAXS=36 in several parameter statements in SPINV. On change to MAXS=72 the predictive calculation for CH2I2 went without hiccups up to J=9. Maximum dimension of hamiltonian=576,J=10 (@ HEAPLN=550000), calc. taking over 1 hour on 486/33+NDP Fortran. However J=10 was not calculated.

        --> this was caused by insufficient value of NDMX in SPINV. On change from NDMX=350 to NDMX=700 in SPINV and also on change in CALCAT.FOR from NDEGY=350 to NDEGY=700 and to HEAPLN=1500000 then calculation for CH2I2 appeared possible up to at least J=18.

   
   

APPENDIX V Useful modifications to FORTRAN versions distributed by HMP

               --> In routine FILGET in SLIB.FOR change

         80  FIL='\SPECTRA\'//FNAME
to 80 FIL='..\'//FNAME for DOS or
to 80 FIL='../'//FNAME for UNIX

        so that programs can be used by opening branches for molecules in the program catalogue and no other subdirs are necessary