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