Instructions for program ERHAM (Version ErhamV16g-R3 5/20/13)
-------------------------------------------------------------
Reference: P. Groner, J. Chem. Phys. (1997) 107, 4483-98.
ERHAM sets up and solves the "Effective rotational Hamiltonian for molecules
with two periodic large-amplitude motions". It allows to fit spectroscopic
constants to observed transition frequencies and to predict the spectrum.
Contents (* indicates changed or new sections)
--------
A. Using Erham *
B. Modifications of previous versions *
C. Input in general
D. Screen Input *
E. File Input *
F. Comment on Output *
G. Relative intensity problem
H. Using ERHAM for molecules with one periodic large-amplitude motion
I. Comment on symmetry parameter ISCD and direction cosine parameter NC
J. S-reduction Hamiltonian and higher order distortion constants *
K. JPL catalog file *
L. Parameter scaling *
M. References *
A. Using Erham
--------------
The ERHAM executable can be started in MS-DOS or by double-clicking the file
icon in the folder where it is located (Windows95 or higher). A DOS window
opens and asks for the names of input and output files. Enter the name (up
to 20 characters) and hit Return.
A version of ERHAM that is compatible with CAAARS 4.0 (I. R. Medvedev et al,
(J. Mol. Struct. 742 (2005) 229-236) is currently tested. Contact the author
(gronerp@umkc.edu) for details.
All input/output files of both versions must be or are text (ASCII) files but
do not need the .txt extension. All input files must be in the same folder/
directory as the executable code. The output files will be placed in the same
folder. The name of the output files are not protected: an existing file with
the same name will be overwritten. The prediction routine opens two scratch
files to store temporary information. They should disappear from the
directories on normal termination of the program.
Least-squares fitting and predictions:
The least-squares fit is performed if the number of iterations is positive,
at least one spectroscopic parameter is declared a variable, and observed
transitions with non-zero weights are provided. (The program does not compare
the number of variables with the number of non-zero weight frequencies.
However, an error message is generated by the Singular Value Decomposition
subroutine.)
Predictions are made after (if any) least-squares fit if JMAX >= JMIN
(Input 4) and at least 1 symmetry block is defined (Input 6). However, no
predictions will be printed if FMX < FMN, THRES is too large (Input 4) or all
dipole components (Input 2) are zero. The parameters marked by *** are only
relevant for the predictions. They have to be entered anyway because of the
free-format input. The partition function is calculated with the predictions,
but only if JMIN = 0 and only up to JMAX (see also Section F).
B. Modifications of previous versions:
--------------------------------------
Version V16g
13 apr 05: correction ISCD test in subroutine INPUT
Modifications in preparation for use with CAAARS
28 jun 05: unified output and format change for predictions (ORDER)
format change fitted frequencies (ITER)
format change JMIN, JMAX (INPUT)
1 jul 05: problem fix in calc. of int. rot. parameters (single rotor)
3 jul 05: dimensions changed to accommodate J < 120
4 & 7 jul 05: changes to allow for S-reduction Hamiltonians
7 jul 05: verified code for variable IMG
problem fix in calc. of int. rot. parameters (en. diff.)
problem fix in calc. of int. rot. parameters (alphaq error,
created on 3 jul 05)
6 feb 06: changed dimensions to allow up to 8191 transitions
some tests for tunneling parameters updated (INPUT)
corrected code for fitting of blends (INPUT)
copied summary of ISCD cases from defunct V16b
11 feb 06: corrections for S-reduction Hamiltonian
31 may 06: changes already made in V15hd transferred into this version:
use of vib. state quantum numbers (PREDIC, ORDER, INPUT)
31 may 06: integrated KVQ into this code
Distinction between CAAARS-compatible and non-CAAARS versions:
At the beginning of ERHAM, ORDER, and INPUT are sections marked
by c&&&& where the versions differ.
1 jun 06: correction print format (ORDER)
5 jun 06: new code for blends (ITER)
6 jun 06: built in new definition for blends for CAAARS version (INPUT)
12 jun 06: removed the c&&&& code pertaining to CAAARS
Rev. R1:
7 jan 09: corrected tunneling parameter matrix elements of Px, Pz
12 oct 09: corrected matrix elements of Pz^n in CD constants
Rev. R3:
20 may 13: added option to print predictions in JPL catalog file format as
separate file (section K), corrected formula for intensities,
threshold now applies to calculated relative intensity =
S*mu^2*B*nu (mu = dipole moment, B = Boltzmann factor, nu =
transition frequency)[previously to S*mu^2], changed dimensions
in ORDER to 140000 transitions, added calculation of partition
function (with respect to lowest torsional sublevel), increased
number of tunneling parameters per state to 37, introduced
parameter scaling for least-squares fit (section L)
C. Input in general
-------------------
ERHAM performs several tests on many input variables. This process often
results in error messages; in other instances, incompatible input is simply
ignored. A very important parameter is ISCD (input 2). Acceptable values for
many variables depend on it. Some input incompatible with ISCD (e.g. non-zero
ALPHA1, ALPHA2 for ISCD=-1,-2) is ignored or overwritten. Also the allowed
combinations of IQ1, IQ2 (input 5), IS1, IS2 (input 6, 7) depend on ISCD. A
comparison of the input file with print of the input data in the output file
shows the overwrites by the program.
D. Screen Input
---------------
1. Respond to the prompt 'Enter input file name !' by typing: FILEIN
FILEIN = name of input file (unit 5, up to 20 characters)
2. Respond to the prompt 'Enter output file name !' by typing: FILEOUT
FILEOUT = name of output file (unit 6, up to 20 characters)
3. If IFPR = 4 (Input 2 below for printing JPL catalog file), catalog file
information is also entered as screen input (section K)
E. File Input
-------------
Note: all input is in file (filename: FILEIN, unit: 5)
1. WT(I),I=1,10
problem information (up to 80 characters, read in (1X,10A8) format; the
first character on this line is not read)
Everything from here on is input in free format.
2. ISCD,NC,NIV,NIT,IFPR,NUNC,(DIP(I),I=1,3),TEMP,KVQ
ISCD = symmetry parameter
ISCD = 1 => C'1 (non-equivalent) = -1 => C's (non-equivalent)
2 => C2 (equivalent) = -2 => C2v (equivalent)
3 => Cs (equivalent) = -3 => C2v (equivalent)
4 => Ci (equivalent) = -4 => C2h (equivalent)
0 => C's (non-equivalent, in bc plane)
-3 => C2v (equivalent, in bc plane)
-5 => C2h (equivalent, in bc plane)
NC = direction cosine parameter (with respect to principal a axis)
NC = +1 ==> direction cosines of internal rotors are equal
NC = -1 ==> direction cosines of internal rotors are opposite
NIV = number of vibrational states (min=1, max=6)
NIT = number of iterations (NIT.LE.0 => no iterations)
IFPR = print option parameter
IFPR = 0 no additional printout
= 1 prints derivatives, a LOT of printout !!!!
= 2 prints eigenvectors of rotational energy levels
= 3 prints eigenvectors and derivatives
= 4 predictions also in JPL catalog file format
*** NUNC = uncertainties parameter
NUNC.NE.0 : calculate uncertainties of predicted frequencies
(possible only after a least squares fit has been
performed in the same run)
NUNC = 0 : no uncertainties of predicted frequencies
*** DIP = components of electric dipole moment (mu(a),mu(b),mu(c)) (Debye)
*** TEMP = temperature in K(elvin) for Boltzmann population factor
KVQ = switch for vibrational label (only in CAAARS-compatible version)
KVQ = 0 : don't read vibrational states label, exp. uncertainty
AFTER blend coefficient in frequency input (old style)
KVQ = 1 : read vibrational states label, exp. uncertainty BEFORE
blend coefficient in frequency input
3a. N1,RHO1,IVRHO1,BETA1,IVBET1,ALPHA1,IVALP1
3b. N2,RHO2,IVRHO2,BETA2,IVBET2,ALPHA2,IVALP2
Note: Input 3b. is needed only if rotors are NOT equivalent (ISCD=-1,0,1)
N1,N2 = periodicities of internal rotors (0 < N1, 0 < N2)
program is aborted if N1 or N2 is 0 or negative
RHO1, RHO2 = internal rotation parameters (dimensionless)
IVRHO1, IVRHO2 = variation parameter for RHO
IVRHO = 0 : keep RHO constant
IVRHO = 1 : RHO is a variable in least-squares fit
IVRHO2 = 2 : RHO2 is equal to RHO1 (non-equivalent rotors)
BETA1, BETA2 = RHO axis polar angle between RHO axis and
"a" principal axis (degrees)
IVBET1, IVBET2 = variation parameter for BETA
IVBET = 0 : keep BETA constant
IVBET = 1 : BETA is a variable in least-squares fit
IVBET2 = 2 : BETA2 is equal to BETA1 (non-equivalent rotors)
ALPHA1, ALPHA2 = RHO axis angle of RHO axis w.r.t. "ab" principal plane
IVALP1, IVALP2 = variation parameter for ALPHA
IVALP = 0 : keep ALPHA constant
IVALP = 1 : ALPHA is a variable in least-squares fit
IVALP2 = 2 : ALPHA2 is equal to +/- ALPHA1 (non-equivalent rotors)
The remaining input (4, 5, 6, 7) is repeated for each vibrational state IV
4. JMIN(IV),JMAX(IV),FMN(IV),FMX(IV),THRES(IV),(A(I,IV),I=5,21),
(IVR(I,IV),I=5,21)
*** JMIN, JMAX = range of J quantum number for predictions ( < 121 )
If no predictions are wanted, use JMIN > JMAX.
*** FMN, FMX = lower and upper limit of frequencies for predictions (MHz)
*** THRES = intensity threshold for predictions (applies to line strength S)
array A : 3 rotational constants (MHz)
5 quartic centrifugal distortion constants (kHz) (A-reduction)
7 sextic centrifugal distortion constants (Hz) (A-reduction)
array IVR : variation parameters for rotational and distortion constants
= 0 : keep constant at input value
> 0 : variable during least-squares fit
< 0 : constant is kept identical to corresponding parameter of
the vibrational state for IV = ABS(IVR)
Input of tunneling parameters (less than 38 per state):
5. IQ1,IQ2,MEG,KAP,JP,KP,PAR,IVAR,SCPP
IQ1,IQ2 = indicate to which localized state a tunneling parameter is
related to
max(IQ1,IQ2) = 7, min(IQ1) = 0
min(IQ2) = 0 for equiv, min(IQ2) = -7 for nonequiv motions
MEG = omega (omega = +1 for even order of angular momenta,
-1 odd ,
0 terminates input of tunneling parameters)
KAP = distance of matrix element from main diagonal (KAPPA = KQ - K)
KAP>0 real, imaginary part of matrix element for MEG = +1, -1, resp.
KAP<0 imaginary, real part of matrix element for MEG = +1, -1, resp.
JP = exponent of P operator (total angular momentum)
KP = exponent of Pz operator (component of P along z)
PAR = value of tunneling parameter (MHz)
IVAR = variation of tunneling parameter
=+1 : yes
=-1 : yes, same value as tunneling parameter with IQ1 and IQ2
interchanged
else: no
SCPP = scaling parameter (approx. uncertainty of parameter, see Section L)
comment: the vibrational energy tunneling parameters have
KAP = 0, JP = 0, KP = 0, MEG = +1
Repeat this input for each parameter you wish to define. Input 5 is
terminated by entering 0 for each of the variables.
Input of symmetry blocks for predictions:
6. IS1,IS2,ISP1,ISP2
*** IS1,IS2 = symmetry labels of symmetry block
*** ISP1,ISP2 = spin multiplicity for (sym) and (asym), resp.
For two-rotor molecules with C2 or C2v symmetry, rotational wavefunctions
are either symmetric w.r.t. the C2 symmetry axis (J=0), or antisymmetric.
ISP1 and ISP2 refer to the spin multiplicities of the symmetric and anti-
symmetric rotational wave functions, respectively.
Repeat this input for each symmetry block you need predictions for. Input
6 is terminated by entering -1 for IS1 and 0 for the other variables.
Input of transitions:
7. IS1,IS2,JQ,NQ,J,NN,FREQ,BLE,ER
IS1, IS2 = symmetry labels of transition
JQ, NQ = J quantum number and label of energy level of upper state
J, NN = J quantum number and label of energy level of lower state
FREQ = frequency of transition (MHz)
BLE = coefficient in linear combinations for blends
BLE = 0. ==> no blend
BLE > 0. ==> transition(s) of blended line
BLE < 0. ==> last transition of blended line
Only one frequency (the weighted average with weight equal
to ABS(BLE)/sum(ABS(BLE)) ) is fit.
ER = estimated uncertainty of transition (MHz) (its inverse square
becomes the weight of the transition in the least-square fit)
ER = 0. ==> zero weight in fit for this transition
Repeat this input for each transition. Input 7 is terminated by entering
-1 for IS1 and 0 for the other variables.
Maximum number of transitions: 8190
F. Comment on Output
--------------------
1. Sometimes, the error message
"0NO CONVERGENCE TOOK PLACE IN COMPUTING EIGENVECTOR xx"
is generated by the diagonalization routine. This is usually not something
to be worried about. It likely occurs when there are near-degenerate
eigenvalues. So far, I have not encountered serious problems. But the
potential for problems is there. Should you suspect such a problem, please
notify me. Eventually, I might have to switch to a different
diagonalization routine.
2. At this point, the information in the predictions is not labeled.
Therefore, the missing information follows here:
For each symmetry block (pair of symmetry labels), the values of the
rotational energy levels are given for each value of J in increasing order
(in MHz). Then the following information is printed for the R- and Q-
transitions:
- J', N', J", N",
- frequency,
- spin weight,
- line strengths (S) for the mu(a), mu(b) and mu(c) components,
- the relative intensity (S * mu(x)^2 * spin weight * Boltzmann factor) * nu
- one-sigma estimated uncertainty.
All transitions are printed with the frequencies in increasing order with
the following information:
- symmetry labels, v', J', N', Ka', Kc', v", J", N", Ka", Kc",
v' and v" refer to the vibrational state labels IV (Input 4)
- frequency,
- one-sigma uncertainty,
- spin weight,
- line strengths (S) for the mu(a), mu(b) and mu(c) components of the dipole
moments
- the relative intensity (S * mu(x)^2 * spin weight * Boltzmann factor) * nu,
- upper state energy (cm^-1)]
Note, that the Ka, Kc label are based on the (assumed) order of energy levels
for rigid asymmetric rotors. These labels may be incorrect when serious
rotational-torsional interactions occur.
Some labels N' or N" are printed as negative numbers in some cases, depending
on the torsional substate and the presence of symmetry operators. A - sign
indicates that the eigenfunction (EF) is antisymmetric under a C_2 operation
which may be a proper C_2 symmetry axis or a C_2 axis perpendicular to a
proper symmetry plane.
* For ISCD = 2,-2,-3,-4,-5 in torsional substates with IS1 = IS2 (including
IS1 = IS2 = 0), - signs mark EFs antisymmetric w.r.t. the proper C_2
symmetry axis.
* For ISCD = 3 (IS2 = N1 - IS1 and IS1 = IS2 = 0 substates) and for ISCD = -1
(IS1 = IS2 = 0 substate), - signs mark EFs antisymmetric w.r.t. the C_2
axis perpendicular to the symmetry plane.
* For ISCD = -2, -3 (IS2 = N1 - IS1 substates but not IS1 = IS2 = 0), - signs
mark EFs antisymmetric w.r.t. the C_2 axis perpendicular to the symmetry
plane that interchanges the internal rotors.
3. The partition function (sum of states) is calculated during predictions
but only if JMIN=0 and for levels up to JMAX. Contributions from rotational
energy levels and torsional substates, but not from excited torsional or
vibrational states. To assess convergence, two accumulated sums are printed
for each symmetry block and each J (beginning with J=0, up to JMAX). The
1st sum is from the levels without the (-) sign mentioned above, the 2nd sum
(which may be zero) from the levels with this sign. After the blank line
follow the sums for JMAX adjusted for the lowest energy level. The final
partition function (up to JMAX) is the sum of all sums in this last row.
G. Relative intensity problem
-----------------------------
Particularly for high J, rotational energy levels may be essentially
degenerate (Ka- or Kc-degeneracy). If symmetrization and factorization are
not allowed, the diagonalization routine generates random phase differences
between the degenerate eigenfunctions. As a result, the (identical) line
strengths of two degenerate transitions is spread randomly among four
degenerate transitions, two with the proper Ka,Kc selection rules, and two
with incorrect selection rules. I don't have an easy solution to the problem.
H. Using ERHAM for molecules with one periodic large-amplitude motion
---------------------------------------------------------------------
ERHAM can be used for molecules with only one internal rotor.
Input 2: Use ISCD = -1 or +1
Use either
Input 3a: define the internal rotation parameters
Input 3b: use 1, .00000001, 0, 0.00005, 0, 0.0005, 0
(small nonzero values prevent problems in calc. of derived parameters)
Input 5: IQ2 must be 0
Input 6: IS2 must be 0
Input 7: IS2 must be 0
Or
Input 3a: use 1, .00000001, 0, 0.00005, 0, 0.0005, 0
(small nonzero values prevent problems in calc. of derived parameters)
Input 3b: define the internal rotation parameters
Input 5: IQ1 must be 0
Input 6: IS1 must be 0
Input 7: IS1 must be 0
I. Comment on symmetry parameter ISCD and direction cosine parameter NC
-----------------------------------------------------------------------
Less than half of all cases programmed have been used. The following cases
seem to work (references in parentheses):
ISCD NC Example
-2 C2v(equivalent) -1 C2=b dimethylether, acetone (2,3,8,9)
+1 C2=a 3-methyl-1,2-butadiene (4)
-1 C's(non-equivalent) -1 Cs=ab ethylmethylether (5), acetone-13C (8)
+1 should also work for ethylmethylether
-1 Cs=ab methylcarbamate (one rotor) (7)
2 C2(equivalent) -1 C2=b dimethyldiselenide (6)
3 Cs(equivalent) -1 Cs=bc 2-fluoropropane (unpublished)
1 C'1(non-equivalent) -1 dimethyldiselenide (asym) (6)
+1 should also work in these cases
None of the other cases have been tried yet.
When NC = -1, the values used and shown for BETA2 are in fact
BETA2' = pi - BETA2
When NC = +1, the values used and shown for ALPHA2 are in fact
ALPHA2' = pi + ALPHA2
J. S-reduction Hamiltonian and higher order distortion constants
----------------------------------------------------------------
This feature has been tested and used successfully.(10,11)
An S-reduction Hamiltonian can be used as follows:
Enter 0 for the A-reduction parameters deltaK, phiJ, phiJK, phiK. Use the
tunneling parameter input (Input 5.) to define parameters in the S-reduction:
d2 : IQ1 = IQ2 = 0, MEG = 1, KAP = 4, JP = KP = 0
h2 : IQ1 = IQ2 = 0, MEG = 1, KAP = 4, JP = 2, KP = 0
h1 : IQ1 = IQ2 = 0, MEG = 1, KAP = 6, JP = KP = 0
Then the parameters printed as DeltaJ, DeltaJK, DeltaK, deltaJ, PhiJ, PhiJK,
PhiKJ, PhiK and phiJ are in fact DJ, DJK, DK, -dJ, HJ, HJK, HKJ, HK and h1
instead, respectively.
Input 5 can also be used to define higher order distortion constants by
using IQ1 = IQ2 = 0, MEG = 1, and appropriate values for KAP, JP and KP.
K. JPL catalog file
-------------------
The catalog file in JPL file format is requested using IFPR = 4 (Input 2)
In this case, additional information is entered interactively during the
prediction as follows:
-Respond to prompt 'Enter catalog file name !' by typing: FILECAT
FILECAT = name of catalog file (unit 7, up to 20 characters)
-Respond to prompt 'Enter catalog ID !' by typing: ICATID
ICATID = catalog ID of molecule
-Respond to prompt 'Enter Partition function !' by typing: Q
Q = partition function at selected temperature TEMP (Input 2)
-Respond to the prompt 'Enter first LOG Intensity Cutoff (LOGSTR0) !' by
typing: LOGSTR0 (as in catdoc.pdf from JPL web site)
-Respond to the prompt 'Enter second LOG Intensity Cutoff (LOGSTR1) !' by
typing: LOGSTR1 (as in catdoc.pdf from JPL web site)
ERHAM prints a summary of this input just before the list of ordered
transitions.
-Please note: The information about degrees of freedom (DR = 3) and quantum
number format (QNFMT = 1404) is coded directly into the program. The 4th
quantum numbers for upper and lower state are used to enter the symmetry
numbers (IS1 & IS2 = 4th QN of upper & lower state, respectively). Q
calculated by ERHAM is NOT used to calculate the catalog entries (allowing
adjustment of Q for vibrational contributions).
L. Parameter scaling
--------------------
Sometimes, "error code = 33" is printed by subroutines LSVDB or LSVDF in the
Singular Value Decomposition, indicating a singular or ill-conditioned normal
matrix [e.g. when a) parameters are linearly dependent or b) when the ratio
of largest singular value to smallest singular value exceeds the precision of
a double precision floating point number (example: rotational and octic CD
constants when expressed in MHz]. Case b) can be avoided by proper scaling of
parameters. Scaling factors don't need to be precise, orders of magnitude are
sufficient. Example: If uncertainties of rotational constants are of the
order 10^(-3) MHz and the uncertainty of an octic constant or a tunneling
parameter is of the order 10^(-12) MHz, a scaling factor in the range 10^(-6)
to 10^(-12) is sufficient.
Scaling factors need to be entered only for parameters defined through the
tunneling parameter input (Input 5). They are set by the program for the
other spectroscopic constants. PS: The entry of quartic and sextic distortion
constants in kHz and Hz, resp. (as opposed to MHz for rotational constants)
accomplishes scaling automatically.
M. References
-------------
1. Theory: P Groner, J. Chem. Phys. (1997) 107, 4483-4498
Contains an example where frequencies from several non-interacting
vibrational states are used in a global fit to determine the global RHO
and BETA parameters with separate rotational, distortion and tunneling
parameters for each state.
Progress report & and mini-review: P Groner, J. Mol. Spectrosc. (2012)
278, 52-67
2. Dimethyl ether: P Groner, S Albert, E Herbst, F C De Lucia, Astrophys. J.
(1998) 500, 1059-1063
Torsional ground state, >1600 frequencies, 22 parameters, fit to
experimental uncertainty
3. Acetone: P Groner, S Albert, E Herbst, F C De Lucia, F J Lovas, B J
Drouin, J C Pearson, Astrophys. J. Suppl. Ser. (2002) 142, 145-151
Torsional ground state, >1000 frequencies, 28 parameters, fit to almost
experimental uncertainty
4. 3-methyl-1,2-butadiene: S Bell, P Groner, G A Guirgis, J R Durig, J. Phys.
Chem. A (2000) 104, 514-520
Another example with a multistate fit
5. Ethyl methyl ether: U Fuchs, G Winnewisser, P Groner, F C De Lucia, E.
Herbst, Astrophys. J. Suppl. Ser. (2003) 144, 277-286
Non-equivalent rotors, torsional ground state, >1000 frequencies, 18
parameters, fit to experimental uncertainty
6. Dimethyl diselenide: P Groner, C W Gillies, J Z Gillies, Y Zhang, E Block,
J. Mol. Spectrosc. (2004) 226, 161-181
C2 symmetry (equivalent rotors) and C1 symmetry (non-equivalent rotors)
7. Methyl carbamate: P Groner, M Winnewisser, I Medvedev, F C De Lucia, E
Herbst, K V L N Sastry, Astrophys. J. Suppl. Ser. 169 (2007) 28-36
One internal rotor, >6000 frequencies, 28 parameters, fit to
experimental uncertainty
8. Acetone-13C: F J Lovas, P Groner, J. Mol. Spectrosc. (2006) 236, 173-177
C2v symmetry (equivalent) and Cs symmetry (non-equivalent)
9. Acetone 1st torsional excited state: P Groner, E Herbst, F C De Lucia, B J
Drouin, H Mäder, J. Mol. Struct. (2006) 795 173-178
10. Octic CD constants: B J Drouin, J C Pearson, A Walters, V Lattanzi, J.
Mol. Spectrosc. (2006) 240, 227-237
11. S-reduction: G S Grubbs II, P Groner, S E Novick, S A Cooke, J. Mol.
Spectrosc. (2012) 280, 21-26.