Procedures

Add partial cross-sections to accumulated cross-sections

check the status after allocation

calculates the (R*2)centrifugal matrix from the second term of Eq. 3 in "What are coupled equations?" section; Matrix elements are given in Eq. 4 and 6 of "Coupling Matrix" secion

Combines the contribution from the interaction potential, total and and internal energy (pes_matrix_) with centrifugal matrix \( W_{\mathrm{N}} = V_{\mathrm{N}} + 1/R^{2} L^{2} \)

calculates diagonal element of the centrifgual matrix, see Eq. 4 in "Coupling Matrix" section

calculates the K-matrix from log-derivative matrix using Eq. (4) in "Solution of coupled equations"

Calculates the largest wavevector in the block; called only if there are any open channels

calculates the log-derivative matrix from $$\begin{equation} {Y}_{\rm N} = \frac{1}{h} \Biggl(\Bigl(\frac{1}{2}\mathbf{I}-{T}_{\rm{N}+1}\Bigr) \Bigl(\mathbf{I}-{T}_{\rm{N}+1}\Bigr)^{-1} {R}_{\rm{N}+1} - \Bigl(\frac{1}{2}\mathbf{I}-{T}_{\rm{N}-1}\Bigr) \Bigl(\mathbf{I}-\mathbf{T}_{\rm{N}-1}\Bigr)^{-1}\mathbf{R}_{\rm{N}}^{-1} \Biggr)\Bigl(\mathbf{I}-{T}_{\rm{N}}\Bigr) \end{equation}$$ called by numerov at the end of the propagation

calculate number of channels in even and odd parity blocks in the body-fixed frame; in principle, \(\bar{\Omega}\in \langle 0, \mathrm{min}(j, J))), but number of channels additionally depends on the sign of \( p (-1)^{J} \): channels with \(\bar{\Omega}=0\) values enter blocks with \( p (-1)^{J} = + 1 \) only;

calculate number of channels in even and odd parity blocks in the space-fixed frame based on available values of orbital angular momentum: \( l \in \langle |j-J|, j+J \rangle \); parity is defined as \(p= (-1)^{j+l}\)

Calculates the number of steps on the intermolecular (R) grid in current block. This is done either directly (if r_step > 0) or through the number of steps per half de Broglie wavelength

calculates off-diagonal element of the centrifgual matrix, see Eq. 5 in "Coupling Matrix" section

calculates the contribution to the coupling matrix from the the interaction potential (PES); see Eq. 1 in "Coupling Matrix" section; diagonal contribution from wavevectors (see the last term in Eq. 3 of "What are coupled equations" section) is added

calculates S-matrix from open-open portion of the K-matrix using Eq. (12) in "Solution of coupled equations"

takes as an input matrix in the body-fixed frame and transforms it to the spec-fixed frame; iterates over all matrix elements and calls calculate_single_SF_element

Implementation of Eq. 1 in "Coupling Matrix" section; diagonal contribution from wavevectors (see the last term in Eq. 3 of "What are coupled equations" section) is added

calculates a single space-fixed matrix element from Eq. (2) in "Solution of coupled equations"

Calculates all state-to-state cross-sections.

calculates the sum \( \sum_{\gamma'} \Bigl|{S}^{Jp}_{\gamma, \gamma'}\Bigr|^{2} \) for all \(\gamma\) channels

Calculates the T-matrix from the coupling matrix at grid point N: \( T_{\mathrm{N}} = h^{2}/12 W_{\mathrm{N}} \)

Calculates the U-matrix from T-matrix at grid point N: \(U_{\mathrm{N}} = 12(\mathbf{I} - 10 T_{\mathrm{N}})^{-1} - 10 \mathbf{I}\)

forces lowercase on a single character

Checks if the elastic_xs_threshold (threshold for elastic XS) and inelastic_xs_threshold (threshold for inelastic XS) conditions are already fulfilled.

Check variables read from namelist "basis"

Check variables read from namelist "input"

Check variables read from namelist "potential"

checks the number of non-zero PES matrix elements due to the \bar{\Omega} = \bar{\Omega}' condition, "number_of_nonzero_pes_matrix_elements", and the total number of non-zero algebraic coefficients, \( g_{{\lambda},\gamma,\gamma'}^{Jp} \), in the whole matrix, "number_of_nonzero_algebraic_coefficients".

check if the number of channels is the same in body-fixed and space-fixed frames

checks if the calculated sum of squares equals 1 for each channel

Computes the imaginary part of the cross-section contribution for given indices.

Calculates cross-section for a given initial and final state.

Computes the real part of the cross-section contribution for given indices.

counts the energetically accessible levels in the basis

