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