print_thermochemistry Subroutine

   subroutine print_thermochemistry(result, electronic_energy, unit)
      !! Print thermochemistry results.
      type(thermochemistry_result_t), intent(in) :: result
      real(dp), intent(in) :: electronic_energy  !! Electronic energy in Hartree
      integer, intent(in), optional :: unit      !! Output unit (default: stdout)

      real(dp) :: H_total_cal, Cv_total, S_total, S_total_J
      real(dp) :: H_vib_cal, H_rot_cal, H_trans_cal
      real(dp) :: H_T, TS, G_T
      real(dp) :: total_energy, total_enthalpy, total_free_energy
      real(dp) :: H0_HT_PV, pV_cal
      character(len=512) :: line

      ! pV term = RT for ideal gas (in cal/mol)
      pV_cal = R_CALMOLK*result%temperature

      ! Compute enthalpy contributions (what xtb shows in "enthalpy" column)
      ! VIB: thermal vibrational energy only (ZPE is reported separately)
      H_vib_cal = result%E_vib*HARTREE_TO_CALMOL
      ! ROT: just internal energy (3/2)RT for nonlinear, RT for linear - no pV term
      H_rot_cal = result%E_rot*HARTREE_TO_CALMOL
      ! TR: includes pV term, so (5/2)RT for translation enthalpy
      H_trans_cal = result%E_trans*HARTREE_TO_CALMOL + pV_cal
      ! TOT: sum of thermal contributions (VIB + ROT + TR), NOT including ZPE
      H_total_cal = H_vib_cal + H_rot_cal + H_trans_cal

      ! Use Cp for translation (Cv + R) as is standard for thermochemistry output
      Cv_total = (result%Cv_trans + R_CALMOLK) + result%Cv_rot + result%Cv_vib
      S_total = result%S_trans + result%S_rot + result%S_vib + result%S_elec
      S_total_J = S_total*CAL_TO_J  ! cal to J

      ! Thermodynamic quantities in Hartree
      ! H(0)-H(T)+PV is the thermal correction WITHOUT ZPE (just E_trans + E_rot + E_vib + RT)
      H0_HT_PV = result%E_trans + result%E_rot + result%E_vib + R_HARTREE*result%temperature
      ! H(T) is full thermal correction including ZPE
      H_T = result%thermal_correction_enthalpy
      TS = result%thermal_correction_enthalpy - result%thermal_correction_gibbs
      G_T = result%thermal_correction_gibbs

      total_energy = electronic_energy
      total_enthalpy = electronic_energy + result%thermal_correction_enthalpy
      total_free_energy = electronic_energy + result%thermal_correction_gibbs

      ! Print header
      call logger%info(" ")
      call logger%info("Thermochemistry (RRHO)")
      call logger%info("======================")
      call logger%info(" ")

      ! Setup section - simple list
      write (line, '(A,F10.4,A)') "  Temperature:       ", result%temperature, " K"
      call logger%info(trim(line))
      write (line, '(A,F10.4,A)') "  Pressure:          ", result%pressure, " atm"
      call logger%info(trim(line))
      write (line, '(A,F10.4,A)') "  Molecular mass:    ", result%total_mass, " amu"
      call logger%info(trim(line))
      write (line, '(A,I6)') "  Vibrational modes: ", result%n_real_freqs
      call logger%info(trim(line))
      if (result%n_imag_freqs > 0) then
         write (line, '(A,I6,A)') "  Imaginary freqs:   ", result%n_imag_freqs, " (skipped)"
         call logger%info(trim(line))
      end if
      if (result%is_linear) then
         call logger%info("  Linear molecule:   yes")
      else
         call logger%info("  Linear molecule:   no")
      end if
      write (line, '(A,I6)') "  Symmetry number:   ", result%symmetry_number
      call logger%info(trim(line))
      call logger%info(" ")

      ! Contribution table
      call logger%info("  temp (K)       q        H(cal/mol)  Cp(cal/K/mol)  S(cal/K/mol)  S(J/K/mol)")
      call logger%info("  -------------------------------------------------------------------------")

      write (line, '(F8.2,A,ES10.3,F12.3,F14.3,F14.3,F12.3)') &
         result%temperature, "  VIB", result%q_vib, H_vib_cal, result%Cv_vib, &
         result%S_vib, result%S_vib*CAL_TO_J
      call logger%info(trim(line))

      write (line, '(A,ES10.3,F12.3,F14.3,F14.3,F12.3)') &
         "          ROT", result%q_rot, H_rot_cal, result%Cv_rot, &
         result%S_rot, result%S_rot*CAL_TO_J
      call logger%info(trim(line))

      write (line, '(A,ES10.3,F12.3,F14.3,F14.3,F12.3)') &
         "          INT", result%q_rot*result%q_vib, H_vib_cal + H_rot_cal, &
         result%Cv_vib + result%Cv_rot, result%S_vib + result%S_rot, &
         (result%S_vib + result%S_rot)*CAL_TO_J
      call logger%info(trim(line))

      ! For TR, report Cp = Cv + R (constant pressure heat capacity for ideal gas)
      write (line, '(A,ES10.3,F12.3,F14.3,F14.3,F12.3)') &
         "          TR ", result%q_trans, H_trans_cal, result%Cv_trans + R_CALMOLK, &
         result%S_trans, result%S_trans*CAL_TO_J
      call logger%info(trim(line))

      call logger%info("  -------------------------------------------------------------------------")
      write (line, '(A,F12.3,F14.3,F14.3,F12.3)') &
         "          TOT           ", H_total_cal, Cv_total, S_total, S_total_J
      call logger%info(trim(line))

      call logger%info(" ")

      ! Thermal corrections table
      call logger%info(" ")
      call logger%info("Thermal Corrections (Hartree)")
      call logger%info("-----------------------------")
      write (line, '(A,F18.12)') "  Zero-point energy:              ", result%zpe_hartree
      call logger%info(trim(line))
      write (line, '(A,F18.12)') "  Thermal correction to Energy:   ", result%thermal_correction_energy
      call logger%info(trim(line))
      write (line, '(A,F18.12)') "  Thermal correction to Enthalpy: ", result%thermal_correction_enthalpy
      call logger%info(trim(line))
      write (line, '(A,F18.12)') "  Thermal correction to Gibbs:    ", result%thermal_correction_gibbs
      call logger%info(trim(line))
      call logger%info(" ")

      ! Final totals
      call logger%info("Total Energies (Hartree)")
      call logger%info("------------------------")
      write (line, '(A,F20.12)') "  Electronic energy:            ", total_energy
      call logger%info(trim(line))
      write (line, '(A,F20.12)') "  Electronic + ZPE:             ", total_energy + result%zpe_hartree
      call logger%info(trim(line))
      write (line, '(A,F20.12)') "  Electronic + thermal (E):     ", total_energy + result%thermal_correction_energy
      call logger%info(trim(line))
      write (line, '(A,F20.12)') "  Electronic + thermal (H):     ", total_enthalpy
      call logger%info(trim(line))
      write (line, '(A,F20.12)') "  Electronic + thermal (G):     ", total_free_energy
      call logger%info(trim(line))
      call logger%info(" ")

   end subroutine print_thermochemistry