counts the energetically accessible channels in the given block

Checks if the input value equals 0; used in the calculation of off-diagonal elements of the centrifugal matrix; see Eq. 5 in "Coupling Matrix" section

Determine the largest partial elastic and inelastic cross-sections in a given set of cross-sections.

utility function used by MSTA1 and MSTA2.

check the status during various io operations on files

Locates the correct quantum number that describes the v/j coupling. Note that coupling terms are symmetric with respect to the change of pre- and post-collisional quantum numbers.

Locates given \(\lambda\) value in legendre_indices.

Finds and returns the reduced term for the given indices.

Converts a floating-point number to a character string.

Returns indices from channel_indices_ that match the specified quantum state.

Returns the interpolated value of a specific radial coupling term at a given distance.

Handles error when the appropriate coupling term is not found.

Handles propagation at the last two grid points: R_{N-1} and R_{N}: provides T-matrix at N-1, N and N+1 points and the Ratio matrix at N and N+1 points

Handles error when \(\lambda\) is not found in legendre_indices.

handle printing messages depending on the outcome of unitarity check

interface for the following message: incorrect value encountered: variable_name = variable_value

incorrect value encountered: variable_name = variable_value

incorrect value encountered: variable_name = variable_value

incorrect value encountered: variable_name = variable_value

incorrect value encountered: variable_name = variable_value

Initial setup of the propagator: call centrifugal matrix (kept throughotu the propagation) and other matrices at \(R_{\mathrm{min}}\)

allocate arrays keeping accumulated and partial cross-sections in each jtot and parity block

launches "check_nonzero_pes_matrix_elements" and "prepare_pes_matrix_elements" subroutines; called from the main program

summarize the input parameters for the current run

transfers integer to a character

Interpolates the radial coupling terms using cubic spline functions. The resulting spline coefficients for each coupling term are stored in coupling_terms_b_coeffs, coupling_terms_c_coeffs, and coupling_terms_d_coeffs matrices.

checks if a channel/level is energetically accessible (open) by comparing energy with total_energy

checks if the sum of 3 integers is an even integer

returns interpolated value at guven u_ point number of points and ascending order of x is not checked since ispline is called after "spline" where these checks are done

calculates the ratio of the modified Bessel function of the second kind K_{l_ + 1/2}(x) and its first derivative (Eq. 8 in the "Solution of the coupled equations" section) Uses Temme's algorithm [N. M. Temme, J. Comput. Phys. 19 (1975) 324], implemented in "modified_bessel_temme_algorithm" subroutine; Unfortunately, the "ikv" function from special_functions library failed at large x_ values.

Implementation of the Temme's algorithm [N. M. Temme, J. Comput. Phys. 19 (1975) 324] to calculating modified Bessel functions of the second kind. This is a direct modernization of the "mbessk" subroutine in MOLSCAT: https://github.com/molscat/molscat/blob/master/source_code/rbessk.f

determines a backward recurrence starting point for Jn(x).

determines a backward recurrence starting point for Jn(x).

print a short message that there are no open channels in given block

Renormalized Numerov propagator

calculates the P coefficients from Eq. (3) in "Solution of coupled equations"

