Using HITEMP line lists

A guide on how to download HITEMP line lists can be found in the downloading line list data tutorial. To download HITEMP data, you need to be logged in on the HITRAN website. Therefore, we will download the demonstration data from a cloud, just to show how the calculation works. If you want to use the HITEMP data for real, please download the line list files from the HITRAN website.

We will calculate the HITEMP CO line list here. Credit for the data goes to G. Li et al. 2015. First we need to import the relevant packages and set up the folder structure.

import numpy as np
import matplotlib.pyplot as plt
import line_racer.line_racer as lr
import urllib.request
import zarr.storage
import os

# if you want to use the conversion to the pRT format
import petitRADTRANS

lr.LineRacer.check_installation()

os.makedirs('line_list/CO', exist_ok=True)

Running line racer installation test...
line racer installation test passed successfully.

After that, we are downloading the line list data. As mentioned, it is downloaded from a cloud using urllib.request. If that is not working for you, please download it manually from HITRAN or here. However, the partition function can be downloaded from the HITRAN website directly. For HITRAN data the .par file is needed, which contains all line information. Additionally, the partition function data is needed, which is stored in a .txt file.

urllib.request.urlretrieve('https://keeper.mpdl.mpg.de/f/be281779a65542db814e/?dl=1', 'line_list/CO/05_HITEMP2019.par')
urllib.request.urlretrieve('https://hitran.org/data/Q/q26.txt', 'line_list/CO/q26.txt')

('line_list/CO/q26.txt', <http.client.HTTPMessage at 0x13713b790>)

As for ExoMol, we need to define the pressure and temperature grid. Here, we are only using one pressure and temperature point for demonstration purposes. If you want to use a pre-defined petitRADTRANS grid, you can use the function lr.LineRacer.prt_pressure_temperature_grid() as shown in the ExoMol example.

pressures = [1.0]  # in bar
temperatures = [1500.0]  # in K

After that, we can define the line racer object, which already contains most of the relevant information for our opacity calculation. In general, everything is written in lower case letters, except for the line list name, since it is only used for naming the output files.

CO_racer = lr.LineRacer(resolution=1e6,
                        cutoff=100.0,  # in 1/cm
                        lambda_min=1.1e-5,  # in cm
                        lambda_max=2.5e-2,  # in cm
                        grid_type='log',
                        sublorentzian='Hartmann',
                        database='hitemp',
                        input_folder='line_list/CO/',  # path to folder with input files
                        temperatures=temperatures,  # in K
                        pressures=pressures,  # in bar
                        species_isotope_dict={'12C-16O': 1.0},
                        line_list_name='HITEMP2019',  # you can specify the version here
                        broadening_species_dict={'air': 1.0},
                        broadening_type='hitran_table',  # could also be 'sharp_burrows', 'hitran_table' or 'constant'

                        # Change the following only if you know what you do!
                        force_molliere2015_below_wn=40,  # use the Mollière et al. 2015 method for lines with effective wavenumbers smaller than this threshold.
                        force_molliere2015_method=False  # whether to force using Mollière et al. (2015) for line profile calculation, recommended only if warning appears or for small line lists
                        )

A detailed description of the class can be found in the ExoMol section. Here, we will just explain the differences.

The input folder where the line list files (.par/.out, .txt) are stored is now the one containing the HITEMP files. The naming convention of the dictionary with the isotopes to calculate is still in the ExoMol naming convention and will be translated internally. Since the HITEMP line lists do not have specific names, we just call it ‘HITEMP2019’. For the dictionary with the broadening species air and self is possible for .par files. The broadening species is air and the broadening type HITRAN table, since HITRAN is the provider of the HITEMP line lists. The molecular mass is not required for HITEMP line lists, since it is read from the molparam.txt file automatically.

The explanation for force_molliere2015_methodand force_molliere2015_below_wn can be found in the ExoMol section.

HITEMP intensities are weighted by their natural abundance by default. In line racer, this is automatically corrected for and the opacities are scaled to abundance 1. If you want a different abundance, such as the natural abundance, you can set the desired abundance in the species_isotope_dict.

