These
are
loose notes based only on applications of Pickett's SPFIT/SPCAT (or
CALFIT/CALCAT) program package tried by kisiel@ifpan.edu.pl. These notes
should only be used once you are familiar with the documentation of the package: - 'official' documentation of the package is contained in http://spec.jpl.nasa.gov/ftp/pub/calpgm/spinv.pdf.
- annotated version of the above is available on the PROSPE website at http://info.ifpan.edu.pl/~kisiel/asym/pickett/spinv15_annotated.pdf.
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. last update: 25.III.2016 Notes: - The use of SPFIT/SPCAT should be
acknowledged in publications by citing: H.M.Pickett,
*J.Mol.Spectrosc.***148**,371-377(1991).
- Current versions of the programs are available from http://spec.jpl.nasa.gov and are collected at http://spec.jpl.nasa.gov/ftp/pub/calpgm/
- Tested 64-bit executables for both Windows and Unix are available in the PICKETT section of the PROSPE website
- The jpl site contains data files for many molecules, which can be treated as worked examples of the use of these programs and are collected in http://spec.jpl.nasa.gov/ftp/pub/catalog/archive/archive_index.html
- You can also browse the
*Fitting Spectra*section of the CDMS site at http://www.astro.uni-koeln.de/cdms/pickett for more information/clarifications/examples on SPFIT/SPCAT
SPFIT
in
the 1994 version worked only for
Third
line
in the I=3/2):s -4 1 0 ,0 , , , , , , , (three
quadrupolar
nuclei, s -334 1 0 ,0 , , , , , , , (three
quadrupolar
nuclei, 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 A s -33 1 0 , 0 , , -1 , 1. , 0. ,,,,,,,,,,, <---- A1 Important constants: 100 B or 20000 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 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 J' F' J'' F'' 0 0 0 0 0 0 0 0 freq error weight
Third
line
in the s -1 1 0 ,,,,,,,,,,,,,,,,,, Important constants (prolate): 10000 A Alternatively: 1000 A-B :prolate Line
in
the J' K' J'' K'' 0 0 0 0 0 0 0 0 freq error weight Note
that
J+1 -3 J 3 0 0 0 0 0 0 0 0 The
sign
of
Treatment
of
Prolate top: s -1, 3, , , , 6, 2, 2, -1, 1, 1, ; this is the g.s. Oblate top: s -1,-3, , , , 6, 2, 2, -1, 1, 1, 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 200011 -2 A zeta.EE 600011 -2 C zeta.EE for
off-diagonal
elements in the form 0.25[ Transitions
which
are normally degenerate are best specified by using positive Prolate
(positive
J' -1 1 J'' 1 1 ....... A+ i.e. upper freq. component Prolate
(positive
J' 1 1 J'' -1 1 ....... A+ i.e. upper freq. component Oblate
(positive
J' -1 1 J'' -1 1 ....... A+ Oblate
(positive
J' -1 1 J'' -1 1 ....... A- If
the
prolate oblate 600001 200001 w = 2 root2 Omega B zeta.EA where Ω = [sqrt( ω)
+
sqrt(_{e}ω/_{e}ω)]/2 and _{a}ζ^{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.
Third
line
in the Reduction- a 1 1 0 ,,,,,,,,,,,,,,,,,, Reduction- a 1 -1 0 ,,,,,,,,,,,,,,,,,, Reduction- s 1 1 0 ,,,,,,,,,,,,,,,,,, Reduction- s 1 -1 0 ,,,,,,,,,,,,,,,,,, Standard
constants
are as given in APPENDIX I. The
difference in the third B+C), (B−C) declared
using:20000 value /(B+C)/2 or, in an alternative way, using: 20000 value /(B+C)/2 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. 10000 value /(A0+A1)/2 Line
in
the J' K-1' K+1' J'' K-1'' K+1'' 0 0 0 0 0 0 freq error weight
The
lower
of the two levels is to be given the identifier 10000 A rotational constant for state v=0 etc. Third
line
in the A-type,
prolate):a 1 2 0 ,,,,,,,,,,,,,,,,,, Line
in
the J' K-1' K+1' v' J'' K-1'' K+1'' v'' 0 0 0 0 freq error weight
Specification
is
identical to the case above with the exception that only 11 Delta E is used. Transition frequency and its error can be specified
in cm
Three
states
need to be specified, ground state a 1 3 0 ,,,,,,,,,,,,,,,,,, and some pertinent constants are: 20000 B for the ground state etc.
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 110030000 1.5 Chi.cc :Oblate 110610000 Chi.ab :Prolate and oblate Third
line
in the a 4 1 0 ,,,,,,,,,,,,,,,,,, Third
line
in the s 4 -1 0 ,,,,,,,,,,,,,,,,,, Line
in
the J' K-1' K+1' F' J'' K-1'' K+1'' F'' 0 0 0 0 freq error weight
Notes from the single nucleus case above are applicable. Prolate rotor constants: 110010000 1.5 Chi.aa \ Third
line
in the a 34 1 0 , , , -1 ,,,,,,,,,,,,,,, Line
in
the I,F
coupling):J' K-1' K+1' I' F' J'' K-1'' K+1'' I'' F'' 0 0 freq error weight
The
D coefficients are for 120010000 1.5 D.aa The first two digits define interacting nuclei. Note
that the constant
The
10010000 M.aa for nucleus 1 For linear molecules the use of: 10000000 M for nucleus 1 is equivalent to the declaration: 10020000 M.bb for nucleus 1 Note
that,
following Flygare, a reverse Hamiltonian definition is often used
i.e.
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 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,
The two RMS quantities in the (MICROWAVE RMS) = sqrt ( Sum[(obs-calc)**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!
In the standard application of
weighted least-squares (section E5 of ASZ) only 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: In this case 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 D 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.
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 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 D 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 _{2}S 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 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 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 - when ERPAR>1.E+015 on
__all__fitted parameters then NDFREE=NLINE-(the number of free parameters)+1 - when ERPAR<=1.E+015 on
__all__fitted parameters then NDFREE=NLINE - 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:
- - - 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: 2/ The THRESH
parameter (the fifth parameter in the second line of the 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 - - - 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 file are negative relative to
the usual convention, as you can check in this
table. The only exception are the two
off-diagonal .parS-reduction
constants, d_{1} and d_{2}. The recommended
solution is to write the parameter descriptions explicitly in
the
file:.par200 -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). 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..parIf the fit goes haywire make sure to recover the original by copying from the .par 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 .bak and .par files under different names such as xxx.good.par.fit In order to fix a value of a constant set its It
also possible to actually zero the value of a constant by using the
following for the value and
the error (PAR
and ERPAR): 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
) but make sure that the (obs-calc)/error for such lines will not fall below ERRTST**.par** - 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. To run a prediction from known constants but without any
lines set up a dummy file. Then run SPFIT, which will generate .par and .bin files, so that SPCAT can then be run..var The
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. .par Uncertain
lines/additional
information can be kept at the end of the PIFORM program. The 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..parWARNINGS: - In the current SPFIT
only 82 columns of the
file are recognised, and lines longer than 82 columns will cause the program to crash.**.par** - An accidental empty line in the
file will cause a loss of all of the information that was stored below it in the file.**.par**
SPFIT no longer requires access to
descriptive names of constants that were taken from the files. Note that this comment has to be separated by a space
and a slash, and cannot be longer than 10 characters:.par1200 2.095827797711319E-007 1.00000000E+000 /HJK The use of a negative parameter identifier in the A second feature is that it is also possible to construct a
dependent parameter in such a way that it is a In fitting of problems with strong intercorrelations between
several constants it is useful to reduce their search range with
suitably chosen
Compilation with NDP-F-486 3.1: mf77 -exp -phar2 -OLM -onrc -nomap calfit.for subfit.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 Compilation with MS Power-Fortran: fl32 -G4 -Ox calfit.for subfit.for ulib.for blas.for slib.for
The Q_{r}
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. It is now
also possible to calculate intensities for a specified temperature with
SPCAT, but remember that, apart from specifying the temperature in the .int file, it is also necessary to calculate Q_{r}
for this temperature. Some commonly used approximations for Q_{r}
are given below.For a linear top (Gordy & Cook, ed.3 p.57, Eq.3.64):
For a (prolate) symmetric top (G&C, Eq.3.68):
For an asymmetric top (G&C, Eq.3.69):
The
line
intensity (i.e. LGSTR.) tabulated
in the SPCAT output is the logarithm of intensity α where Δν is the HWHH in MHz at 1
Torr pressure. If the symmetry factor σ and
the nuclear spin weight Dipole moments are to be specified as: 1 x.xxx mua/D for state v=0
file to also be in MHz (and not in the default cm.egy^{-}^{1}) then in the file set FLAGS=1101, for example. The trick is to set the most significant
digit, IRFLG=1. Even though this involves pretending that the constants
are in cm.int^{-}^{1}, the result
of associated conversions is that numerical values for the energies
come out as if they were in MHz.Since it is the file and not the .var file that is used by SPCAT it is possible to change predictions by editing some control parameters in the .par file without the need for a new fit..varPartition function, hyperfine structure and state multiplicity: The SPFIT/SPCAT package provides several ways of dealing with predictions for higher- I+1). Each rotational level is also declared to have multiplicity of 1 (WTPL=1, WTMN=1 in the file). The partition function value QROT in the .var file is then compatible with the linear top formula discussed above and is .intQ(300K)=1374.9098.
Q(300K)=4124.495. Nevertheless, the calculated transition intensities as given by the LGINT values in the file are identical to those from CASE 1..catCASE 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 with NDP-F-486 3.1: mf77 -exp -phar2 -OLM -onrc -nomap calcat.for ulib.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.
projection type TYP_ _KSQ = power of K^2 - 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 the constant for all vibrational states. If
a given constant is then declared with a different specification of
V1,V2, the declared value is the
*difference*relative to the global specification e.g.*B*_{v}−*B* - 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
---------------------------------------------------------------------- 10000 A 10000 A spin 1: 11 E1 500 PJ 500 PJ
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 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 - 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.
(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.
--> In routine FILGET in SLIB.FOR change 80 FIL='\SPECTRA\'//FNAME so that programs can be used by opening branches for molecules in the program catalogue and no other subdirs are necessary |