calculates Percival-Seaton coefficients (body-fixed variant) $$\begin{equation} \label{eq:algebraic_coeffs} g_{{\lambda},\gamma,\gamma'}^{Jp} = \delta_{\bar{\Omega},\bar{\Omega}'} (-1)^{\bar{\Omega}} \sqrt{(2j+1)(2j'+1)} \begin{pmatrix} j & j' & \lambda \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} j & j' & \lambda \\ \bar{\Omega} & -\bar{\Omega} & 0 \end{pmatrix}. \end{equation}$$

-- nonzero_terms_per_element - number of non-vanishing terms in the sum over \(\lambda\) in Eq. 1 in the "Coupling Matrix" section -- nonzero_legendre_indices - corresponding \(\lambda\) value for each non-vanishing coefficient is saved as an index to "legendre_indices" -- nonzero_algebraic_coefficients -- holds all non-vanishing \( g_{{\lambda},\gamma,\gamma'}^{Jp} \) coefficients

Prepare an array of wavevectors in a given block (in A^2) which are saved in the S-matrix file

Prints information about cross-sections from provided "cross_sections_" array

Prints basic information about the largest elastic and inelastic state-to-state xs (print_level <= 2)

prints information about body-fixed channels on screen

Prints information about cross-sections at the end of each total angular momentum (jtot) block

Prints detailed information about the largest elastic and inelastic state-to-state xs (print_level >= 3)

Prints information about cross-sections at the end of the program

print the final warning that the unitarity check failed in given blocks

Print the largest partial elastic and inelastic state-to-state cross-sections in a given block.

print a shor summary on the number of non-zero matrix elements of the PES matrix

Prints quantum numbers describing radial coupling terms of the PES based on the provided set type and column count.

print S-matrix on screen

calculates the non-zero algebraic coefficients \( g_{{\lambda},\gamma,\gamma'}^{Jp} \) for a single matrix element - see Eq. (1) in the "Coupling matrix" section; algebraic coefficients are saved to nonzero_algebraic_coefficients array; corresponding indices to legendre_indices are saved to nonzero_legendre_indices array

Print a simple message after the propagation is finished

computes Riccati-Bessel function of the first kind, and derivatives.

computes Riccati-Bessel function of the second kind, and derivatives.

Reads the value of lambda and compares with expected value.

reads the input file prepared by the user using NAMELIST feature the code uses 3 namelists: "input", "basis" and "potential"

Reads the intermolecular distance and radial coupling terms formatted in columns by iterating over number of tabulated ]](R\) points. Immediately converts \(R\) and radial coupling terms to a.u.

Reads the radial coupling terms from the external file. The file is assumed to be formatted as described in "Supplying radial terms" section. The read radial coupling terms are kept in "tabulated_coupling_terms"

Reduces the coupling terms based on the existence of couplings.

Reduces the tabulated_coupling_terms matrix to retain only the necessary coupling terms.

Calculates 1/Gamma(1-X); modernized version of Molscat's rgamma function; see: https://github.com/molscat/molscat/blob/36fa8f93a92f851e9d84245dd6a972e2910541c5/source_code/rbesjy.f --------------------------------------------------------------------!

calculates the Riccati-Bessel function of the first kind and its first derivative. Calls the rctj function from special_functions.f90

calculates the Riccati-Bessel function of the second kind and its first derivative. Calls the rcty function from special_functions.f90

saves indices to open levels in the basis and corresponding wavevectors (in A^2)

save "header" of the partial cross-sections file: -- label, "itype", number of levels in the basis, reduced mass of the system -- vibrational and rotational quantum numbers -- rovibrational energies -- index pointing to the initial level and the kinetic/total energy

save information about current block -- total angular momentum, parity exponent, number of open channels in the current block -- array of indices pointing to the basis arrays, array holding \(l\) values, wavenumbers -- real part of the S-matrix -- imaginary part of the S-matrix

save "header" of the S-matrix file: -- label, "itype", number of levels in the basis, reduced mass of the system -- vibrational and rotational quantum numbers -- rovibrational energies -- index pointing to the initial level and the kinetic/total energy

Prepares the channel_indices array which holds indices that refer to the basis arrays: v1level/j1level/internal_energies, and channels_omega_values which holds values of \bar{\Omega}

determine the number of scattering channels in each parity block for given total angular momentum in both body-fixed and space-fixed frames

Prepares the channel_l_values array which holds values of orbital angular momentum, \(l\), a space-fixed-frame quantum number.

Skips the first n_skip_lines (read on input) lines in the coupling_terms_file

determines b, c and d coefficients of the cubic spline function y(x) = y_i + b_i * dx + c_i * dx^2 + d_i * dx^3, where dx = x - x_i, and x_i <= x < x_i+1. The algorithm is based on Gerald, C., and Wheatley, P., "Applied Numerical Analysis", Addison-Wesley, 1994.

Sum the contributions to the cross-section from S-matrix components

print the message about the time it took to complete a single task

forces lowercase on given string

returns the total energy

check if the triangle inequality for 3 variables hols

checks the unitarity of the S-matrix (Eq. (13) in "Solution of coupled equations")

converts all physical quantities to atomic units

update channel_indices array which holds indices within the loop over level_index_ in set_body_fixed_channels

updates number_of_channels_even and number_of_channels_odd in the body-fixed frame for given \(J\) and \(\bar{Omega}_{max}\)

updates number_of_channels_even and number_of_channels_odd in the space-fixed frame for given range of orbital angular momentum

Checks if read R values are consistent with r_min and r_max.

returns the squared wavevector, \(k_{a}^{2}\), given the energy of a given state, \(E_{a}\); calls etot() function; atomic units in the whole function \( k_{a} = \sqrt(2 \mu (E_{tot} - E_{a}) \) since it uses reduced_mass and total_energy(), the function checks if units are already converted

writes an error message on a chosen unit

writes headers on screen

writes a message on a chosen unit

writes a warning message on a chosen unit

checks the condition for nonvanishing 3-j symbol with zero projections: triangle inequality on x,y,z and if the sum x+y+z is an even integer