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