C ************************************************** C * Program BARRIER Version 22 jun 2019 * C ************************************************** C C Program to determine the potential barrier for a single simple C internal rotor from torsional transitions or splittings (e.g. from C ERHAM). It is a simpler offspring of program ASTOR described in C P. Groner, et al., J. Mol. Struct. 142 (1986) 363-366. C C Developed by Peter Groner, University of Missouri - Kansas City C C * Solves Hamiltonian of an asymmetric internal rotor (ASymmetric C TORsion). C * Calculates energy levels, eigenvectors, transition frequencies, C simple torsional transition moments C * Fits energy level differences (torsional frequencies and/or C splittings) to observed ones by adjusting the initial estimates C of the potential coefficients by an iterative weighted least C squares fit. C C This version of the program has been created primarily to derive C barriers to internal rotation for a single rotor from the tunneling C splittings obtained by program ERHAM. The tunneling energy splittings C are printed by ERHAM after the least-squares fit, and labeled by C TORSIONAL ENERGY DIFFERENCES which are relative to A (sigma = 0) C level. For a CH3 internal rotor, it is always equal to E(E) - E(A). C These instructions are intentionally vague about asymmetric internal C rotors (AIR). BARRIER could, but should not, be used for an AIR C because (1) such a problem needs much more info than tunneling C splittings; (2) some features specific for asymmetric rotors are C missing; (3) asymmetric internal rotor potentials are MUCH MUCH more C difficult to determine, particularly if different conformations are C present. C C C A Theory and methods C ********************** C C Hamiltonian H = p F(tau) p + V(tau) C C internal rotation coordinate tau. C F(tau) is the "F-number" function (= internal rotation constant) as C function of tau. C V(tau) is the potential energy function. C C F(tau) = F0 + sum [ F_k*cos(k*tau) + F'_k*sin(k*tau) ], k = 1,6 C V(tau) = (1/2) sum [ V_k*(1-cos(k*tau)) + V'_k*sin(k*tau) ], k = 1,6 C C The coefficients F0, F_k and F'_k have to be determined from a model C of the structure as a function of tau. They cannot be varied C during the least-squares fit. C The coefficients V_k and V'_k can be adjusted during the fit. C C The Hamiltonian is solved in the complex exponential basis set C |k> = (1/SQRT(2*PI)) exp(i*k*tau) C by diagonalization of the Hermitian matrix (Givens-Householder, C bisection and inverse iteration (Wieland) of the eigenvectors). C C C B Input/Output C **************** C C The input and output files are text files; their names are entered C once the program is running. C C C INPUT DESCRIPTION FOR PROGRAM ASTOR C *************************************** C C 1) Format(9A8): Text information (72 columns) C C Note: All other input lines are free format C C 2) IPRT,ISM,ITM,IT1,NEIG,IMAX C IPRT: print parameter C IPRT = 0 : standard printout C IPRT = 1 : extended printout including eigenvectors C and derivatives of energy levels w.r.t. C potential coefficients C ISM: symmetry parameter ( IABS(ISM) < 4 ) (--> remarks) C ISM < 0 : frame and internal rotor have symmetry planes C ISM > 0 : frame or internal rotor has no symmetry plane C IABS(ISM) = P : internal rotor is P-fold, P = 1,2,3 C ITM: number of iterations ( 0 < ITM < 6 ) (--> remarks) C IT1: iteration parameter (--> remarks) C IT1 = 1 : calculated changes to potential parameters C are divided by two (slows convergence, may C prevent divergent fit) C else : regular weighted least-squares process C NEIG: number of calculated energy levels ( 0 < NEIG < 25 ) C IMAX: number of basis functions ( 14 < IMAX < 82 ) C Note: Default values are used for some of these parameters if input C is out of range !!!!!!!!!! C C 3) F0, F_1, F_2, F_3, F_4, F_5, F_6, C F'_1, F'_2, F'_3, F'_4, F'_5, F'_6: C Coefficients of F-number expansion (all in cm-1) C C 4) I1,I2,FR,S,M: transition information for iterations C I1 : label number of lower energy level C I2 : label number of upper energy level C FR : transition frequency or energy level difference C S : estimated error of FR C If S = 0, transition/difference is not used in the fit C M : unit parameter (see remarks) C M = 1 : FR and S in MHz C else : FR and S in cm-1 C C Note : Repeat Input 4) until all transitions are in. (Maximum 24) C If there are LESS THAN 24 transitions and/or splittings, C a dummy transition with zeros for each of I1,I2,FR,S,M marks C the end of the list of transitions C C 5) V_1, V_2, V_3, V_4, V_5, V_6 C V'_1, V'_2, V'_3, V'_4, V'_5, V'_6: C Coefficients of potential function expansion (all in cm-1) C C 6) IC_k (k = 1,6 for V_k; k = 7,12 for V'_(k-6) C IC(k) = 1 : the k-th coefficient is varied C IC(k) = 0 : the k-th coefficient is kept constant C C C D General C *********** ....5...10....5...20....5...30....5...40....5...50....5...60....5...70.. C C At the beginning, one has to predict the transitions or differences C for any given set of potential coefficients before transitions or C differences are provided to determine potential coefficients because C one has to establish first the correct assignment of the transitions C (particularly true for problems with more than one conformation). C This is done by just using zeros for each of I1,I2,FR,S,M for a dummy C transition. This must be done also any time a new set of initial C potential parameters is used. C If one needs to continue an iteration process from a previous run, C the potential coefficients printed during the last (previous) cycle C have to be used as well as the correct assignment. Just changing the C potential coefficients won't do because assignments may have changed C since your last input. C Remember: GIGO (garbage in, garbage out). C C C E Printout C ************ C C The first part shows what's been put in: problem information, IPRT, C ISM, ITM, IT1, NEIG, IMAX. The cosine, sine coefficients of the F C number are reprinted in the order they are put in. C If provided on input, transitions are reprinted (energy levels C involved, frequency, estimated error, unit) C C The following is printed during each iteration cycle. This even C applies when no iteration is performed. C C The coefficients of the potential energy are printed with their C current values and whether they are varied(1) or constant(0). The C energy eigenvalues are given with numbering and value. Numbering C starts at the lowest level with 1. For IPRT = 1, eigenvectors with C respect to the normalized exponential basis set |m> are printed in C double columns, the first and second column containing the real and C imaginary parts, resp. The first coefficient belongs to C m = -(IMAX-1)/2, with m increasing down the column. C If both frame and top have symmetry planes, the symmetry of the C eigenfunctions can be deduced from the components of the eigen- C vectors: their real parts should be symmetric or antisymmetric with C respect to |0> in the middle of each column. C C Transitions are printed in a table. For each pair of energy levels, C the entry consists of the predicted transition frequency (cm-1) and C the squares of the transition moments C and , with n = abs(ISM), C in that order. Use it to (i) read predicted transition frequencies; C (ii) check the relative symmetry of the energy levels (for ISM < 1, C the nonzero cos and sine transition moments indicate that the wave C functions have identical or opposite symmetries (parities), resp., C and that transitions are Raman or IR allowed, resp.) C C For IPRT = 1, derivatives of energy levels with respect the varied C potential coefficients are printed next. C Then follows the table of transitions with energy levels numbers C (check them against the table with the transition moments), observed C and calculated frequencies, residuals, weighted residuals. Next are C the eigenvalues and -vectors of the normal matrix (inverse of C variance-covariance matrix) and weighted standard deviation, followed C by current values of potential coefficients, their standard errors, C calculated changes and the precision for quotation required to keep C the weighted differences within the (weighted) standard deviation. C (The new parameter values appear at the beginning of the next cycle.) C C F Remarks C *********** C C Parameter ISM C ------------- C ISM overrides any incompatible information. C If ISM < 0, all V'_k, F'_k are set to 0 (the structure with tau = 0 C has a symmetry plane). C If IABS(ISM) = L, all coefficients V_k, V'_k, F_k and F'_k where C k is not an integer multiple of L are set to 0. C C ISM = 0: Program may blow up (depends on how computer handles ISIGN C function with 0 argument!) C ISM = 1 or -1: 1-fold internal rotor (not very interesting) C ISM = 2 or -2: 2-fold internal rotor C ISM = 3 or -3: 3-fold internal rotor C If ISM = 3 or -3, the (degenerate) E levels are calculated only once. C Thus, for a common 3-fold rotor, the levels 1,2,3,4,5,6,7,8 etc. are C A,E,E,A,A,E,E,A, respectively, instead of A,E,E,E,E,A,A,E,E,E,E,A. C C 4-fold internal rotor: Use ISM = 2 or -2, dominant terms: F0, V_4. C (also allowed F_4, and if ISM = 2: V'_4, F'_4). C 5-fold internal rotor: Use ISM = 1 or -1, dominant terms: F0, V_5. C (also allowed F_5, and if ISM = 1: V'_5, F'_5). C 6-fold internal rotor: Use ISM = 3 or -3, dominant terms: F0, V_6. C (also allowed F_6, and if ISM = 3: V'_6, F'_6). C C Parameter ITM C ------------- C The maximum number of iteration cycles is 5. This limit was set low C on purpose because asymmetric internal rotor fits, particularly if C different conformations are present, can be pathologically non- C linear and the LS fit is prone to blow up. (see GIGO under General) C C Parameter IT1 C ------------- C It was introduced as a quick fix to a convergence problem (slows C the iteration of course). C