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.