SCATTERING Program

program SCATTERING
   !---------------------------------------------------------------------------!
   !!                      the SCATTERING code - version 0.0.1
   !!                    optimized for diatom-atom calculations
   !!                                 25-01-2024
   !!                                 H. Jozwiak
   !---------------------------------------------------------------------------!
   !! Please, refer to this version of the code by citing the following paper
   !! H. Jozwiak, F. Thibault, A. Viel, P. Wcislo, F. Lique,
   !! Rovibrational (de-)excitation of H2 by He revisited
   !! https://doi.org/10.48550/arXiv.2311.09890
   !---------------------------------------------------------------------------!
   use, intrinsic :: iso_fortran_env, only: int32, sp => real32, dp => real64
   use utility_functions_mod, only: write_header,                              &
      write_message, float_to_character, integer_to_character,                 &
      time_count_summary, no_open_channels_message
   use array_operations_mod, only: append, allocate_1d, allocate_2d
   use global_variables_mod
   use physics_utilities_mod, only: count_open_basis_levels,                   &
      save_open_basis_levels, units_conversion
   use input_reader_mod, only: read_input_file
   use radial_coupling_terms_mod, only: read_radial_coupling_terms,            &
      reduce_radial_coupling_terms, interpolate_radial_coupling_terms
   use channels_mod, only: set_number_of_channels, set_body_fixed_channels,    &
      set_space_fixed_channels, count_open_channels_in_block,                  &
      calculate_largest_wavevector, calculate_number_of_steps,                 &
      prepare_wavevector_array, print_short_block_summary, print_channels
   use pes_matrix_mod, only: initialize_pes_matrix
   use propagator_mod, only: numerov
   use boundary_conditions_mod, only: calculate_sf_matrix_from_bf_matrix,      &
      calculate_k_matrix, calculate_s_matrix
   use unitarity_check_mod, only: unitarity_check, print_final_unitarity_warning
   use state_to_state_cross_sections_mod, only:initialize_cross_section_arrays,&
      calculate_state_to_state_cross_section,  add_cross_sections,             &
      print_largest_partial_cross_sections, print_cross_sections_for_jtot,     &
      print_final_cross_sections,  determine_largest_cross_sections,           &
      check_cross_section_thresholds, save_partial_xs_file_header,             &
      save_partial_xs_single_block
   use save_s_matrix_mod, only: save_s_matrix_file_header, save_s_matrix_block_info
   !---------------------------------------------------------------------------!
   implicit none
   !---------------------------------------------------------------------------!
   integer(int32) :: number_of_channels, size_even, size_odd,                  &
      number_of_open_basis_levels, jtot_, parity_exponent, number_of_steps,    &
      number_of_open_channels, max_elastic_index, max_inelastic_index_1,       &
      max_inelastic_index_2
   !---------------------------------------------------------------------------!
   integer(int32) :: count_blocks = 0
   integer(int32) :: consecutive_blocks_thresholddiag = 0
   integer(int32) :: consecutive_blocks_thresholdoff = 0
   !---------------------------------------------------------------------------!
   real(dp) :: largest_wavevector, k_potential_depth,                          &
      max_elastic_cross_section, max_inelastic_cross_section, time_total_start,&
      time_total_stop, time_total, time_init_stop, time_init, time_jtot_start, &
      time_jtot_stop, time_jtot, time_parity_start,time_parity_stop,           &
      time_parity
   logical :: unitarity_block_check
   logical :: terminate = .false.
   integer, allocatable :: channel_indices(:), channels_omega_values(:),&
      channel_l_values(:), open_basis_levels(:), nonzero_terms_per_element(:), &
      nonzero_legendre_indices(:), smatcheckarr(:)
   real(dp), allocatable :: basis_wavevectors(:), block_wavevectors(:),        &
      nonzero_algebraic_coefficients(:), accumulated_cross_sections(:),        &
      partial_cross_sections_block(:), partial_cross_sections_jtot(:)
   real(dp), allocatable :: BF_log_der_matrix(:,:), SF_log_der_matrix(:,:),    &
      k_matrix(:,:), s_matrix_real(:,:), s_matrix_imag(:,:)
   !---------------------------------------------------------------------------!
   !---------------------------------------------------------------------------!
   ! Initizalization: start the time count
   !---------------------------------------------------------------------------!
   call cpu_time(time_total_start)
   !---------------------------------------------------------------------------!
   ! Initialize fwigxjpf library
   !---------------------------------------------------------------------------!
   call fwig_table_init(2*100, 9)
   call fwig_temp_init(2*100)
   !---------------------------------------------------------------------------!
   ! Print the header
   !---------------------------------------------------------------------------!
   call write_header("main")
   !---------------------------------------------------------------------------!
   ! Read the input file
   !---------------------------------------------------------------------------!
   call read_input_file
   !---------------------------------------------------------------------------!
   ! S-matrix file: write input parameters and basis levels
   !---------------------------------------------------------------------------!
   call save_s_matrix_file_header
   !---------------------------------------------------------------------------!
   ! Prepare the file with the partial XS
   !---------------------------------------------------------------------------!
   if (print_partial_cross_sections) then
      call save_partial_xs_file_header
   endif
   !---------------------------------------------------------------------------!
   ! Convert units: starting now, everything is in atomic units
   !---------------------------------------------------------------------------!
   call units_conversion
   !---------------------------------------------------------------------------!
   ! Read the radial terms of the potential from external file
   !---------------------------------------------------------------------------!
   call read_radial_coupling_terms
   !---------------------------------------------------------------------------!
   ! Reduce matrix elements that are not needed
   !---------------------------------------------------------------------------!
   call reduce_radial_coupling_terms
   !---------------------------------------------------------------------------!
   ! Interpolate radial terms                                            
   !---------------------------------------------------------------------------!
   call interpolate_radial_coupling_terms
   !---------------------------------------------------------------------------!
   ! Search for energetically accessible levels and prepare the arrays that are
   ! needed in the calculations of the state-to-state XS                       
   !---------------------------------------------------------------------------!
   number_of_open_basis_levels = count_open_basis_levels()
   call save_open_basis_levels(number_of_open_basis_levels, open_basis_levels, &
      basis_wavevectors)
   !---------------------------------------------------------------------------!
   ! Initialize arrays that save the state-to-state cross-sections
   !---------------------------------------------------------------------------!
   call initialize_cross_section_arrays(number_of_open_basis_levels,           &
      partial_cross_sections_block, partial_cross_sections_jtot,               &
      accumulated_cross_sections)
   !---------------------------------------------------------------------------!
   ! Initialization is finished                                                
   !---------------------------------------------------------------------------!
   call cpu_time(time_init_stop)
   if (print_level.ge.2) call time_count_summary(time_total_start,             &
      time_init_stop, time_init, "Initialization completed in ")
   !---------------------------------------------------------------------------!
   ! Loop over total angular momentum
   !---------------------------------------------------------------------------!
   call write_header("jtot_loop")
   !---------------------------------------------------------------------------!
   do jtot_ = jtot_min,jtot_max,jtot_step
      !------------------------------------------------------------------------!
      call write_header("block", opt_integer_ = jtot_)
      !------------------------------------------------------------------------!
      call cpu_time(time_jtot_start)
      !------------------------------------------------------------------------!
      partial_cross_sections_jtot = 0
      call set_number_of_channels(jtot_, size_even, size_odd)
      !------------------------------------------------------------------------!
      do parity_exponent = 0,1
         !---------------------------------------------------------------------!
         call cpu_time(time_parity_start)
         !---------------------------------------------------------------------!
         select case(parity_exponent)
            case(0)
               number_of_channels = size_even
            case(1)
               number_of_channels = size_odd
         end select
         !---------------------------------------------------------------------!
         if (number_of_channels == 0) cycle
         !---------------------------------------------------------------------!
         ! Summary of the current block
         !---------------------------------------------------------------------!
         count_blocks = count_blocks+1
         if (print_level.ge.1) then
            call print_short_block_summary(jtot_, parity_exponent,             &
               count_blocks, number_of_channels)
         endif
         !---------------------------------------------------------------------!
         ! Prepare of the basis for each J/p block:
         ! channels_omega_values holds all values of omega (BF_)
         ! channel_l_values holds all values of l (SF_)
         ! channel_indices holds the indices which refer to the basis arrays:
         !   --   v1level/j1level/internal_energies
         !---------------------------------------------------------------------!
         call allocate_1d(channels_omega_values,number_of_channels)
         call allocate_1d(channel_l_values,number_of_channels)
         call allocate_1d(channel_indices,number_of_channels)
         !---------------------------------------------------------------------!
         ! Prepare channels_omega_values, channel_indices and channel_l_values
         !---------------------------------------------------------------------!
         call set_body_fixed_channels(jtot_, parity_exponent, channel_indices, &
            channels_omega_values)
         call set_space_fixed_channels(jtot_, parity_exponent, channel_l_values)
         !---------------------------------------------------------------------!
         ! Print the BF quantum numbers on screen
         !---------------------------------------------------------------------!
         if (print_level.ge.1) call print_channels(parity_exponent,            &
            channel_indices, channels_omega_values)
         !---------------------------------------------------------------------!
         ! Determine the number of open (energetically accessible) channels
         !---------------------------------------------------------------------!
         number_of_open_channels = count_open_channels_in_block(channel_indices)
         !---------------------------------------------------------------------!
         ! If there are no open channels, skip this block
         !---------------------------------------------------------------------!
         if (number_of_open_channels == 0) then
            call no_open_channels_message(count_blocks)
            cycle
         endif
         !---------------------------------------------------------------------!
         ! Array of wavevectors in given block (to save in the S-matrix file)
         !---------------------------------------------------------------------!
         call allocate_1d(block_wavevectors, number_of_open_channels)
         call prepare_wavevector_array(channel_indices, block_wavevectors)
         !---------------------------------------------------------------------!
         ! Determine the largest wavevector in the block
         !---------------------------------------------------------------------!
         largest_wavevector = calculate_largest_wavevector(channel_indices)
         !---------------------------------------------------------------------!
         ! Determine the number of steps on the intermolecular (R) grid
         !---------------------------------------------------------------------!
         k_potential_depth = dsqrt(2*reduced_mass*potential_depth)
         number_of_steps = calculate_number_of_steps(largest_wavevector,       &
            k_potential_depth)
         !---------------------------------------------------------------------!
         ! Prepare the PES matrix
         !---------------------------------------------------------------------!
         call initialize_pes_matrix(channel_indices, channels_omega_values,    &
            nonzero_terms_per_element, nonzero_legendre_indices,               &
            nonzero_algebraic_coefficients)
         !---------------------------------------------------------------------!
         ! Allocate asymptotic body-fixed (BF) log-derivative matrix
         ! and call the propagator
         !---------------------------------------------------------------------!
         call allocate_2d(BF_log_der_matrix, number_of_channels,               &
            number_of_channels)
         call numerov(number_of_channels, channel_indices,                     &
            channels_omega_values, nonzero_terms_per_element,                  &
            nonzero_legendre_indices, nonzero_algebraic_coefficients,          &
            number_of_steps, jtot_, BF_log_der_matrix)
         !---------------------------------------------------------------------!
         ! Allocate asymptotic space-fixed (SF) log-derivative matrix and
         ! transform the BF result to the SF frame
         !---------------------------------------------------------------------!
         call allocate_2d(SF_log_der_matrix, number_of_channels, number_of_channels)
         call calculate_sf_matrix_from_bf_matrix(number_of_channels, jtot_,    &
            channel_indices, channels_omega_values, channel_l_values,          &
            BF_log_der_matrix, SF_log_der_matrix)
         !---------------------------------------------------------------------!
         ! Get the K-matrix from log-derivative matrix
         ! (Eq. 4 in "Solution of the coupled equations")
         !---------------------------------------------------------------------!
         call allocate_2d(k_matrix, number_of_open_channels, number_of_open_channels)
         call calculate_k_matrix(number_of_channels, SF_log_der_matrix,        &
            number_of_open_channels, channel_indices, channel_l_values,        &
            r_max, k_matrix)
         !---------------------------------------------------------------------!
         ! Get the S-matrix from the K-matrix
         ! (Eq. 12 in "Solution of the coupled equations")
         !---------------------------------------------------------------------!
         call allocate_2d(s_matrix_real,number_of_open_channels,number_of_open_channels)
         call allocate_2d(s_matrix_imag,number_of_open_channels,number_of_open_channels)
         call calculate_s_matrix(number_of_open_channels,k_matrix,s_matrix_real,s_matrix_imag)
         !---------------------------------------------------------------------!
         ! S-matrix is written to the binary S-matrix file
         !---------------------------------------------------------------------!
         call save_s_matrix_block_info(jtot_, parity_exponent,                 &
            number_of_open_channels, channel_indices, channel_l_values,        &
            block_wavevectors, s_matrix_real, s_matrix_imag)
         !---------------------------------------------------------------------!
         ! Check if the S-matrices are unitary
         !---------------------------------------------------------------------!
         call unitarity_check(number_of_open_channels, s_matrix_real,          &
            s_matrix_imag, unitarity_block_check)
         !---------------------------------------------------------------------!
         ! If the unitary is not fulfilled, keep the information about this block
         !---------------------------------------------------------------------!
         if (.not.(unitarity_block_check)) then
            call append(smatcheckarr, jtot_)
         endif
         !---------------------------------------------------------------------!
         ! Calculate all available cross-sections
         !---------------------------------------------------------------------!
         call calculate_state_to_state_cross_section(jtot_, open_basis_levels, &
            basis_wavevectors,s_matrix_real,s_matrix_imag,channel_indices,     &
            channel_l_values,partial_cross_sections_block)
         !---------------------------------------------------------------------!
         ! Print the results from this parity block to the partial XS file
         ! and add the calculated partial XS to the partial_cross_sections_jtot array
         !---------------------------------------------------------------------!
         if (print_partial_cross_sections) then
            call save_partial_xs_single_block(jtot_, count_blocks,             &
               open_basis_levels, partial_cross_sections_block)
         endif
         call add_cross_sections(number_of_open_basis_levels,                  &
            partial_cross_sections_block, partial_cross_sections_jtot)
         !---------------------------------------------------------------------!
         ! Check the time after each parity block:
         !---------------------------------------------------------------------!
         call cpu_time(time_parity_stop)
         if (print_level.ge.2) call time_count_summary(time_parity_start,      &
            time_parity_stop, time_parity, "Parity block completed in ")
         !---------------------------------------------------------------------!
         ! ... end of the loop over parity                                      
         !---------------------------------------------------------------------!
         call write_message(repeat("-", 90))
      enddo
      !------------------------------------------------------------------------!
      ! Add the cross-sections from this Jtot block:
      !------------------------------------------------------------------------!
      call add_cross_sections(number_of_open_basis_levels,                     &
            partial_cross_sections_jtot, accumulated_cross_sections)
      !------------------------------------------------------------------------!
      ! Determine the largest partial elastic/inelastic XS in this Jtot block:
      !------------------------------------------------------------------------!
      call determine_largest_cross_sections(open_basis_levels,                 &
         partial_cross_sections_jtot, max_elastic_cross_section,               &
         max_inelastic_cross_section, max_elastic_index,                       &
         max_inelastic_index_1, max_inelastic_index_2)
      !------------------------------------------------------------------------!
      call print_largest_partial_cross_sections(jtot_,                         &
         max_elastic_cross_section, max_inelastic_cross_section,               &
         max_elastic_index, max_inelastic_index_1, max_inelastic_index_2,      &
         open_basis_levels)
      !------------------------------------------------------------------------!
      if (jtot_max == 999999) then
         call check_cross_section_thresholds(max_elastic_cross_section,        &
            max_inelastic_cross_section, consecutive_blocks_thresholddiag,     &
            consecutive_blocks_thresholdoff, terminate)
      endif
      !------------------------------------------------------------------------!
      ! Check the time after each JTOT block:                                  
      !------------------------------------------------------------------------!
      call cpu_time(time_jtot_stop)
      !------------------------------------------------------------------------!
      ! Print all the XS after current JTOT block                              
      !------------------------------------------------------------------------!
      if (print_level.ge.3) then
         call print_cross_sections_for_jtot(jtot_, open_basis_levels,          &
            accumulated_cross_sections)
      endif
      !------------------------------------------------------------------------!
      if (print_level.ge.2) call time_count_summary(time_jtot_start,           &
         time_jtot_stop, time_jtot, "Total angular momentum block completed in ")
      !------------------------------------------------------------------------!
      ! terminate the loop if elastic_xs_threshold/inelastic_xs_threshold
      ! condition is satisfied
      !------------------------------------------------------------------------!
      if (terminate) exit
   enddo
   !---------------------------------------------------------------------------!
   call write_header("loop_terminated")
   call write_header("summary")
   !---------------------------------------------------------------------------!
   ! if for some JTOTs the S-matrix did not fulfill the unitary check,
   ! these are listed here
   !---------------------------------------------------------------------------!
   if (allocated(smatcheckarr)) then
      call print_final_unitarity_warning(smatcheckarr)
   endif
   !---------------------------------------------------------------------------!
   ! Print all the calculated XS                                                
   !---------------------------------------------------------------------------!
   call print_final_cross_sections(open_basis_levels, accumulated_cross_sections)
   !---------------------------------------------------------------------------!
   call fwig_temp_free();
   call fwig_table_free();
   !---------------------------------------------------------------------------!
   ! Stop the time count
   !---------------------------------------------------------------------------!
   call cpu_time(time_total_stop)
   call time_count_summary(time_total_start, time_total_stop, time_total,      &
      "Total CPU time: ")
   !---------------------------------------------------------------------------!
   close(s_matrix_unit)
   !---------------------------------------------------------------------------!
   if (print_partial_cross_sections) then
      close(partial_file_unit)
   endif
   !---------------------------------------------------------------------------!
end program SCATTERING