We also need to prepare the calculation by setting up the wavelength grid. Finally, we can provide the transition files list directly here, or let line racer search for all .par files in the input folder.

transition_files_list = CO_racer.prepare_opacity_calculation()

print("The transition file(s) to be used for the calculation are: ", transition_files_list)

The transition file(s) to be used for the calculation are:  ['line_list/CO/05_HITEMP2019.par']

After preparing the calculation, we can start the actual opacity calculation. Here, we can also choose whether to use MPI for parallelization or not. If you are on a cluster and have MPI set up, you can set use_mpi=True to use it. Otherwise, just set it to False and line racer will use multiprocessing for parallelization. In the case of normal multiprocessing (use_mpi=False) you can choose the number of cores by setting n_cores to the number of cores you want to use.

If you are using the line_racer_obj.calculate_opacity in a script make sure to put it in a if __name__ == "__main__": block to avoid problems with multiprocessing.

You can either directly use the petitRADTRANS format by setting prt_format=True or first calculate the cross-sections in line racer’s own format, which is a zarr compressed with zip. More about that, and which further information are needed for the petitRADTRANS format can be found in the ExoMol section. Here, they are set to None. We use the flag verbose=True, to get more information about the line calculation. In addition to the intensity correction files, HITEMP calculations also need the molparam.txt file, which is downloaded automatically if not present in the data/ folder.

# calculate the final cross-sections
final_cross_section_file_name = CO_racer.calculate_opacity(transition_files_list,
                                                           use_mpi=False,
                                                           n_cores=1,
                                                           sampling_boost=1.0,
                                                           coarse_grid_switch=True,

                                                           # pRT output options
                                                           prt_format=False,
                                                           doi=None,
                                                           additional_information=None,
                                                           use_prt_input_file_path=False,

                                                           # more output
                                                           verbose=True)

print("The cross-sections are stored in the file: ", final_cross_section_file_name)

Using cutoff and Hartmann intensity correction grid and interpolated to 100.0 1/cm
Reading HITRAN transition file...
Finished reading HITRAN transition file in 3.14 seconds.


 Line parameter calculation time: 0.010356903076171875 s
Starting the line profile calculation of 125496 lines

 Starting internal lines!
Line profile calculation done Internal lines done!
 Starting external lines!
Progress: 770/773 subgrids
 External lines done!

Total time for opacity calculation:  52.764594078063965 s

The cross-sections are stored in the file:  cross-sections/cross-section_12C-16O__HITEMP2019__40-90909cm-1.zarr.zip

The resulting cross-section file is stored in the folder cross-sections/ by default. The file contains the pressure grid in bar and temperature grid in Kelvin, the wavenumber grid in 1/cm and the cross-sections in \(\rm cm^2/molecule\) for each pressure and temperature point. If you want to convert the cross-sections to the petitRADTRANS format, you could directly set prt_format=True in the calculate_opacity function or just convert it afterward, as explained in the ExoMol section.

store = zarr.storage.ZipStore(final_cross_section_file_name, mode='a')
z = zarr.group(store=store)

# Access datasets
pressures = z['pressures'][:]
temperatures = z['temperatures'][:]
wavenumbers = z['wavenumbers'][:]
cross_section = z['cross-sections']

print("Pressures (bar): ", pressures)
print("Temperatures (K): ", temperatures)

# Extract slice and plot
cross_section_slice = cross_section[:]

Pressures (bar):  [1.]
Temperatures (K):  [1500.]

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(1/wavenumbers * 1e4, cross_section_slice[0, 0, :])  # convert wavenumber to wavelength in micron

ax.set_xlabel('Wavelength (micron)')
ax.set_ylabel('Cross-section (cm$^2$/molecule)')
ax.set_title(f'CO HITEMP at P={pressures[0]} bar, T={temperatures[0]} K')
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylim(1e-37, 1e-17)
ax.set_xlim(0.5, 10);