Repository: radis/radis (language: python)

Issue:
HITEMP-Hdf5 improvement
- renamed "hdf5" dbformat (a fileformat not giving any idea about the line database used) to "hitemp-radisdb"
- fix nonequilibrium calculations that were blocked on "hitemp-radisdb" for CO
(misc)
- add local files in Spectrum output
- improved error message when missing 'branch' (not fully parsed HITRAN file)
- add non-eq calculation with direct-HITEMP download [fails with 0.9.28]
- auto-download HITEMP CO database for noneq tests (14 Mb ~ 14 gr CO2 🌳 ).
- parse HITEMP local/global quanta before storing the fetched file to disk. Fixes previous test.

Add/vaex for hitemp
<!-- Please be sure to check out our contributing guidelines,
https://github.com/radis/radis/blob/develop/CONTRIBUTING.md . -->

### Description
<!-- Provide a general description of what your pull request does. -->

Adding vaex to load HDF5 files. Proof-of-concept : [test_vaex.ipynb](https://github.com/radis/radis-benchmark/blob/master/manual_benchmarks/test_vaex.ipynb)

Built on top of #281 : 

- [x] implemented different engines in HDF5 cache files & generation : vaex and h5py
- [x] added vaex filtering when loading database with wmin, wmax, or isotpes
- [ ] vaex is the default library to load HITEMP .hdf5 files  👉  "pytables" remains the default for the moment
- [x] auto update of your existing .hdF5 databases

```python
from radis.io.hitemp import fetch_hitemp
fetch_hitemp("OH", engine="pytables")       #  (if testing : make sure there is no vaex-formatted hdf5 already)
import radis; radis.AUTO_UPDATE_DATABASE = True
fetch_hitemp("OH", engine="vaex")   #  auto-converts to vaex-compatible HDF5

```

- [x] re-add comparison with already registered database (if exists)



Will be merged / rebased on #281 . 

Fixes #239 #295 


### TODO : 

- [ ] add multiple isotope support when querying vaex files
- [ ] add an auto-update test 
- [ ] write the parsed HITEMP files in Vaex-compatible / H5py format (see e3a1d734) 
- [x] Feature idea : for CO2 / H2O (multiple files) : only parse & download the files needed. Use the metadata keys wavenum_min; wavenum_max (?) in the database to see if new files should be downloaded ?  Could be extended for CDSD-4000 later.

Identified problems : 
- process sometimes stuck in a Spyder Ipython console (while converting to VAEX). File is converted properly, though, but only appears after the console is closed.   Similar to #338 and https://github.com/spyder-ide/spyder/issues/16183
- Some HDF5 files not closed properly if a crash happens : https://github.com/pandas-dev/pandas/issues/8692




### Performance improvement : 

Old (full loading and partial loading of HITEMP CO)

```python

In [2]: %timeit fetch_hitemp('CO', engine='pytables')
4.84 s ± 77 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [3]: %timeit fetch_hitemp('CO', engine='pytables', load_wavenum_min=2000, load_wavenum_max=3000)
814 ms ± 28 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

```

New:

```python
In [2]: %timeit fetch_hitemp('CO', engine='vaex')
1.42 s ± 35.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [3]: %timeit fetch_hitemp('CO', engine='vaex', load_wavenum_min=2000, load_wavenum_max=3000)
519 ms ± 27.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

```

So 1.5 to 3.5x faster. 

refactor : MolParams
### 🎯 Todo
<!--- A clear and concise description of what you want to happen. -->

MolParams (molecular parameters) is currently reinitialized at several places through a calculation. A same object should be used. 
Will also make it possible to use #301 

### 💡 Possible implementations
<!--- Not obligatory, but suggest an idea for implementing addition or change -->

- Store as an attribute of SpectrumFactory ? 



---
<sub>*Tip : you can vote for [New Features](https://github.com/radis/radis/issues/79) on https://feathub.com/radis/radis*</sub>

ExoMol Support
### 🎯 Describe what you'd like
<!--- A clear and concise description of what you want to happen. -->

ExoMol support. Asked for in #34 

### Implementation (see https://github.com/radis/radis/pull/320)

```python
from radis import SpectrumFactory

sf = SpectrumFactory(
        4310,
        4320,
        molecule="H2O",
        isotope="1",
        verbose=3,
    )
sf.fetch_databank("exomol")
s = sf.eq_spectrum(500)
s.plot()
```

![image](https://user-images.githubusercontent.com/16088743/126080498-5f437d32-c274-47bf-a17a-995ad7f3d4a7.png)

### 💡 Possible implementations
<!--- Not obligatory, but suggest an idea for implementing addition or change -->

Borrowing and adapting some codes from the [ExoJax](https://github.com/HajimeKawahara/exojax) project from @HajimeKawahara (which you should also have a look at !) Hopefully we can end up maintaining a unified code base for local database formats (currently : ExoJax uses Feather format; which is extremelly efficient but cannot be masked/filtered as easily as HDF5)





---
<sub>*Tip : you can vote for [New Features](https://github.com/radis/radis/issues/79) on https://feathub.com/radis/radis*</sub>



Discussion / hints:

Fix patch:
diff --git a/docs/features/features.rst b/docs/features/features.rst
index 34082a06d..249722f94 100644
--- a/docs/features/features.rst
+++ b/docs/features/features.rst
@@ -13,8 +13,8 @@ RADIS is both an infrared :ref:`line-by-line code <label_line_by_line>`
 and a :ref:`post-processing library <label_spectrum>`.
 It includes:
 
-- Absorption and emission spectra of all [HITRAN-2016]_ species under equilibrium calculations (:py:data:`~radis.db.MOLECULES_LIST_EQUILIBRIUM`)
-- Absorption and emission spectra of CO2 and CO for nonequilibrium calculations (see :py:data:`~radis.db.MOLECULES_LIST_NONEQUILIBRIUM` )
+- Absorption and emission spectra of all [HITRAN-2016]_ and [ExoMol-2020]_ species under equilibrium calculations (:py:data:`~radis.db.MOLECULES_LIST_EQUILIBRIUM`)
+- Absorption and emission spectra of CO2 and CO for non-LTE calculations (see :py:data:`~radis.db.MOLECULES_LIST_NONEQUILIBRIUM` )
 - Different Line Databases: support of [HITRAN-2016]_, [HITEMP-2010]_ and [CDSD-4000]_ line databases (see :py:data:`~radis.lbl.loader.KNOWN_DBFORMAT`)
 - Calculation of :ref:`Rovibrational Energies of molecules <label_examples_rovibrational_energies>`.
 - Calculation of equilibrium and nonequilibrium :ref:`Partition Functions <label_examples_partition_functions>`.
@@ -92,6 +92,7 @@ List of supported line databases formats: :py:data:`~radis.lbl.loader.KNOWN_DBFO
 - [HITRAN-2016]_
 - [HITEMP-2010]_
 - [CDSD-4000]_
+- [ExoMol-2020]_
 
 For download and configuration of line databases, see the :ref:`Line Databases section <label_line_databases>`
 
diff --git a/docs/references/_bibliography.rst b/docs/references/_bibliography.rst
index 677cb96fd..c070e6183 100644
--- a/docs/references/_bibliography.rst
+++ b/docs/references/_bibliography.rst
@@ -52,4 +52,8 @@ The latest HITEMP database version is automatically downloaded if using ``databa
                spectroscopic databank, Journal of Quantitative Spectroscopy and Radiative Transfer 112 (9) (2011)
                1403–1410, ISSN 00224073, `doi:10.1016/j.jqsrt.2011.03.005 <https://www.sciencedirect.com/science/article/pii/S0022407311001154>`__
 
+.. [ExoMol-2020] Tennyson et al., The 2020 release of the ExoMol database: Molecular line lists for
+                exoplanet and other hot atmospheres, Journal of Quantitative Spectroscopy and Radiative Transfer 255,
+                (2020), 107228,  `doi:https://doi.org/10.1016/j.jqsrt.2020.107228 <https://www.sciencedirect.com/science/article/abs/pii/S002240732030491X>`__
+
 For download and configuration of line databases, see the :ref:`Line Databases section <label_line_databases>`
diff --git a/environment.yml b/environment.yml
index 0dadd0791..89c38f278 100644
--- a/environment.yml
+++ b/environment.yml
@@ -25,6 +25,9 @@ dependencies:
 - psutil # for getting user RAM
 - pip
 - pip:
+  - beautifulsoup4 # parse ExoMol website
+  - lxml           # parser used for ExoMol website
+  - pyarrow       # for the feather format (temporarily needed for ExoMol)
   - publib>=0.3.2  # Plotting styles for Matplotlib
   - hitran-api     # HAPI, used to access TIPS partition functions
   - json-tricks>=3.15.0  # to deal with non jsonable formats
diff --git a/examples/plot_exomol_spectrum.py b/examples/plot_exomol_spectrum.py
new file mode 100644
index 000000000..801accedb
--- /dev/null
+++ b/examples/plot_exomol_spectrum.py
@@ -0,0 +1,25 @@
+# -*- coding: utf-8 -*-
+"""
+================================
+Calculate a spectrum from ExoMol
+================================
+
+Auto-download and compute a SiO spectrum from the ExoMol database ([ExoMol-2020]_)
+
+"""
+
+from radis import calc_spectrum
+
+s = calc_spectrum(
+    1080,
+    1320,  # cm-1
+    molecule="SiO",
+    isotope="1",
+    pressure=1.01325,  # bar
+    Tgas=1000,  # K
+    mole_fraction=0.1,
+    path_length=1,  # cm
+    databank="exomol",  # uses the recommended database. Use ('exomol', "EBJT") for a specific database ("EBJT")
+)
+s.apply_slit(1, "cm-1")  # simulate an experimental slit
+s.plot("radiance")
diff --git a/examples/plot_hitemp_spectrum.py b/examples/plot_hitemp_spectrum.py
new file mode 100644
index 000000000..46ac2db92
--- /dev/null
+++ b/examples/plot_hitemp_spectrum.py
@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+"""
+================================
+Calculate a spectrum from HITEMP
+================================
+
+Auto-download and compute a CO spectrum from HITEMP ([HITEMP-2010]_).
+
+Compare with HITRAN-2016.
+
+"""
+
+from radis import calc_spectrum
+
+s = calc_spectrum(
+    2135,
+    2170,  # cm-1
+    molecule="CO",
+    isotope="1",
+    pressure=3,  # bar
+    Tgas=2000,  # K
+    mole_fraction=0.1,
+    path_length=1,  # cm
+    databank="hitemp",  # latest version is fetched
+    name="HITEMP-2020",
+)
+
+from radis.test.utils import setup_test_line_databases
+
+setup_test_line_databases()
+
+s2 = calc_spectrum(
+    2135,
+    2170,  # cm-1
+    molecule="CO",
+    isotope="1",
+    pressure=3,  # bar
+    Tgas=2000,  # K
+    mole_fraction=0.1,
+    path_length=1,  # cm
+    databank="HITRAN-CO-TEST",  # we could also have fetched the latest HITRAN with "databank='hitran'"
+    name="HITRAN-2016",
+)
+
+# Compare:
+from radis import plot_diff
+
+plot_diff(s, s2)
diff --git a/radis/__init__.py b/radis/__init__.py
index e53ffed24..b426f823e 100644
--- a/radis/__init__.py
+++ b/radis/__init__.py
@@ -186,41 +186,3 @@ def get_version(verbose=False, add_git_number=True):
 __all__.extend(phys.__all__)
 __all__.extend(spectrum.__all__)
 __all__.extend(tools.__all__)
-
-# from .db import (  # database of molecules
-#     MOLECULES_LIST_EQUILIBRIUM,
-#     MOLECULES_LIST_NONEQUILIBRIUM,
-#     Molecules,
-#     getMolecule,
-# )
-# from .io import cdsd2df, hit2df  # input / output
-# from .lbl import LevelsList, SpectrumFactory, calc_spectrum  # line-by-line module
-
-# # from .levels import *  # rovibrational energies and partition functions
-# from .los import MergeSlabs, SerialSlabs  # line-of-sight module
-# from .phys import planck, planck_wn, sPlanck  # conversion functions, blackbody objects
-# from .spectrum import Spectrum  # Spectrum object
-# from .spectrum import (
-#     Radiance,
-#     Radiance_noslit,
-#     Transmittance,
-#     Transmittance_noslit,
-#     calculated_spectrum,
-#     experimental_spectrum,
-#     get_baseline,
-#     get_diff,
-#     plot_diff,
-#     sub_baseline,
-#     transmittance_spectrum,
-# )
-
-# # from .test import *  # test
-# from .tools import (  # get_FWHW, get_effective_FWHW,        # doesnt work for some reason? retry. TODO; slit, database, line survey, etc.
-#     SpecDatabase,
-#     crop_slit,
-#     get_eq_mole_fraction,
-#     load_spec,
-#     plot_slit,
-#     plot_spec,
-#     recenter_slit,
-# )
diff --git a/radis/db/__init__.py b/radis/db/__init__.py
index 17630eba9..2a9542144 100644
--- a/radis/db/__init__.py
+++ b/radis/db/__init__.py
@@ -12,9 +12,10 @@ def _get_supported_molecules_equilibrium():
     """Molecules supported in RADIS equilibrium calculations Basically, all
     [HITRAN-2016]_ species."""
     # Add all HITRAN species
-    from .classes import HITRAN_MOLECULES
+    # Add ExoMol species
+    from .classes import EXOMOL_MOLECULES, HITRAN_MOLECULES
 
-    return HITRAN_MOLECULES
+    return list(set(HITRAN_MOLECULES).union(set(EXOMOL_MOLECULES)))
 
 
 def _get_supported_molecules_nonequilibrium():
@@ -32,9 +33,10 @@ def _get_supported_molecules_nonequilibrium():
     __supported_molecules_equilibrium__
 ) = _get_supported_molecules_equilibrium()
 """ list: molecules that can be calculated in RADIS at equilibrium.
-All [HITRAN-2016]_ species are available.
+All [HITRAN-2016]_ and [ExoMol-2020]_ species are available.
+
 Absorption coefficient calculated with RADIS at 300 K, 1 atm are shown for all
-these molecules in the :ref:`Examples page <label_examples_hitran_spectra>` (HITRAN spectra).
+[HITRAN-2016]_ molecules in the :ref:`HITRAN spectra page <label_examples_hitran_spectra>` .
 
 - 1 	``'H2O'`` : 	Water 	(`spectrum <https://raw.githubusercontent.com/radis/radis-examples/master/hitran_spectra/out/0%20-%20H2O%20infrared%20spectrum.png>`__)
 - 2 	``'CO2'`` : 	Carbon Dioxide    (`spectrum <https://raw.githubusercontent.com/radis/radis-examples/master/hitran_spectra/out/1%20-%20CO2%20infrared%20spectrum.png>`__)
diff --git a/radis/db/classes.py b/radis/db/classes.py
index aa584311a..60aef7a42 100644
--- a/radis/db/classes.py
+++ b/radis/db/classes.py
@@ -212,25 +212,21 @@
 
 def get_molecule_identifier(molecule_name):
     r"""
-    For a given input molecular formula, return the corresponding HITRAN molecule
-    identifier number [1]_.
-
+    For a given input molecular formula, return the corresponding
+    :py:data:`~radis.db.classes.HITRAN_MOLECULES` identifier number [1]_.
 
     Parameters
     ----------
     molecular_formula : str
         The string describing the molecule.
 
-
     Returns
     -------
     M: int
         The HITRAN molecular identified number.
 
-
     References
     ----------
-
     .. [1] `HITRAN 1996, Rothman et al., 1998 <https://www.sciencedirect.com/science/article/pii/S0022407398000788>`__
 
     Function is from from https://github.com/nzhagen/hitran/blob/master/hitran.py
@@ -252,20 +248,17 @@ def get_molecule_identifier(molecule_name):
 
 def get_molecule(molecule_id):
     r"""
-    For a given input molecular identifier, return the corresponding HITRAN
-    molecule name [1]_.
-
+    For a given input molecular identifier, return the corresponding
+    :py:data:`~radis.db.classes.HITRAN_MOLECULES` name [1]_.
 
     Parameters
     ----------
-
     molecular_id: str
         Hitran identifier of the molecule.
 
 
     References
     ----------
-
     .. [1] `HITRAN 1996, Rothman et al., 1998 <https://www.sciencedirect.com/science/article/pii/S0022407398000788>`__
 
     """
@@ -281,6 +274,19 @@ def get_molecule(molecule_id):
         )
 
 
+# %% ExoMol molecules
+
+
+KNOWN_EXOMOL_ISOTOPES_NAMES = {
+    ("FeH", 1): "56Fe-1H",
+    ("SiO", 1): "28Si-16O",  # Placeholder until all molecules parsed from website TODO
+}
+"""All :py:data:`~radis.db.classes.HITRAN_MOLECULES` are also converted to their ExoMol full-name format
+in :py:func:`~radis.io.exomol.get_exomol_full_isotope_name`
+"""
+EXOMOL_MOLECULES = list(set([M for M, iso in KNOWN_EXOMOL_ISOTOPES_NAMES.keys()]))
+
+
 # %% Molecule class
 
 
@@ -289,19 +295,16 @@ class Molecule(object):
 
     Parameters
     ----------
-
     name: str, int
         molecule name, or HITRAN identifier
 
     Other Parameters
     ----------------
-
     verbose: boolean
         more blabla
 
     See Also
     --------
-
     :py:class:`~radis.db.classes.Isotope`, :py:class:`~radis.db.classes.ElectronicState`,
     :py:func:`~radis.db.molecules.getMolecule`
 
diff --git a/radis/db/molecules.py b/radis/db/molecules.py
index 4a516337d..f930ab468 100644
--- a/radis/db/molecules.py
+++ b/radis/db/molecules.py
@@ -146,7 +146,9 @@
 """
 
 
-def getMolecule(molecule, isotope=None, electronic_state=None, verbose=True):
+def getMolecule(
+    molecule, isotope=None, electronic_state=None, verbose=True
+) -> ElectronicState:
     """Get an :py:class:`~radis.db.classes.ElectronicState` object in the
     RADIS :py:data:`~radis.db.molecules.Molecules` list, which use the defaults
     :ref:`spectroscopic constants <label_db_spectroscopic_constants>`.
@@ -166,8 +168,7 @@ def getMolecule(molecule, isotope=None, electronic_state=None, verbose=True):
 
     Returns
     -------
-    ElectronicState
-        an :py:class:`~radis.db.classes.ElectronicState` object.
+    ElectronicState: an :py:class:`~radis.db.classes.ElectronicState` object.
 
     .. minigallery:: radis.db.molecules.getMolecule
         :add-heading:
diff --git a/radis/db/molparam.py b/radis/db/molparam.py
index 9ae3bc94e..bf94e34fe 100644
--- a/radis/db/molparam.py
+++ b/radis/db/molparam.py
@@ -7,11 +7,143 @@
 """
 
 
+import re
+
 import pandas as pd
 
-from radis.db.classes import get_molecule_identifier
 from radis.db.utils import getFile
 
+# Isotope definitions in HITRAN
+isotope_name_dict = {
+    (1, 1): "H2(16O)",
+    (1, 2): "H2(18O)",
+    (1, 3): "H2(17O)",
+    (1, 4): "HD(16O)",
+    (1, 5): "HD(18O)",
+    (1, 6): "HD(17O)",
+    (1, 7): "D2(16O)",
+    (2, 1): "(12C)(16O)2",
+    (2, 2): "(13C)(16O)2",
+    (2, 3): "(16O)(12C)(18O)",
+    (2, 4): "(16O)(12C)(17O)",
+    (2, 5): "(16O)(13C)(18O)",
+    (2, 6): "(16O)(13C)(17O)",
+    (2, 7): "(12C)(18O)2",
+    (2, 8): "(17O)(12C)(18O)",
+    (2, 9): "(12C)(17O)2",
+    (2, 10): "(13C)(18O)2",
+    (2, 11): "(18O)(13C)(17O)",
+    (2, 12): "(13C)(17O)2",
+    (3, 1): "(16O)3",
+    (3, 2): "(16O)(16O)(18O)",
+    (3, 3): "(16O)(18O)(16O)",
+    (3, 4): "(16O)(16O)(17O)",
+    (3, 5): "(16O)(17O)(16O)",
+    (4, 1): "(14N)2(16O)",
+    (4, 2): "(14N)(15N)(16O)",
+    (4, 3): "(15N)(14N)(16O)",
+    (4, 4): "(14N)2(18O)",
+    (4, 5): "(14N)2(17O)",
+    (5, 1): "(12C)(16O)",
+    (5, 2): "(13C)(16O)",
+    (5, 3): "(12C)(18O)",
+    (5, 4): "(12C)(17O)",
+    (5, 5): "(13C)(18O)",
+    (5, 6): "(13C)(17O)",
+    (6, 1): "(12C)H4",
+    (6, 2): "(13C)H4",
+    (6, 3): "(12C)H3D",
+    (6, 4): "(13C)H3D",
+    (7, 1): "(16O)2",
+    (7, 2): "(16O)(18O)",
+    (7, 3): "(16O)(17O)",
+    (8, 1): "(14N)(16O)",
+    (8, 2): "(15N)(16O)",
+    (8, 3): "(14N)(18O)",
+    (9, 1): "(32S)(16O)2",
+    (9, 2): "(34S)(16O)2",
+    (10, 1): "(14N)(16O)2",
+    (11, 1): "(14N)H3",
+    (11, 2): "(15N)H3",
+    (12, 1): "H(14N)(16O)3",
+    (12, 2): "H(15N)(16O)3",
+    (13, 1): "(16O)H",
+    (13, 2): "(18O)H",
+    (13, 3): "(16O)D",
+    (14, 1): "H(19F)",
+    (14, 2): "D(19F)",
+    (15, 1): "H(35Cl)",
+    (15, 2): "H(37Cl)",
+    (15, 3): "D(35Cl)",
+    (15, 4): "D(37Cl)",
+    (16, 1): "H(79Br)",
+    (16, 2): "H(81Br)",
+    (16, 3): "D(79Br)",
+    (16, 4): "D(81Br)",
+    (17, 1): "H(127I)",
+    (17, 2): "D(127I)",
+    (18, 1): "(35Cl)(16O)",
+    (18, 2): "(37Cl)(16O)",
+    (19, 1): "(16O)(12C)(32S)",
+    (19, 2): "(16O)(12C)(34S)",
+    (19, 3): "(16O)(13C)(32S)",
+    (19, 4): "(16O)(12C)(33S)",
+    (19, 5): "(18O)(12C)(32S)",
+    (20, 1): "H2(12C)(16O)",
+    (20, 2): "H2(13C)(16O)",
+    (20, 3): "H2(12C)(18O)",
+    (21, 1): "H(16O)(35Cl)",
+    (21, 2): "H(16O)(37Cl)",
+    (22, 1): "(14N)2",
+    (22, 2): "(14N)(15N)",
+    (23, 1): "H(12C)(14N)",
+    (23, 2): "H(13C)(14N)",
+    (23, 3): "H(12C)(15N)",
+    (24, 1): "(12C)H3(35Cl)",
+    (24, 2): "(12C)H3(37Cl)",
+    (25, 1): "H2(16O)2",
+    (26, 1): "(12C)2H2",
+    (26, 2): "(12C)(13C)H2",
+    (26, 3): "(12C)2HD",
+    (27, 1): "(12C)2H6",
+    (27, 2): "(12C)H3(13C)H3",
+    (28, 1): "(31P)H3",
+    (29, 1): "(12C)(16O)(19F)2",
+    (29, 2): "(13C)(16O)(19F)2",
+    (30, 1): "(32S)(19F)6",
+    (31, 1): "H2(32S)",
+    (31, 2): "H2(34S)",
+    (31, 3): "H2(33S)",
+    (32, 1): "H(12C)(16O)(16O)H",
+    (33, 1): "H(16O)2",
+    (34, 1): "(16O)",
+    (36, 1): "(14N)(16O)+",
+    (37, 1): "H(16O)(79Br)",
+    (37, 2): "H(16O)(81Br)",
+    (38, 1): "(12C)2H4",
+    (38, 2): "(12C)H2(13C)H2",
+    (39, 1): "(12C)H3(16O)H",
+    (40, 1): "(12C)H3(79Br)",
+    (40, 2): "(12C)H3(81Br)",
+    (41, 1): "(12C)H3(12C)(14N)",
+    (42, 1): "(12C)(19F)4",
+    (43, 1): "(12C)4H2",
+    (44, 1): "H(12C)3(14N)",
+    (45, 1): "H2",
+    (45, 2): "HD",
+    (46, 1): "(12C)(32S)",
+    (46, 2): "(12C)(34S)",
+    (46, 3): "(13C)(32S)",
+    (46, 4): "(12C)(33S)",
+    (47, 1): "(32S)(16O)3",
+    (48, 1): "(12C)2(14N)2",
+    (49, 1): "(12C)(16O)(35Cl)2",
+    (49, 2): "(12C)(16O)(35Cl)(37Cl)",
+}
+# TODO : add test in test_molparams comparing dictionary above to
+# latest values from HAPI; to make sure we never have any problem
+# even if HITRAN eventually changes the conventions and labels.
+
 
 class MolParams:
     def __init__(self, file=None):
@@ -24,7 +156,7 @@ def __init__(self, file=None):
 
         Examples
         --------
-        Get abundance of CO2, isotope 1::
+        Get earth abundance of CO2, isotope 1::
 
             molpar = Molparams()
             molpar.get(2, 1, 'abundance')         # 2 for CO2, 1 for isotope 1
@@ -36,6 +168,7 @@ def __init__(self, file=None):
 
         References
         ----------
+
         http://hitran.org/media/molparam.txt
         """
 
@@ -45,6 +178,8 @@ def __init__(self, file=None):
         df = pd.read_csv(file, comment="#", delim_whitespace=True)
         df = df.set_index(["id", "iso"])
 
+        df["isotope_name_exomol"] = _add_exomol_name(df)
+
         self.df = df
 
         # ------
@@ -54,36 +189,58 @@ def __init__(self, file=None):
             pass  # old Python version
 
     def get(self, M, I, key):
-        """
+        """Get attribute of molecule, isotope.
+
         Parameters
         ----------
-        M: int or str
-            molecule id or molecule name
+        M: int
+            molecule id
         I: int
             molecule isotope #
-        key: ``'abundance'``, ``'mol_mass'``
+        key: ``'abundance'``, ``'mol_mass'``, ``'isotope_name_exomol'``,
             parameter
+
         """
-        if isinstance(M, str):
+        try:
+            float(M)
+        except ValueError:
+            # not a number
+            from radis.db.classes import get_molecule_identifier
+
             M = get_molecule_identifier(M)
         return self.df.loc[(M, I), key]
 
 
-def _test(verbose=True, *args, **kwargs):
-    molpar = MolParams()
+def _add_exomol_name(df):
+    """Convert HITRAN isotope name from :py:data:`~radis.db.molparam.isotope_name_dict`
+    in ExoMol format
+    """
+
+    df["isotope_name"] = df.index.map(isotope_name_dict)
 
-    b = True
+    # add exomol full name:
+    def replace_parenthesis(s):
+        s = str(s)
+        # add parenthesis for hydrogen : Eg. H2 -> (1H)2,   H->(1H)
+        s = re.sub("H([1-9])*", r"(1H)\1", s)
+        # add parenthesis for deuterium : Eg. D2 -> (2H)2,   D->(2H)
+        s = re.sub("D([1-9])*", r"(2H)\1", s)
 
-    b *= molpar.get(2, 2, "abundance") == 0.0110574
-    if verbose:
-        print("CO2-636 abundance:", molpar.get(2, 2, "abundance"))
+        # Cut elements and replace parenthesis with - :
+        elements = re.findall("\(.*?\)[1-9]*", s)
+        elements = [e.replace("(", "").replace(")", "") for e in elements]
+        s = "-".join(elements)
 
-    if verbose:
-        print("Testing molparams.py:", b)
+        return s
 
-    return True
+    return df.isotope_name.apply(replace_parenthesis)
 
 
 if __name__ == "__main__":
+    from radis.test.io.test_exomol import test_exomol_parsing_functions
+
+    test_exomol_parsing_functions()
+
+    from radis.test.db.test_molparams import test_molparams
 
-    _test()
+    test_molparams()
diff --git a/radis/db/molparam.txt b/radis/db/molparam.txt
index d35069389..ab0d74200 100644
--- a/radis/db/molparam.txt
+++ b/radis/db/molparam.txt
@@ -10,7 +10,7 @@
 # Note
 # ---------
 # Isotope number was derived manually assuming the isonames were ordered in the database
-# The isotope name (ex: CO2 626) is kept for comparison if ever needed
+# The isotope AFGL code (ex: CO2 626) is kept for comparison if ever needed
 #
 # -----------------------------------------------------
 id iso isoname abundance     Q_296K      gj    mol_mass
@@ -183,4 +183,4 @@ id iso isoname abundance     Q_296K      gj    mol_mass
 46 3  32  .105565E-01    5.3750E+02    2     44.974368
 46 4  23  .741668E-02    1.0230E+03    4     44.970399
 #   SO3 (47)
-47 5  26  .943400E+00    7.9638E+03    1     79.956820
+47 1  26  .943400E+00    7.9638E+03    1     79.956820
diff --git a/radis/io/exomol.py b/radis/io/exomol.py
new file mode 100644
index 000000000..20b52826e
--- /dev/null
+++ b/radis/io/exomol.py
@@ -0,0 +1,709 @@
+"""Molecular database (MDB) class
+
+   * MdbExomol is the MDB for ExoMol
+
+Initial code borrowed from the `Exojax <https://github.com/HajimeKawahara/exojax>`__
+code (which you should also have a look at !), by @HajimeKawahara, under MIT License.
+
+"""
+import pathlib
+import warnings
+
+import numpy as np
+
+# from exojax.spec import hapi, exomolapi, exomol
+# from exojax.spec.hitran import gamma_natural as gn
+import pandas as pd
+
+try:
+    from . import exomolapi
+    from .exomol_utils import e2s
+except ImportError:  # if local import
+    from radis.io import exomolapi
+    from radis.io.exomol_utils import e2s
+
+from radis.db.classes import get_molecule_identifier
+
+EXOMOL_URL = "http://www.exomol.com/db/"
+
+from urllib.request import HTTPError, urlopen
+
+from bs4 import BeautifulSoup
+
+
+def get_exomol_full_isotope_name(molecule, isotope):
+    """Get full isotope name for ``molecule`` and ``isotope`` number
+
+    Parameters
+    ----------
+    molecule: str
+    isotope: int
+        terrestrial abundance
+
+    Examples
+    --------
+    ::
+        get_exomol_full_isotope_name("CH4", 1)
+        >>> '12C-1H4'
+
+    See Also
+    --------
+    :py:func:`~radis.io.exomol.get_exomol_database_list`"""
+
+    from radis.db.classes import KNOWN_EXOMOL_ISOTOPES_NAMES
+
+    if (molecule, isotope) in KNOWN_EXOMOL_ISOTOPES_NAMES:
+        return KNOWN_EXOMOL_ISOTOPES_NAMES[(molecule, isotope)]
+    else:
+        # Read and convert from HITRAN molecules
+        from radis.db.molparam import MolParams
+
+        mp = MolParams()
+        return mp.get(molecule, isotope, "isotope_name_exomol")
+
+
+def get_exomol_database_list(molecule, isotope_full_name):
+    """Parse ExoMol website and return list of available databases, and recommended database
+
+    Parameters
+    ----------
+    molecule: str
+    isotope_full_name: str
+        isotope full name (ex. ``12C-1H4`` for CH4,1). Get it from
+        :py:func:`radis.io.exomol.get_exomol_full_isotope_name`
+
+    Returns
+    -------
+
+    Examples
+    --------
+    ::
+        databases, recommended = get_exomol_database_list("CH4", "12C-1H4")
+        >>> ['xsec-YT10to10', 'YT10to10', 'YT34to10'], 'YT34to10'
+
+    See Also
+    --------
+    :py:func:`~radis.io.exomol.get_exomol_full_isotope_name`
+    """
+
+    url = f"https://exomol.com/data/molecules/{molecule}/{isotope_full_name}"
+    try:
+        response = urlopen(url).read()
+    except HTTPError as err:
+        raise ValueError(f"HTTPError opening url={url}") from err
+
+    soup = BeautifulSoup(
+        response, features="lxml"
+    )  # make soup that is parse-able by bs
+
+    # Recommended database
+    rows = soup.find_all(
+        "a", {"class": "list-group-item link-list-group-item recommended"}
+    )
+    databases_recommended = [r.get_attribute_list("title")[0] for r in rows]
+
+    # All others
+    rows = soup.find_all("a", {"class": "list-group-item link-list-group-item"})
+    databases = [r.get_attribute_list("title")[0] for r in rows]
+
+    assert len(databases_recommended) <= 1
+
+    databases = databases + databases_recommended
+
+    return databases, databases_recommended[0]
+
+
+def fetch_exomol(
+    molecule,
+    database=None,
+    local_databases=r"~/.radisdb/exomol/",
+    databank_name="EXOMOL-{molecule}",
+    isotope=None,
+    load_wavenum_min=None,
+    load_wavenum_max=None,
+    cache=True,
+    verbose=True,
+    clean_cache_files=True,
+    return_local_path=False,
+    return_partition_function=False,
+):
+    """Stream ExoMol file from EXOMOL website. Unzip and build a HDF5 file directly.
+
+    Returns a Pandas DataFrame containing all lines.
+
+    Parameters
+    ----------
+    molecule: str
+        ExoMol molecule
+    database: str
+        database name. Ex:: ``POKAZATEL`` or ``BT2`` for ``H2O``. See
+        :py:data:`~radis.io.exomol.KNOWN_EXOMOL_DATABASE_NAMES`. If ``None`` and
+        there is only one database available, use it.
+    local_databases: str
+        where to create the RADIS HDF5 files. Default ``"~/.radisdb/exomol"``
+    databank_name: str
+        name of the databank in RADIS :ref:`Configuration file <label_lbl_config_file>`
+        Default ``"EXOMOL-{molecule}"``
+    isotope: str
+        load only certain isotopes : ``'2'``, ``'1,2'``, etc. Default ``1``.
+    load_wavenum_min, load_wavenum_max: float (cm-1)
+        load only specific wavenumbers.
+
+    Other Parameters
+    ----------------
+    cache: bool, or ``'regen'``
+        if ``True``, use existing HDF5 file. If ``False`` or ``'regen'``, rebuild it.
+        Default ``True``.
+    verbose: bool
+    clean_cache_files: bool
+        if ``True`` clean downloaded cache files after HDF5 are created.
+    return_local_path: bool
+        if ``True``, also returns the path of the local database file.
+    return_partition_function: bool
+        if ``True``, also returns a :py:class:`~radis.levels.partfunc.PartFuncExoMol` object.
+
+    Returns
+    -------
+    df: pd.DataFrame
+        Line list
+        A HDF5 file is also created in ``local_databases`` and referenced
+        in the :ref:`RADIS config file <label_lbl_config_file>` with name
+        ``databank_name``
+    local_path: str
+        path of local database file if ``return_local_path``
+
+    Notes
+    -----
+    if using ``load_only_wavenum_above/below`` or ``isotope``, the whole
+    database is anyway downloaded and uncompressed to ``local_databases``
+    fast access .HDF5 files (which will take a long time on first call). Only
+    the expected wavenumber range & isotopes are returned. The .HFD5 parsing uses
+    :py:func:`~radis.io.hdf5.hdf2df`
+
+    See Also
+    --------
+    :py:func:`~radis.io.hdf5.hdf2df`
+
+    """
+    full_molecule_name = get_exomol_full_isotope_name(molecule, isotope)
+    known_exomol_databases, recommended_database = get_exomol_database_list(
+        molecule, full_molecule_name
+    )
+
+    _exomol_use_hint = f"Select one of them with `radis.fetch_exomol(DATABASE_NAME)`, `SpectrumFactory.fetch_databank('exomol', exomol_database=DATABASE_NAME')`, or `calc_spectrum(..., databank=('exomol', DATABASE_NAME))` \n"
+    if database is None:
+        if len(known_exomol_databases) == 1:
+            database = known_exomol_databases[0]
+            if verbose:
+                print(f"Using ExoMol database {database} for {full_molecule_name}")
+        elif recommended_database:
+            database = recommended_database
+            if verbose:
+                print(
+                    f"Using ExoMol database {database} for {full_molecule_name} (recommended by the ExoMol team). All available databases are {known_exomol_databases}. {_exomol_use_hint}"
+                )
+        else:
+            raise KeyError(
+                f"Choose one of the several databases available for {full_molecule_name} in ExoMol: {known_exomol_databases}. ({recommended_database} is recommended by the ExoMol team). {_exomol_use_hint}"
+            )
+    else:
+        if database not in known_exomol_databases:
+            raise KeyError(
+                f"{database} is not of the known available ExoMol databases for {full_molecule_name}. Choose one of : {known_exomol_databases}. ({recommended_database} is recommended by the ExoMol team). {_exomol_use_hint}"
+            )
+
+    local_path = (
+        pathlib.Path(local_databases).expanduser()
+        / molecule
+        / full_molecule_name
+        / database
+    )
+
+    mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max])
+
+    df = mdb.to_df()
+
+    out = [df]
+
+    if return_local_path:
+        out.append(str(mdb.path))
+    if return_partition_function:
+        assert return_local_path
+        out.append(mdb.to_partition_function_tabulator())
+
+    return out
+
+
+class MdbExomol(object):
+    """molecular database of ExoMol
+
+    MdbExomol is a class for ExoMol.
+
+    Attributes:
+        nu_lines (nd array): line center (cm-1)
+        Sij0 (nd array): line strength at T=Tref (cm)
+        dev_nu_lines (np array): line center in device (cm-1)
+        logsij0 (np array): log line strength at T=Tref
+        A (np array): Einstein A coefficient
+        elower (np array): the lower state energy (cm-1)
+        gpp (np array): statistical weight
+        jlower (np array): J_lower
+        jupper (np array): J_upper
+        n_Tref (np array): temperature exponent
+        alpha_ref (np array): alpha_ref (gamma0)
+        n_Tref_def: default temperature exponent in .def file, used for jlower not given in .broad
+        alpha_ref_def: default alpha_ref (gamma0) in .def file, used for jlower not given in .broad
+
+    """
+
+    def __init__(
+        self,
+        path,
+        nurange=[-np.inf, np.inf],
+        margin=1.0,
+        crit=-np.inf,
+        bkgdatm="H2",
+        broadf=True,
+    ):
+        """Molecular database for Exomol form
+
+        Args:
+           path: path for Exomol data directory/tag. For instance, "/home/CO/12C-16O/Li2015"
+           nurange: wavenumber range list (cm-1) or wavenumber array
+           margin: margin for nurange (cm-1)
+           crit: line strength lower limit for extraction
+           bkgdatm: background atmosphere for broadening. e.g. H2, He,
+           broadf: if False, the default broadening parameters in .def file is used
+
+        Note:
+           The trans/states files can be very large. For the first time to read it, we convert it to the feather-format. After the second-time, we use the feather format instead.
+
+        """
+        explanation = "Note: Couldn't find the feather format. We convert data to the feather format. After the second time, it will become much faster."
+
+        self.path = pathlib.Path(path)
+        t0 = self.path.parents[0].stem
+        molec = t0 + "__" + str(self.path.stem)
+        self.bkgdatm = bkgdatm
+        print("Background atmosphere: ", self.bkgdatm)
+        molecbroad = t0 + "__" + self.bkgdatm
+
+        self.crit = crit
+        self.margin = margin
+        self.nurange = [np.min(nurange), np.max(nurange)]
+        self.broadf = broadf
+        # Where exomol files are
+        self.states_file = self.path / pathlib.Path(molec + ".states.bz2")
+        self.pf_file = self.path / pathlib.Path(molec + ".pf")
+        self.def_file = self.path / pathlib.Path(molec + ".def")
+        self.broad_file = self.path / pathlib.Path(molecbroad + ".broad")
+
+        if not self.def_file.exists():
+            self.download(molec, extension=[".def"])
+        if not self.pf_file.exists():
+            self.download(molec, extension=[".pf"])
+        if not self.states_file.exists():
+            self.download(molec, extension=[".states.bz2"])
+        if not self.broad_file.exists():
+            self.download(molec, extension=[".broad"])
+
+        # Add molecule name
+        tag = molec.split("__")
+        self.isotope_fullname = tag[0]
+        self.molecule = e2s(tag[0])
+        self.isotope = 1  # Placeholder. TODO : impement parsing of other isotopes.
+
+        # load def
+        # @minou: I feel this next line should just update the self.stuff
+        dic_def = exomolapi.read_def(self.def_file)  # approx. 3 ms
+        self.n_Texp_def = dic_def["n_Texp"]
+        self.alpha_ref_def = dic_def["alpha_ref"]
+        self.molmass = dic_def["molmass"]
+
+        #  default n_Texp value if not given
+        if self.n_Texp_def is None:
+            warnings.warn(
+                Warning(
+                    """
+                    No default broadening exponent in def file. Assigned n = 0.5
+                    """
+                )
+            )
+            self.n_Texp_def = 0.5
+        #  default alpha_ref value if not given
+        if self.alpha_ref_def is None:
+            warnings.warn(
+                Warning(
+                    """
+                    No default broadening in def file. Assigned alpha_ref = 0.07
+                    """
+                )
+            )
+            self.alpha_ref_def = 0.07
+
+        # load states
+        if self.states_file.with_suffix(".feather").exists():
+            states = pd.read_feather(self.states_file.with_suffix(".feather"))
+        else:
+            print(explanation)
+            states = exomolapi.read_states(self.states_file, dic_def)
+            states.to_feather(self.states_file.with_suffix(".feather"))
+        # load pf
+        pf = exomolapi.read_pf(self.pf_file)
+        self.gQT = pf["QT"].to_numpy()  # grid QT
+        self.T_gQT = pf["T"].to_numpy()  # T forgrid QT
+
+        # trans file(s)
+        print("Reading transition file")
+        if dic_def["numinf"] is None:
+            self.trans_file = self.path / pathlib.Path(molec + ".trans.bz2")
+            if not self.trans_file.exists():
+                self.download(molec, [".trans.bz2"])
+
+            if self.trans_file.with_suffix(".feather").exists():
+                trans = pd.read_feather(self.trans_file.with_suffix(".feather"))
+            else:
+                print(explanation)
+                trans = exomolapi.read_trans(self.trans_file)
+                trans.to_feather(self.trans_file.with_suffix(".feather"))
+
+            #!!TODO:restrict NOW the trans size to avoid useless overload of memory and CPU
+            # trans = trans[(trans['nu'] > self.nurange[0] - self.margin) & (trans['nu'] < self.nurange[1] + self.margin)]
+            # compute gup and elower
+            (
+                self._A,
+                self.nu_lines,
+                self._elower,
+                self._gpp,
+                self._jlower,
+                self._jupper,
+                self._quantumNumbers,
+            ) = exomolapi.pickup_gE(states, trans, dic_def)
+        else:
+            imin = (
+                np.searchsorted(dic_def["numinf"], nurange[0], side="right") - 1
+            )  # left side
+            imax = (
+                np.searchsorted(dic_def["numinf"], nurange[1], side="right") - 1
+            )  # left side
+            self.trans_file = []
+            for k, i in enumerate(range(imin, imax + 1)):
+                trans_file = self.path / pathlib.Path(
+                    molec + "__" + dic_def["numtag"][i] + ".trans.bz2"
+                )
+                if not trans_file.exists():
+                    self.download(
+                        molec, extension=[".trans.bz2"], numtag=dic_def["numtag"][i]
+                    )
+                if trans_file.with_suffix(".feather").exists():
+                    trans = pd.read_feather(trans_file.with_suffix(".feather"))
+                else:
+                    print(explanation)
+                    trans = exomolapi.read_trans(trans_file)
+                    trans.to_feather(trans_file.with_suffix(".feather"))
+                    #!!TODO:restrict NOW the trans size to avoid useless overload of memory and CPU
+                    # trans = trans[(trans['nu'] > self.nurange[0] - self.margin) & (trans['nu'] < self.nurange[1] + self.margin)]
+
+                self.trans_file.append(trans_file)
+
+                # compute gup and elower
+                if k == 0:
+                    (
+                        self._A,
+                        self.nu_lines,
+                        self._elower,
+                        self._gpp,
+                        self._jlower,
+                        self._jupper,
+                        self._quantumNumbers,
+                    ) = exomolapi.pickup_gE(states, trans, dic_def)
+                else:
+                    (
+                        Ax,
+                        nulx,
+                        elowerx,
+                        gppx,
+                        jlowerx,
+                        jupperx,
+                        _quantumNumbersx,
+                    ) = exomolapi.pickup_gE(states, trans, dic_def)
+                    self._A = np.hstack([self._A, Ax])
+                    self.nu_lines = np.hstack([self.nu_lines, nulx])
+                    self._elower = np.hstack([self._elower, elowerx])
+                    self._gpp = np.hstack([self._gpp, gppx])
+                    self._jlower = np.hstack([self._jlower, jlowerx])
+                    self._jupper = np.hstack([self._jupper, jupperx])
+                    self._quantumNumbers = np.hstack(
+                        [self._quantumNumbers, _quantumNumbersx]
+                    )
+
+        self.Tref = 296.0
+        self.QTref = np.array(self.QT_interp(self.Tref))
+
+        Ia = 1  #  TODO    Add isotope abundance
+
+        from radis.lbl.base import linestrength_from_Einstein  # TODO: move elsewhere
+
+        self.Sij0 = linestrength_from_Einstein(
+            A=self._A,
+            gu=self._gpp,
+            El=self._elower,
+            Ia=Ia,
+            nu=self.nu_lines,
+            Q=self.QTref,
+            T=self.Tref,
+        )
+
+        ### MASKING ###
+        mask = (
+            (self.nu_lines > self.nurange[0] - self.margin)
+            * (self.nu_lines < self.nurange[1] + self.margin)
+            * (self.Sij0 > self.crit)
+        )
+
+        self.masking(mask)
+
+    def masking(self, mask):
+        """applying mask and (re)generate jnp.arrays
+
+        Args:
+           mask: mask to be applied. self.mask is updated.
+
+        """
+        # TODO : replace with HDF5 masking ?
+
+        # numpy float 64 Do not convert them jnp array
+        self.nu_lines = self.nu_lines[mask]
+        self.Sij0 = self.Sij0[mask]
+        self._A = self._A[mask]
+        self._elower = self._elower[mask]
+        self._gpp = self._gpp[mask]
+        self._jlower = self._jlower[mask]
+        self._jupper = self._jupper[mask]
+
+        # jnp arraysgamma_
+        self.dev_nu_lines = np.array(self.nu_lines)
+        self.logsij0 = np.array(np.log(self.Sij0))
+        self.A = np.array(self._A)
+        # self.gamma_natural = gn(self.A)   #  Natural Broadening neglected in RADIS
+        self.elower = np.array(self._elower)
+        self.gpp = np.array(self._gpp)
+        self.jlower = np.array(self._jlower, dtype=int)
+        self.jupper = np.array(self._jupper, dtype=int)
+        ##Broadening parameters
+        self.set_broadening()
+
+        for key, array in self._quantumNumbers.items():
+            self._quantumNumbers[key] = array[mask]
+
+    def set_broadening(self, alpha_ref_def=None, n_Texp_def=None):
+        """setting broadening parameters
+
+        Args:
+           alpha_ref: set default alpha_ref and apply it. None=use self.alpha_ref_def
+           n_Texp_def: set default n_Texp and apply it. None=use self.n_Texp_def
+        """
+        if alpha_ref_def:
+            self.alpha_ref_def = alpha_ref_def
+        if n_Texp_def:
+            self.n_Texp_def = n_Texp_def
+
+        if self.broadf:
+            try:
+                print(".broad is used.")
+                bdat = exomolapi.read_broad(self.broad_file)
+                codelv = exomolapi.check_bdat(bdat)
+                print("Broadening code level=", codelv)
+                if codelv == "a0":
+                    j2alpha_ref, j2n_Texp = exomolapi.make_j2b(
+                        bdat,
+                        alpha_ref_default=self.alpha_ref_def,
+                        n_Texp_default=self.n_Texp_def,
+                        jlower_max=np.max(self._jlower),
+                    )
+                    self.alpha_ref = np.array(j2alpha_ref[self._jlower])
+                    self.n_Texp = np.array(j2n_Texp[self._jlower])
+                elif codelv == "a1":
+                    j2alpha_ref, j2n_Texp = exomolapi.make_j2b(
+                        bdat,
+                        alpha_ref_default=self.alpha_ref_def,
+                        n_Texp_default=self.n_Texp_def,
+                        jlower_max=np.max(self._jlower),
+                    )
+                    jj2alpha_ref, jj2n_Texp = exomolapi.make_jj2b(
+                        bdat,
+                        j2alpha_ref_def=j2alpha_ref,
+                        j2n_Texp_def=j2n_Texp,
+                        jupper_max=np.max(self._jupper),
+                    )
+                    self.alpha_ref = np.array(jj2alpha_ref[self._jlower, self._jupper])
+                    self.n_Texp = np.array(jj2n_Texp[self._jlower, self._jupper])
+            except:
+                print(
+                    "Warning: Cannot load .broad. The default broadening parameters are used."
+                )
+                self.alpha_ref = np.array(
+                    self.alpha_ref_def * np.ones_like(self._jlower)
+                )
+                self.n_Texp = np.array(self.n_Texp_def * np.ones_like(self._jlower))
+
+        else:
+            print("The default broadening parameters are used.")
+            self.alpha_ref = np.array(self.alpha_ref_def * np.ones_like(self._jlower))
+            self.n_Texp = np.array(self.n_Texp_def * np.ones_like(self._jlower))
+
+    def QT_interp(self, T):
+        """interpolated partition function
+
+        Args:
+           T: temperature
+
+        Returns:
+           Q(T) interpolated in jnp.array
+
+        """
+        return np.interp(T, self.T_gQT, self.gQT)
+
+    def qr_interp(self, T):
+        """interpolated partition function ratio
+
+        Args:
+           T: temperature
+
+        Returns:
+           qr(T)=Q(T)/Q(Tref) interpolated in jnp.array
+
+        """
+        return self.QT_interp(T) / self.QT_interp(self.Tref)
+
+    def download(self, molec, extension, numtag=None):
+        """Downloading Exomol files
+
+        Args:
+           molec: like "12C-16O__Li2015"
+           extension: extension list e.g. [".pf",".def",".trans.bz2",".states.bz2",".broad"]
+           numtag: number tag of transition file if exists. e.g. "11100-11200"
+
+        Note:
+           The download URL is written in exojax.utils.url.
+
+        """
+        import os
+        import urllib.request
+
+        tag = molec.split("__")
+        molname_simple = e2s(tag[0])
+
+        for ext in extension:
+            if ext == ".trans.bz2" and numtag is not None:
+                ext = "__" + numtag + ext
+
+            if ext == ".broad":
+                pfname_arr = [
+                    tag[0] + "__H2" + ext,
+                    tag[0] + "__He" + ext,
+                    tag[0] + "__air" + ext,
+                ]
+                url = EXOMOL_URL + molname_simple + "/" + tag[0] + "/"
+            else:
+                pfname_arr = [molec + ext]
+                url = EXOMOL_URL + molname_simple + "/" + tag[0] + "/" + tag[1] + "/"
+
+            for pfname in pfname_arr:
+                pfpath = url + pfname
+                os.makedirs(str(self.path), exist_ok=True)
+                print("Downloading " + pfpath)
+                try:
+                    urllib.request.urlretrieve(pfpath, str(self.path / pfname))
+                except:
+                    print("Error: Couldn't download " + ext + " file and save.")
+
+    def to_partition_function_tabulator(self):
+        from radis.levels.partfunc import PartFuncExoMol
+
+        return PartFuncExoMol(self.isotope_fullname, self.T_gQT, self.gQT)
+
+    def to_df(self):
+        """Export Line Database to a RADIS-friendly Pandas DataFrame"""
+
+        ##Broadening parameters
+        self.set_broadening()
+
+        # TODO : define RADIS database format somewhere else; with description of the column names.
+        df = pd.DataFrame(
+            {
+                "wav": self.nu_lines,
+                "int": self.Sij0,
+                "A": self._A,
+                "airbrd": self.alpha_ref,  # temperature dependance exponent for air broadening
+                # "selbrd": None,                  # self-temperature dependant exponent. Not implementedi in ExoMol
+                "El": self._elower,
+                # "Tdpair": None,
+                # "Pshft": None,
+                "ju": self._jupper,
+                "jl": self._jlower,
+                "gpp": self._gpp,
+                "Tdpair": self.n_Texp,  # temperature dependance exponent. Here we use Tdair; no Tdpsel. TODO
+            }
+        )
+        for key, array in self._quantumNumbers.items():
+            try:
+                df[key] = pd.to_numeric(array, errors="raise")
+            except:
+                df[key] = array
+                df[key] = df[key].astype("str")
+
+        from radis.db.classes import HITRAN_MOLECULES
+
+        if self.molecule in HITRAN_MOLECULES:
+            df.attrs["id"] = get_molecule_identifier(self.molecule)
+        df.attrs["molecule"] = self.molecule
+        df.attrs["iso"] = 1  # TODO : get isotope number while parsing ExoMol name
+
+        assert self.Tref == 296  # default in RADIS
+
+        return df
+
+
+if __name__ == "__main__":
+    # mdb=MdbExomol("/home/kawahara/exojax/data/CO/12C-16O/Li2015/")
+    # mdb=MdbExomol("/home/kawahara/exojax/data/CH4/12C-1H4/YT34to10/",nurange=[6050.0,6150.0])
+
+    # %% Test  by overriding Spectrumfactory's DataFrame df0
+    # mdb = MdbExomol(".database/H2O/1H2-16O/POKAZATEL", [4310.0, 4320.0], crit=1.0e-45)
+    # df = mdb.to_df()
+    # from radis import SpectrumFactory
+    # sf = SpectrumFactory(
+    #     4310,
+    #     4320,
+    #     molecule="H2O",
+    #     isotope="1",
+    # )
+    # sf.fetch_databank(
+    #     "hitran"
+    # )  # placeholder. Load lines (will be replaced), load partition function.
+    # sf.df0 = df  # override.
+    # s = sf.eq_spectrum(500, name="ExoMol")
+    # # sf.fetch_databank('hitran')  # placeholder. Load lines (will be replaced), load partition function.
+    # # s_hit = sf.eq_spectrum(500, name='HITRAN')
+
+    #%% Test by direct caclulation
+    from radis import SpectrumFactory
+
+    sf = SpectrumFactory(
+        4310,
+        4320,
+        molecule="H2O",
+        isotope="1",
+        verbose=3,
+    )
+    sf.fetch_databank("exomol")
+    s = sf.eq_spectrum(500, name="ExoMol")
+    s.plot()
+
+#    mask=mdb.A>1.e-42
+#    mdb.masking(mask)
+#    mdb=MdbExomol("/home/kawahara/exojax/data/exomol/NH3/14N-1H3/CoYuTe/",nurange=[6050.0,6150.0])
+#    mdb=MdbExomol("/home/kawahara/exojax/data/exomol/H2S/1H2-32S/AYT2/",nurange=[6050.0,6150.0])
+#    mdb=MdbExomol("/home/kawahara/exojax/data/exomol/FeH/56Fe-1H/MoLLIST/",nurange=[6050.0,6150.0])
+#    mdb=MdbExomol("/home/kawahara/exojax/data/exomol/NO/14N-16O/NOname/14N-16O__NOname")
diff --git a/radis/io/exomol_utils.py b/radis/io/exomol_utils.py
new file mode 100644
index 000000000..1ca69d4a9
--- /dev/null
+++ b/radis/io/exomol_utils.py
@@ -0,0 +1,67 @@
+"""Utils to read ExoMol
+
+Borrowed from the `Exojax <https://github.com/HajimeKawahara/exojax>`__
+code (which you should also have a look at !), by @HajimeKawahara, under MIT License.
+
+"""
+
+import re
+
+
+def e2s(molname_exact):
+    """convert the exact molname (used in ExoMol) to the simple molname
+
+    Args:
+       molname_exact: the exact molname
+
+    Returns:
+       simple molname
+
+    Examples:
+       >>> print(e2s("12C-1H4"))
+       >>> CH4
+       >>> print(e2s("23Na-16O-1H"))
+       >>> NaOH
+       >>> print(e2s("HeH_p"))
+       >>> HeH_p
+       >>> print(e2s("trans-31P2-1H-2H")) #not working
+       >>> Warning: Exact molname  trans-31P2-1H-2H cannot be converted to simple molname
+       >>> trans-31P2-1H-2H
+
+    """
+
+    try:
+        t = molname_exact.split("-")
+        molname_simple = ""
+        for ele in t:
+            alp = "".join(re.findall(r"\D", ele))
+            num0 = re.split("[A-Z]", ele)[1]
+            if num0.isdigit():
+                num = num0
+            else:
+                num = ""
+            molname_simple = molname_simple + alp + num
+        return molname_simple
+    except:
+        print(
+            "Warning: Exact molname ",
+            molname_exact,
+            "cannot be converted to simple molname",
+        )
+        return molname_exact
+
+
+if __name__ == "__main__":
+    print(e2s("12C-1H4"))
+    print(e2s("23Na-16O-1H"))
+    print(e2s("HeH_p"))
+    print(e2s("trans-31P2-1H-2H"))  # not working
+
+    # TODO : move to unitary tests of exomol_utils
+    assert e2s("12C-1H4") == "CH4"
+    assert e2s("23Na-16O-1H") == "NaOH"
+    assert e2s("HeH_p") == "HeH_p"
+    assert e2s("trans-31P2-1H-2H") == "trans-31P2-1H-2H"  # convert not working
+
+    assert e2s("12C-16O2") == "CO2"
+    assert e2s("13C-16O2") == "CO2"
diff --git a/radis/io/exomolapi.py b/radis/io/exomolapi.py
new file mode 100644
index 000000000..405faf283
--- /dev/null
+++ b/radis/io/exomolapi.py
@@ -0,0 +1,514 @@
+"""API for Exomol molecular database
+
+Borrowed from the `Exojax <https://github.com/HajimeKawahara/exojax>`__
+code (which you should also have a look at !), by @HajimeKawahara, under MIT License.
+
+"""
+import numpy as np
+import pandas as pd
+
+
+def read_def(deff):
+    """Exomol IO for a definition file
+
+    Parameters
+    ----------
+    deff: definition file
+
+    Returns
+    -------
+    temperature exponent n_Texp
+    broadening parameter alpha_ref
+    molecular mass
+    numinf: nu minimum for trans
+    numtag: tag for wavelength range
+
+    Note:
+       For some molecules, ExoMol provides multiple trans files. numinf and numtag are the ranges and identifiers for the multiple trans files.
+
+
+    """
+
+    dat = pd.read_csv(deff, sep="#", names=("VAL", "COMMENT"))
+    alpha_ref = None
+    # texp = None
+    molmasssw = False
+    n_Texp = None
+    ntransf = 1
+    maxnu = 0.0
+    quantum_labels = []
+    for i, com in enumerate(dat["COMMENT"]):
+        if "Default value of Lorentzian half-width" in com:
+            alpha_ref = float(dat["VAL"][i])
+        elif "Default value of temperature exponent" in com:
+            n_Texp = float(dat["VAL"][i])
+        elif "Element symbol 2" in com:
+            molmasssw = True
+        elif "No. of transition files" in com:
+            ntransf = int(dat["VAL"][i])
+        elif "Maximum wavenumber (in cm-1)" in com:
+            maxnu = float(dat["VAL"][i])
+            # maxnu=20000.0
+        elif molmasssw:
+            c = np.unique(dat["VAL"][i].strip(" ").split(" "))
+            c = np.array(c, dtype=np.float)
+            molmass = np.max(c)
+            molmasssw = False
+
+        elif "Lifetime availability" in com:
+            lifetime = dat["VAL"][i] == 1
+        elif "Lande g-factor availability" in com:
+            lande = dat["VAL"][i] == 1
+
+        elif "Quantum label" in com:
+            quantum_labels.append(dat["VAL"][i].strip(" "))
+
+        # SOME DEF FILES CONTAINS ERRORS. THESE ARE THE EXCEPTIONS
+        if deff.stem == "12C-16O2__UCL-4000":
+            ntransf = 20
+        if deff.stem == "14N-1H3__CoYuTe":
+            maxnu = 20000.0
+
+    if ntransf > 1:
+        dnufile = maxnu / ntransf
+        numinf = dnufile * np.array(range(ntransf + 1))
+        numtag = []
+        for i in range(len(numinf) - 1):
+            imin = "{:05}".format(int(numinf[i]))
+            imax = "{:05}".format(int(numinf[i + 1]))
+            numtag.append(imin + "-" + imax)
+    else:
+        numinf = None
+        numtag = ""
+
+    output = {
+        "n_Texp": n_Texp,
+        "alpha_ref": alpha_ref,
+        "molmass": molmass,
+        "numinf": numinf,
+        "numtag": numtag,
+        "quantum_labels": quantum_labels,  # array
+        "Landé": lande,  # bool
+        "lifetime": lifetime,  # bool
+    }
+    return output
+
+
+def read_pf(pff):
+    """Exomol IO for partition file
+
+    Note:
+        T=temperature
+        QT=partition function
+
+    Args:
+        pff: partition file
+    Returns:
+        partition data in pandas DataFrame
+
+    """
+    dat = pd.read_csv(pff, sep="\s+", names=("T", "QT"))
+    return dat
+
+
+def read_trans(transf):
+    """Exomol IO for a transition file
+    Note:
+        i_upper=Upper state counting number
+        i_lower=Lower state counting number
+        A=Einstein coefficient in s-1
+        nu_lines=transition wavenumber in cm-1
+        See Table 12 in https://arxiv.org/pdf/1603.05890.pdf
+
+    Args:
+        transf: transition file
+    Returns:
+        transition data in pandas DataFrame
+
+    """
+    try:
+        dat = pd.read_csv(
+            transf,
+            compression="bz2",
+            sep="\s+",
+            names=("i_upper", "i_lower", "A", "nu_lines"),
+        )
+    except:
+        dat = pd.read_csv(
+            transf, sep="\s+", names=("i_upper", "i_lower", "A", "nu_lines")
+        )
+
+    return dat
+
+
+def read_states(statesf, dic_def):
+    """Exomol IO for a state file
+    Note:
+        i=state counting number
+        E=state energy
+        g=state degeneracy
+        J=total angular momentum
+        See Table 11 in https://arxiv.org/pdf/1603.05890.pdf
+
+    Args:
+        statesf: state file
+        dic_def: Info from def file to read extra quantum numbers
+    Returns:
+        states data in pandas DataFrame
+
+    """
+    # we read first 4 columns for ("i", "E", "g", "J"),
+    # skip lifetime, skip Landé g-factor,
+    # read quantum numbers
+    N = np.size(dic_def["quantum_labels"])
+    quantum_labels = dic_def["quantum_labels"]
+
+    if dic_def["Landé"] and dic_def["lifetime"]:
+        usecol = np.concatenate((np.arange(4), 5 + 2 + np.arange(N)))
+        names = ("i", "E", "g", "J", "Lande", "lifetime") + tuple(quantum_labels)
+    elif dic_def["Landé"]:
+        usecol = np.concatenate((np.arange(4), 5 + 1 + np.arange(N)))
+        names = ("i", "E", "g", "J", "Lande") + tuple(quantum_labels)
+    elif dic_def["lifetime"]:
+        usecol = np.concatenate((np.arange(4), 5 + 1 + np.arange(N)))
+        names = ("i", "E", "g", "J", "lifetime") + tuple(quantum_labels)
+    else:  # no lifetime, nor landé according to the def file
+        usecol = np.concatenate((np.arange(4), 5 + np.arange(N)))
+        names = ("i", "E", "g", "J") + tuple(quantum_labels)
+
+    try:
+        dat = pd.read_csv(
+            statesf, compression="bz2", sep="\s+", usecols=usecol, names=names
+        )
+    except:  #!!!TODO What was the expected error?
+        dat = pd.read_csv(statesf, sep="\s+", usecols=usecol, names=names)
+
+    return dat
+
+
+def pickup_gE(states, trans, dic_def, trans_lines=False):
+    """extract g_upper (gup), E_lower (elower), and J_lower and J_upper from states DataFrame and insert them to transition DataFrame.
+
+    Args:
+       states: states pandas DataFrame
+       trans: transition pandas DataFrame
+       trans_lines: By default (False) we use nu_lines computed using the state file, i.e. E_upper - E_lower. If trans_nuline=True, we use the nu_lines in the transition file. Note that some trans files do not this info.
+       dic_def: Informations about additional quantum labels
+
+
+    Returns:
+       A, nu_lines, elower, gup, jlower, jupper, **quantum_labels
+
+    Note:
+       We first convert pandas DataFrame to ndarray. The state counting numbers in states DataFrame is used as indices of the new array for states (newstates). We remove the state count numbers as the column of newstate, i.e. newstates[:,k] k=0: E, 1: g, 2: J. Then, we can directly use the state counting numbers as mask.
+
+
+    """
+    ### Step 1. Essential quantum number for spectra
+    ndstates = states.to_numpy()[:, :4]  # the i, E, g, J are in the 4 first columns
+    ndtrans = trans.to_numpy()
+
+    iorig = np.array(ndstates[:, 0], dtype=int)
+    maxii = int(np.max(iorig) + 1)
+    newstates = np.zeros((maxii, np.shape(ndstates)[1] - 1), dtype=float)
+    newstates[iorig, :] = ndstates[:, 1:]
+
+    i_upper = np.array(ndtrans[:, 0], dtype=int)
+    i_lower = np.array(ndtrans[:, 1], dtype=int)
+
+    # use the state counting numbers (upper and lower) as masking.
+    elower = newstates[i_lower, 0]
+    eupper = newstates[i_upper, 0]
+    gup = newstates[i_upper, 1]
+    jlower = np.array(newstates[i_lower, 2], dtype=int)
+    jupper = np.array(newstates[i_upper, 2], dtype=int)
+
+    ### Step 2. Extra quantum numbers (e/f parity, vib and rot numbers)
+    #!!! Todo: Awfull management of memory, i know! @minou
+    n0 = (
+        4 + dic_def["Landé"] + dic_def["lifetime"]
+    )  # we need to shift if these quantities exist
+    ndstates_extra = states.to_numpy()[
+        :, n0:
+    ]  # the i, E, g, J are in the 4 first columns
+    a, b = np.shape(ndstates_extra)
+    dico_quantumNumbers = {}
+    if b != 0:
+        for index, label in enumerate(dic_def["quantum_labels"]):
+            if label == "K":  #!!!TODO: allow users to have other quantum numbers
+                dico_quantumNumbers["{}_l".format(label)] = ndstates_extra[
+                    i_lower, index
+                ]
+                # dico_quantumNumbers['{}_u'.format(label)] = ndstates_extra[index][i_upper]
+
+    A = ndtrans[:, 2]
+
+    if trans_lines:
+        nu_lines = ndtrans[:, 3]
+    else:
+        nu_lines = eupper - elower
+
+    # See Issue #16
+    # import matplotlib.pyplot as plt
+    # nu_lines_t=ndtrans[:,3]
+    # plt.plot(nu_lines_t-nu_lines,".",alpha=0.03)
+    # plt.ylabel("diff nu from trans and state (cm-1)")
+    # plt.xlabel("wavenuber (cm-1)")
+    # plt.savefig("nudiff.png", bbox_inches="tight", pad_inches=0.0)
+    # plt.show()
+
+    return A, nu_lines, elower, gup, jlower, jupper, dico_quantumNumbers
+
+
+# def pickup_gEslow(states, trans):
+#     """Slow version to extract g_upper (gup) and E_lower (elower) from states DataFrame and insert them to transition DataFrame.
+
+#     Note:
+#        This code is the (thousands times) slower version of pickup_gE. However, we did not remove this one just because this version is much easier to understand.
+
+#     Args:
+#        states: states pandas DataFrame
+#        trans: transition pandas DataFrame
+
+#     Returns:
+#        A, nu_lines, elower, gup
+
+#     """
+#     import tqdm
+
+#     E = states["E"].values
+#     g = states["g"].values
+
+#     # insert new columns in transition array
+#     trans["gup"] = 0
+#     trans["elower"] = 0.0
+
+#     for k, i in tqdm.tqdm(enumerate(states["i"])):
+#         # transition upper state
+#         mask_upper = trans["i_upper"] == i
+#         trans["gup"][mask_upper] = g[k]
+#         # transition lower state
+#         mask_lower = trans["i_lower"] == i
+#         trans["elower"][mask_lower] = E[k]
+
+#     A = trans["A"].to_numpy()
+#     nu_lines = trans["nu_lines"].to_numpy()
+#     elower = trans["elower"].to_numpy()
+#     gup = trans["gup"].to_numpy()
+#     return A, nu_lines, elower, gup
+
+
+# def read_broad(broadf):
+#     """Reading braodening file (.broad)
+
+#     Args:
+#        broadf: .broad file
+
+#     Return:
+#        broadening info in bdat form (pandas), defined by this instance.
+
+#     Note:
+#        See Table 16 in https://arxiv.org/pdf/1603.05890.pdf
+
+#     """
+#     bdat = pd.read_csv(
+#         broadf,
+#         sep="\s+",
+#         names=(
+#             "code",
+#             "alpha_ref",
+#             "n_Texp",
+#             "jlower",
+#             "jupper",
+#             "kalower",
+#             "kclower",
+#             "kaupper",
+#             "kcupper",
+#             "v1lower",
+#             "v2lower",
+#             "v3lower",
+#             "v1upper",
+#             "v2upper",
+#             "v3upper",
+#         ),
+#     )
+#     return bdat
+
+
+# def check_bdat(bdat):
+#     """cheking codes in .broad
+
+#     Args:
+#        bdat: exomol .broad data given by exomolapi.read_broad
+
+#     Returns:
+#        code level: None, a0, a1, other codes unavailable currently,
+
+#     """
+
+#     def checkcode(code):
+#         cmask = bdat["code"] == code
+#         if len(bdat["code"][cmask]) > 0:
+#             return True
+#         else:
+#             return False
+
+#     codelv = None
+#     for code in ["a0", "a1"]:
+#         if checkcode(code):
+#             codelv = code
+
+#     return codelv
+
+
+# def make_j2b(bdat, alpha_ref_default=0.07, n_Texp_default=0.5, jlower_max=None):
+#     """compute j2b (code a0, map from jlower to alpha_ref)
+
+#     Args:
+#        bdat: exomol .broad data given by exomolapi.read_broad
+#        alpha_ref_default: default value
+#        n_Texp_default: default value
+#        jlower_max: maximum number of jlower
+
+#     Returns:
+#        j2alpha_ref[jlower] provides alpha_ref for jlower
+#        j2n_Texp[jlower]  provides nT_exp for jlower
+
+#     """
+#     # a0
+#     cmask = bdat["code"] == "a0"
+#     jlower_arr = np.array(bdat["jlower"][cmask], dtype=int)
+#     alpha_ref_arr = np.array(bdat["alpha_ref"][cmask])
+#     n_Texp_arr = np.array(bdat["n_Texp"][cmask])
+
+#     if jlower_max is None:
+#         Nblower = np.max(jlower_arr) + 1
+#     else:
+#         Nblower = np.max([jlower_max, np.max(jlower_arr)]) + 1
+#     j2alpha_ref = np.ones(Nblower) * alpha_ref_default
+#     j2n_Texp = np.ones(Nblower) * n_Texp_default
+
+#     j2alpha_ref[jlower_arr] = alpha_ref_arr
+#     j2n_Texp[jlower_arr] = n_Texp_arr
+
+#     Ndef = Nblower - (np.max(jlower_arr) + 1)
+#     if Ndef > 0:
+#         print(
+#             "default broadening parameters are used for ",
+#             Ndef,
+#             " J lower states in ",
+#             Nblower,
+#             " states",
+#         )
+
+#     return j2alpha_ref, j2n_Texp
+
+
+# def make_jj2b(bdat, j2alpha_ref_def, j2n_Texp_def, jupper_max=None):
+#     """compute jj2b (code a1, map from (jlower, jupper) to alpha_ref and n_Texp)
+
+#     Args:
+#        bdat: exomol .broad data given by exomolapi.read_broad
+#        j2alpha_ref_def: default value from a0
+#        j2n_Texp_def: default value from a0
+#        jupper_max: maximum number of jupper
+
+#     Returns:
+#        jj2alpha_ref[jlower,jupper] provides alpha_ref for (jlower, jupper)
+#        jj2n_Texp[jlower,jupper]  provides nT_exp for (jlower, jupper)
+
+#     Note:
+#        The pair of (jlower, jupper) for which broadening parameters are not given, jj2XXX contains None.
+
+#     """
+#     # a1
+#     cmask = bdat["code"] == "a1"
+#     jlower_arr = np.array(bdat["jlower"][cmask], dtype=int)
+#     jupper_arr = np.array(bdat["jupper"][cmask], dtype=int)
+#     alpha_ref_arr = np.array(bdat["alpha_ref"][cmask])
+#     n_Texp_arr = np.array(bdat["n_Texp"][cmask])
+
+#     # Nblower = len(j2alpha_ref_def)
+
+#     if jupper_max is None:
+#         Nbupper = np.max(jupper_arr) + 1
+#     else:
+#         Nbupper = np.max([jupper_max, np.max(jupper_arr)]) + 1
+
+#     jj2alpha_ref = j2alpha_ref_def[:, np.newaxis] * np.ones(Nbupper)
+#     jj2n_Texp = j2n_Texp_def[:, np.newaxis] * np.ones(Nbupper)
+
+#     jj2alpha_ref[jlower_arr, jupper_arr] = alpha_ref_arr
+#     jj2n_Texp[jlower_arr, jupper_arr] = n_Texp_arr
+
+#     return jj2alpha_ref, jj2n_Texp
+
+
+if __name__ == "__main__":
+    pass
+    # import pathlib
+    # import sys
+    # import time
+
+    # deff = pathlib.Path(
+    #     "/home/kawahara/exojax/examples/luhman16/.database/CO2/12C-16O2/UCL-4000/12C-16O2__UCL-4000.def"
+    # )
+    # n_Texp, alpha_ref, molmass, numinf, numtag = read_def(deff)
+    # print(numtag)
+    # sys.exit()
+    # # various broad file
+    # #    broadf="/home/kawahara/exojax/data/broad/12C-16O__H2.broad"
+    # broadf = "/home/kawahara/exojax/data/broad/1H2-16O__H2.broad"
+    # bdat = read_broad(broadf)
+    # codelv = check_bdat(bdat)
+    # print(codelv)
+    # if codelv == "a0":
+    #     j2alpha_ref, j2n_Texp = make_j2b(bdat, jlower_max=100)
+    # elif codelv == "a1":
+    #     j2alpha_ref, j2n_Texp = make_j2b(bdat, jlower_max=100)
+    #     jj2alpha_ref, jj2n_Texp = make_jj2b(bdat, j2alpha_ref, j2n_Texp, jupper_max=100)
+    #     print(jj2alpha_ref[1, 2])
+    #     print(jj2alpha_ref[1, 15])
+
+    # sys.exit()
+    # # broad file
+    # broadf = "/home/kawahara/exojax/data/CO/12C-16O/12C-16O__H2.broad"
+    # bdat = read_broad(broadf)
+    # j2alpha_ref, j2n_Texp = make_j2b(bdat, jlower_max=100)
+
+    # # partition file
+    # pff = "/home/kawahara/exojax/data/exomol/CO/12C-16O/Li2015/12C-16O__Li2015.pf"
+    # dat = read_pf(pff)
+
+    # check = False
+    # if check:
+    #     print("Checking compution of Elower and gupper.")
+    # statesf = (
+    #     "/home/kawahara/exojax/data/exomol/CO/12C-16O/Li2015/12C-16O__Li2015.states.bz2"
+    # )
+    # states = read_states(statesf)
+    # transf = (
+    #     "/home/kawahara/exojax/data/exomol/CO/12C-16O/Li2015/12C-16O__Li2015.trans.bz2"
+    # )
+    # trans = read_trans(transf)
+
+    # ts = time.time()
+    # A, nu_lines, elower, gup, jlower, jupper = pickup_gE(states, trans)
+    # #    for i in range(0,len(A)):
+    # #        print(jlower[i],"-",jupper[i])
+    # te = time.time()
+
+    # tsx = time.time()
+    # if check:
+    #     A_s, nu_lines_s, elower_s, gup_s = pickup_gEslow(states, trans)
+    # tex = time.time()
+    # print(te - ts, "sec")
+    # if check:
+    #     print(tex - tsx, "sec for the slow version")
+    #     print("CHECKING DIFFERENCES...")
+    #     print(np.sum((A_s - A) ** 2))
+    #     print(np.sum((nu_lines_s - nu_lines) ** 2))
+    #     print(np.sum((elower_s - elower) ** 2))
+    #     print(np.sum((gup_s - gup) ** 2))
+
+    # # computing alpha_ref, n_Texp
diff --git a/radis/io/hitemp.py b/radis/io/hitemp.py
index 39f386897..415f08581 100644
--- a/radis/io/hitemp.py
+++ b/radis/io/hitemp.py
@@ -114,7 +114,8 @@ def fetch_hitemp(
     Other Parameters
     ----------------
     cache: bool, or ``'regen'``
-        if ``True``, use existing HDF5 file. If ``False`` or ``'regen'``, rebuild it.
+        if ``True``, use existing HDF5 file. If ``False`` or ``'regen'``, rebuild it..
+        Default ``True``.
     verbose: bool
     chunksize: int
         number of lines to process at a same time. Higher is usually faster
diff --git a/radis/io/query.py b/radis/io/query.py
index 160e6f1b2..bcb479a23 100644
--- a/radis/io/query.py
+++ b/radis/io/query.py
@@ -182,6 +182,7 @@ def fetch_astroquery(
     # Process response
 
     # Rename columns from Astroquery to RADIS format
+    # TODO : define RADIS database format somewhere else; with description of the column names.
     rename_columns = {
         "molec_id": "id",
         "local_iso_id": "iso",
diff --git a/radis/lbl/base.py b/radis/lbl/base.py
index 91fe0077a..4e0ba872e 100644
--- a/radis/lbl/base.py
+++ b/radis/lbl/base.py
@@ -520,6 +520,8 @@ def fillEvibl(r):
             df = df.groupby(by=["polyu", "wangu"]).apply(fillEvibu)
             # ~ 6.6 s (probably faster since neq==0.9.20 (radis<1.0)
             df = df.groupby(by=["polyl", "wangl"]).apply(fillEvibl)
+            # TODO : use map(dict) version (see _add_EvibErot_CDSD_pcN)
+
             #            except KeyError:
             #                import traceback
             #                traceback.print_exc()
@@ -641,21 +643,13 @@ def get_Evib_CDSD_pcN_1iso(df, iso):
 
             # Add lower state energy
             df_pcN = df.set_index(["polyl", "wangl", "rankl"])
-            #            for i in df.index:
-            #                df.loc[i, 'Evibl'] = energies.at[i, 'Evib']
-            # the map below is crazy fast compared to above loop
             df["Evibl"] = df_pcN.index.map(Evib_dict.get).values
             # Add upper state energy
             df_pcN = df.set_index(["polyu", "wangu", "ranku"])
-            #            for i in df.index:
-            #                df.loc[i, 'Evibu'] = energies.at[i, 'Evib']
-            # the map below is crazy fast compared to above loop
             df["Evibu"] = df_pcN.index.map(Evib_dict.get).values
 
             return df.loc[:, ["Evibl", "Evibu"]]
 
-        #        df = df.groupby('iso').apply(lambda x: get_Evib_CDSD_pcN_1iso(x, x.name))
-
         df["Evibl"] = np.nan
         df["Evibu"] = np.nan
 
@@ -696,13 +690,11 @@ def _add_EvibErot_CDSD_pcJN(self, df, calc_Evib_harmonic_anharmonic=False):
 
         Parameters
         ----------
-
         df: DataFrame
             list of transitions
 
         Other Parameters
         ----------------
-
         calc_Evib_harmonic_anharmonic: boolean
             if ``True``, calculate harmonic and anharmonic components of
             vibrational energies (for Treanor distributions)
@@ -743,7 +735,6 @@ def get_Evib_CDSD_pcJN_1iso(df, iso):
 
             Notes
             -----
-
             for devs:
 
             Unlike get_EvibErot_CDSD_pcN_1iso and get_EvibErot_CDSD_pc_1iso,
@@ -820,13 +811,11 @@ def _add_Evib123Erot_CDSD_pc(self, df, calc_Evib_harmonic_anharmonic=False):
 
         Parameters
         ----------
-
         df: DataFrame
             list of transitions
 
         Other Parameters
         ----------------
-
         calc_Evib_harmonic_anharmonic: boolean
             if ``True``, calculate harmonic and anharmonic components of
             vibrational energies (for Treanor distributions)
@@ -968,13 +957,11 @@ def _add_EvibErot_RADIS_cls1(self, df, calc_Evib_harmonic_anharmonic=False):
 
         Parameters
         ----------
-
         df: DataFrame
             list of transitions
 
         Other Parameters
         ----------------
-
         calc_Evib_harmonic_anharmonic: boolean
             if ``True``, calculate harmonic and anharmonic components of
             vibrational energies (for Treanor distributions)
@@ -1068,18 +1055,15 @@ def _add_Evib123Erot_RADIS_cls5(self, df):
 
         Parameters
         ----------
-
         df: DataFrame
             list of transitions
 
         Other Parameters
         ----------------
-
         calc_Evib_harmonic_anharmonic: boolean
             if ``True``, calculate harmonic and anharmonic components of
             vibrational energies (for Treanor distributions)
 
-
         """
         if __debug__:
             printdbg("called _add_Evib123Erot_RADIS_cls5()")
@@ -1110,7 +1094,6 @@ def get_Evib123_RADIS_cls5_1iso(df, iso):
 
             Notes
             -----
-
             We comb the Line Database with groups of same vibrational level,
             and fetch the corresponding vibrational energy from the Energy Level
             Database.
@@ -1204,7 +1187,6 @@ def _add_Evib123Erot_RADIS_cls5_harmonicanharmonic(self, df):
 
         Parameters
         ----------
-
         df: DataFrame
             list of transitions
         """
@@ -1519,50 +1501,30 @@ def _calc_degeneracies(self, df):
 
         See Also
         --------
-
         :func:`~radis.db.degeneracies.gs`, :func:`~radis.db.degeneracies.gi`
         """
         from radis.db.degeneracies import gi, gs
 
-        #        levelsfmt = self.params.levelsfmt
         dbformat = self.params.dbformat
 
         # Rotational + state specific degeneracy (J + isotope dependant)
         df["gju"] = 2 * df.ju + 1
         df["gjl"] = 2 * df.jl + 1
 
-        # multiple-isotopes in database
-        if "iso" in df:
-            isotope_set = df.iso.unique()
-            isotope = isotope_set[0]
-
-        else:
-            isotope = df.attrs["iso"]
-
-        if "id" in df and "iso" in df:
-            dgb = df.groupby(by=["id", "iso"])
-            for (id, iso), idx in dgb.indices.items():
-                _gs = gs(id, iso)
-                if isinstance(_gs, tuple):
-                    # Molecules that have alternating degeneracy.
-                    if id not in [2]:  # CO2
-                        raise NotImplementedError
-                    # normally we should find whether the rovibrational level is symmetric
-                    # or asymmetric. Here we just assume it's symmetric, because for
-                    # symmetric isotopes such as CO2(626), CO2 asymmetric levels
-                    # dont exist (gs=0) and they should not be in the line database.
-                    _gs = _gs[0]
-
-                dg = df.loc[idx]
-                _gi = gi(id, iso)
-                df.loc[idx, "grotu"] = dg.gju * _gs * _gi
-                df.loc[idx, "grotl"] = dg.gjl * _gs * _gi
-
-                if radis.DEBUG_MODE:
-                    assert (df.loc[idx, "id"] == id).all()
-                    assert (df.loc[idx, "iso"] == iso).all()
+        if "id" in df:
+            id_set = df.id.unique()
+            if len(id_set) > 0:
+                raise NotImplementedError("> 1 molecules in same DataFrame")
+            else:
+                self.warn(
+                    "There shouldn't be a Column 'id' with a unique value",
+                    "PerformanceWarning",
+                )
+            df.attrs["id"] = int(id_set)
+            if self.save_memory:
+                del df["id"]
 
-        elif "id" not in df and "iso" in df:
+        if "iso" in df:  # multiple isotopes in database
             id = df.attrs["id"]
             dgb = df.groupby(by=["iso"])
             for (iso), idx in dgb.indices.items():
@@ -1585,8 +1547,9 @@ def _calc_degeneracies(self, df):
                 if radis.DEBUG_MODE:
                     assert (df.loc[idx, "iso"] == iso).all()
 
-        elif ("id" not in df) and ("iso" not in df):
+        else:
             id = df.attrs["id"]
+            isotope = df.attrs["iso"]
             _gs = gs(id, isotope)
             if isinstance(_gs, tuple):
                 # Molecules that have alternating degeneracy.
@@ -1603,28 +1566,6 @@ def _calc_degeneracies(self, df):
             df.loc[:, "grotu"] = dg.gju * _gs * _gi
             df.loc[:, "grotl"] = dg.gjl * _gs * _gi
 
-        elif ("id" in df) and ("iso" not in df):
-            dgb = df.groupby(by=["id"])
-            for (id), idx in dgb.indices.items():
-                _gs = gs(id, isotope)
-                if isinstance(_gs, tuple):
-                    # Molecules that have alternating degeneracy.
-                    if id not in [2]:  # CO2
-                        raise NotImplementedError
-                    # normally we should find whether the rovibrational level is symmetric
-                    # or asymmetric. Here we just assume it's symmetric, because for
-                    # symmetric isotopes such as CO2(626), CO2 asymmetric levels
-                    # dont exist (gs=0) and they should not be in the line database.
-                    _gs = _gs[0]
-
-                dg = df.loc[idx]
-                _gi = gi(id, isotope)
-                df.loc[idx, "grotu"] = dg.gju * _gs * _gi
-                df.loc[idx, "grotl"] = dg.gjl * _gs * _gi
-
-                if radis.DEBUG_MODE:
-                    assert (df.loc[idx, "id"] == id).all()
-
         # %%
 
         if dbformat in [
@@ -1678,34 +1619,36 @@ def calc_weighted_trans_moment(self):
 
         if "id" in df:
             id_set = df.id.unique()
-            id = id_set[0]
-
-        else:
-            id = df.attrs["id"]
+            if len(id_set) > 0:
+                raise NotImplementedError("> 1 molecules in same DataFrame")
+            else:
+                self.warn(
+                    "There shouldn't be a Column 'id' with a unique value",
+                    "PerformanceWarning",
+                )
+            df.attrs["id"] = int(id_set)
+            if self.save_memory:
+                del df["id"]
+        molecule = get_molecule(df.attrs["id"])
 
         iso_set = self._get_isotope_list(self.input.molecule)  # df1.iso.unique()
 
-        Qref_dict = {}
-
-        if (("id" in df and len(id_set) == 1) or ("id" not in df and id)) and len(
-            iso_set
-        ) == 1:
+        if len(iso_set) == 1:
 
             # Shortcut if only 1 isotope. We attribute molar_mass & abundance
             # as attributes of the line database, instead of columns. Much
             # faster!
 
-            molecule = get_molecule(id)
             state = self.input.state
             parsum = self.get_partition_function_calculator(
                 molecule, iso_set[0], state
             )  # partition function
-            df.Qref = parsum.at(
+            df.attrs["Qref"] = parsum.at(
                 Tref, update_populations=False
             )  # stored as attribute, not column
             assert "Qref" not in df.columns
 
-            Qref = df.Qref
+            Qref = df.attrs["Qref"]
 
         else:
 
@@ -1713,7 +1656,7 @@ def calc_weighted_trans_moment(self):
             # still much faster than the groupby().apply() method (see radis<=0.9.19)
             # (tested + see https://stackoverflow.com/questions/44954514/efficient-way-to-conditionally-populate-elements-in-a-pandas-groupby-object-pos)
 
-            molecule = get_molecule(id)
+            Qref_dict = {}
 
             dgb = df.groupby(by=["iso"])
             for (iso), idx in dgb.indices.items():
@@ -1759,6 +1702,77 @@ def calc_weighted_trans_moment(self):
 
         return
 
+    def calc_reference_linestrength(self):
+        """ Calculate reference linestrength from Einstein coefficients"""
+
+        df = self.df0
+        Tref = self.input.Tref
+
+        self.profiler.start("calc_ref_linestrength", 2)
+
+        if "id" in df:
+            id_set = df.id.unique()
+            if len(id_set) > 0:
+                raise NotImplementedError("> 1 molecules in same DataFrame")
+            else:
+                self.warn(
+                    "There shouldn't be a Column 'id' with a unique value",
+                    "PerformanceWarning",
+                )
+            df.attrs["id"] = int(id_set)
+        molecule = get_molecule(df.attrs["id"])
+
+        iso_set = self._get_isotope_list(self.input.molecule)  # df1.iso.unique()
+
+        if len(iso_set) == 1:
+
+            # Shortcut if only 1 isotope. We attribute molar_mass & abundance
+            # as attributes of the line database, instead of columns. Much
+            # faster!
+
+            state = self.input.state
+            parsum = self.get_partition_function_calculator(
+                molecule, iso_set[0], state
+            )  # partition function
+            df.attrs["Qref"] = parsum.at(
+                Tref, update_populations=False
+            )  # stored as attribute, not column
+            assert "Qref" not in df.columns
+            Qref = df.attrs["Qref"]
+
+        else:
+
+            # normal method
+            # still much faster than the groupby().apply() method (see radis<=0.9.19)
+            # (tested + see https://stackoverflow.com/questions/44954514/efficient-way-to-conditionally-populate-elements-in-a-pandas-groupby-object-pos)
+
+            Qref_dict = {}
+
+            dgb = df.groupby(by=["iso"])
+            for (iso), idx in dgb.indices.items():
+                state = self.input.state
+                parsum = self.get_partition_function_calculator(
+                    molecule, iso, state
+                )  # partition function
+                Qref_dict[iso] = parsum.at(Tref, update_populations=False)
+                # ... note: do not update the populations here, so populations in the
+                # ... energy level list correspond to the one calculated for T and not Tref
+
+                if radis.DEBUG_MODE:
+                    if "id" in df:
+                        assert (df.loc[idx, "id"] == id).all()
+                    assert (df.loc[idx, "iso"] == iso).all()
+
+            Qref = df["iso"].map(Qref_dict)
+
+        S0 = linestrength_from_Einstein(
+            A=df.A, gu=df.gu, El=df.El, Ia=df.Ia, nu=df.wav, Q=Qref, T=Tref
+        )
+
+        self.profiler.stop("calc_ref_linestrength", "Calculated reference linestrength")
+
+        return S0
+
     def calc_einstein_coefficients(self):
         """Calculate :math:`A_{ul}`, :math:`B_{lu}`, :math:`B_{ul}` Einstein coefficients from weighted
         transition moments squared :math:`R_s^2`.
@@ -1769,7 +1783,6 @@ def calc_einstein_coefficients(self):
 
         Notes
         -----
-
         Einstein A coefficient already in database under df0.A
         Difference between df0.A and df0.Aul < 0.5%
 
@@ -1857,7 +1870,14 @@ def calc_lineshift(self):
 
         # Calculate
         air_pressure = self.input.pressure_mbar / 1013.25  # convert from mbar to atm
-        df["shiftwav"] = df.wav + (df.Pshft * air_pressure)
+        if "Pshft" in df.columns:
+            df["shiftwav"] = df.wav + (df.Pshft * air_pressure)
+        else:
+            self.warn(
+                "Pressure-shift coefficient not given in database: assumed 0 pressure shift",
+                "MissingPressureShiftWarning",
+            )
+            df["shiftwav"] = df.wav
 
         self.profiler.stop("calc_lineshift", "Calculated lineshift")
 
@@ -1943,49 +1963,40 @@ def _calc_Q(molecule, iso, state):
 
         if "id" in df1.columns:
             id_set = df1.id.unique()
-            id = id_set[0]
-
-        else:
-            id = df1.attrs["id"]
+            if len(id_set) > 1:
+                raise NotImplementedError(">1 molecule.")
+            else:
+                self.warn(
+                    "There shouldn't be a Column 'id' with a unique value",
+                    "PerformanceWarning",
+                )
+                df1.attrs["id"] = int(id_set)
 
         Qref_Qgas_ratio = {}
 
-        if ("id" not in df1 and id) or ("id" in df1 and len(id_set) == 1):
-            molecule = get_molecule(id)
-            iso_set = self._get_isotope_list(molecule)  # df1.iso.unique()
-
-            # Shortcut if only 1 molecule, 1 isotope. We attribute molar_mass & abundance
-            # as attributes of the line database, instead of columns. Much
-            # faster!
-
-            if len(iso_set) == 1:
-                Qref, Qgas = _calc_Q(molecule, iso_set[0], self.input.state)
-                Qref_Qgas_ratio[iso_set[0]] = Qref / Qgas
-
-            # Else, parse for all isotopes. Use df.map that is very fast
-            # Use the fact that isotopes are int, and thus can be considered as
-            # index in an array.
-            # ... in the following we exploit this to use the np.map function,
-            # ... which is amazingly fast
-            # ... Read https://stackoverflow.com/questions/49259580 to understand more
-
-            else:
+        if "molecule" in df1.attrs:
+            molecule = df1.attrs["molecule"]  # used for ExoMol, which has no HITRAN-id
+        else:
+            molecule = get_molecule(df1.attrs["id"])
+        iso_set = self._get_isotope_list(molecule)  # df1.iso.unique()
 
-                iso_arr = list(range(max(iso_set) + 1))
+        # Shortcut if only 1 molecule, 1 isotope. We attribute molar_mass & abundance
+        # as attributes of the line database, instead of columns. Much
+        # faster!
 
-                for iso in iso_arr:
-                    if iso in iso_set:
-                        Qref, Qgas = _calc_Q(molecule, iso, self.input.state)
-                        Qref_Qgas_ratio[iso] = Qref / Qgas
+        if len(iso_set) == 1:
+            Qref, Qgas = _calc_Q(molecule, iso_set[0], self.input.state)
+            Qref_Qgas_ratio[iso_set[0]] = Qref / Qgas
 
+        # Else, parse for all isotopes. Use df.map that is very fast
         else:
-            raise NotImplementedError(
-                ">1 molecule. Can use the np.take trick. Need to "
-                + "fallback to pandas.map(dict)"
-            )
-            # TODO: Implement. Read https://stackoverflow.com/a/51388828/5622825 to understand more
 
-        # Note on performance: few times faster than doing a groupby().apply()
+            iso_arr = list(range(max(iso_set) + 1))
+
+            for iso in iso_arr:
+                if iso in iso_set:
+                    Qref, Qgas = _calc_Q(molecule, iso, self.input.state)
+                    Qref_Qgas_ratio[iso] = Qref / Qgas
 
         # %% Calculate line strength at desired temperature
         # -------------------------------------------------
@@ -2303,22 +2314,24 @@ def calc_populations_noneq(
 
         if "id" in df:
             id_set = df.id.unique()
-            id = id_set[0]
-
-        else:
-            id = df.attrs["id"]
+            if len(id_set) > 0:
+                raise NotImplementedError("> 1 molecules in same DataFrame")
+            else:
+                self.warn(
+                    "There shouldn't be a Column 'id' with a unique value",
+                    "PerformanceWarning",
+                )
+                df.attrs["id"] = int(id_set)
+        molecule = get_molecule(df.attrs["id"])
 
         iso_set = self._get_isotope_list()  # df1.iso.unique()
-        if (("id" not in df and id) or ("id" in df and len(id_set) == 1)) and len(
-            iso_set
-        ) == 1:
+        if len(iso_set) == 1:
             iso = iso_set[0]
 
             # Shortcut if only 1 molecule, 1 isotope. We attribute molar_mass & abundance
             # as attributes of the line database, instead of columns. Much
             # faster!
 
-            molecule = get_molecule(id)
             state = self.input.state
             parsum = self.get_partition_function_calculator(molecule, iso, state)
             Q, Qvib, dfQrot = parsum.at_noneq(
@@ -2333,8 +2346,8 @@ def calc_populations_noneq(
             # ... make sure PartitionFunction above is calculated with the same
             # ... temperatures, rovibrational distributions and overpopulations
             # ... as the populations of active levels (somewhere below)
-            df.Qvib = Qvib
-            df.Q = Q
+            df.attrs["Qvib"] = Qvib
+            df.attrs["Q"] = Q
             assert "Qvib" not in df.columns
             assert "Q" not in df.columns
 
@@ -2352,110 +2365,48 @@ def calc_populations_noneq(
             dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
             df.loc[:, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
 
-            # ... note: the .map() version is about 3x faster than
-            # ... dg.groupby('viblvl_l').apply(lambda x: dfQrot_dict[x.name])
-
-        else:
-
-            # Normal version
-            # TODO: optimize by filling for first isotope first (see equilibrium case).
-            # TODO: use the np.take() trick if 1 molecule only
-
-            if "id" in df:
-                dgb = df.groupby(by=["id", "iso"])
-
-                for (id, iso), idx in dgb.indices.items():
-
-                    # Get partition function for all lines
-                    molecule = get_molecule(id)
-                    state = self.input.state
-                    parsum = self.get_partition_function_calculator(
-                        molecule, iso, state
-                    )
-                    Q, Qvib, dfQrot = parsum.at_noneq(
-                        Tvib,
-                        Trot,
-                        vib_distribution=vib_distribution,
-                        rot_distribution=rot_distribution,
-                        overpopulation=overpopulation,
-                        returnQvibQrot=True,
-                        update_populations=self.misc.export_populations,
-                    )
-
-                    # ... make sure PartitionFunction above is calculated with the same
-                    # ... temperatures, rovibrational distributions and overpopulations
-                    # ... as the populations of active levels (somewhere below)
-                    df.at[idx, "Qvib"] = Qvib
-                    df.at[idx, "Q"] = Q
-
-                    # reindexing to get a direct access to Qrot database
-                    # create the lookup dictionary
-                    # dfQrot index is already 'viblvl'
-                    dfQrot_dict = dict(list(zip(dfQrot.index, dfQrot.Qrot)))
-
-                    dg = df.loc[idx]
-
-                    # Add lower state Qrot
-                    dg_sorted = dg.set_index(["viblvl_l"], inplace=False)
-                    df.loc[idx, "Qrotl"] = dg_sorted.index.map(dfQrot_dict.get).values
-                    # Add upper state energy
-                    dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
-                    df.loc[idx, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
-
-                    # ... note: the .map() version is about 3x faster than
-                    # ... dg.groupby('viblvl_l').apply(lambda x: dfQrot_dict[x.name])
-
-                    if radis.DEBUG_MODE:
-                        assert (df.loc[idx, "id"] == id).all()
-                        assert (df.loc[idx, "iso"] == iso).all()
-
-            else:
-                dgb = df.groupby(by=["iso"])
-                id = df.attrs["id"]
-                molecule = get_molecule(id)
-
-                for (iso), idx in dgb.indices.items():
+        else:  #  multiple isotopes
+            dgb = df.groupby(by=["iso"])
+            for (iso), idx in dgb.indices.items():
 
-                    # Get partition function for all lines
-                    state = self.input.state
-                    parsum = self.get_partition_function_calculator(
-                        molecule, iso, state
-                    )
-                    Q, Qvib, dfQrot = parsum.at_noneq(
-                        Tvib,
-                        Trot,
-                        vib_distribution=vib_distribution,
-                        rot_distribution=rot_distribution,
-                        overpopulation=overpopulation,
-                        returnQvibQrot=True,
-                        update_populations=self.misc.export_populations,
-                    )
+                # Get partition function for all lines
+                state = self.input.state
+                parsum = self.get_partition_function_calculator(molecule, iso, state)
+                Q, Qvib, dfQrot = parsum.at_noneq(
+                    Tvib,
+                    Trot,
+                    vib_distribution=vib_distribution,
+                    rot_distribution=rot_distribution,
+                    overpopulation=overpopulation,
+                    returnQvibQrot=True,
+                    update_populations=self.misc.export_populations,
+                )
 
-                    # ... make sure PartitionFunction above is calculated with the same
-                    # ... temperatures, rovibrational distributions and overpopulations
-                    # ... as the populations of active levels (somewhere below)
-                    df.at[idx, "Qvib"] = Qvib
-                    df.at[idx, "Q"] = Q
+                # ... make sure PartitionFunction above is calculated with the same
+                # ... temperatures, rovibrational distributions and overpopulations
+                # ... as the populations of active levels (somewhere below)
+                df.at[idx, "Qvib"] = Qvib
+                df.at[idx, "Q"] = Q
 
-                    # reindexing to get a direct access to Qrot database
-                    # create the lookup dictionary
-                    # dfQrot index is already 'viblvl'
-                    dfQrot_dict = dict(list(zip(dfQrot.index, dfQrot.Qrot)))
+                # reindexing to get a direct access to Qrot database
+                # create the lookup dictionary
+                # dfQrot index is already 'viblvl'
+                dfQrot_dict = dict(list(zip(dfQrot.index, dfQrot.Qrot)))
 
-                    dg = df.loc[idx]
+                dg = df.loc[idx]
 
-                    # Add lower state Qrot
-                    dg_sorted = dg.set_index(["viblvl_l"], inplace=False)
-                    df.loc[idx, "Qrotl"] = dg_sorted.index.map(dfQrot_dict.get).values
-                    # Add upper state energy
-                    dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
-                    df.loc[idx, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
+                # Add lower state Qrot
+                dg_sorted = dg.set_index(["viblvl_l"], inplace=False)
+                df.loc[idx, "Qrotl"] = dg_sorted.index.map(dfQrot_dict.get).values
+                # Add upper state energy
+                dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
+                df.loc[idx, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
 
-                    # ... note: the .map() version is about 3x faster than
-                    # ... dg.groupby('viblvl_l').apply(lambda x: dfQrot_dict[x.name])
+                if radis.DEBUG_MODE:
+                    assert (df.loc[idx, "iso"] == iso).all()
 
-                    if radis.DEBUG_MODE:
-                        assert (df.loc[idx, "iso"] == iso).all()
+            Qvib = df.Qvib
+            Q = df.Q
 
         self.profiler.stop("part_function", "partition functions")
         self.profiler.start("population", 3)
@@ -2466,18 +2417,14 @@ def calc_populations_noneq(
         # ... vibrational distribution
         if vib_distribution == "boltzmann":
             # equation generated with @pytexit.py2tex > see docstrings.
-            df["nu_vib"] = df.gvibu / df.Qvib * exp(-hc_k * df.Evibu / Tvib)
-            df["nl_vib"] = df.gvibl / df.Qvib * exp(-hc_k * df.Evibl / Tvib)
+            df["nu_vib"] = df.gvibu / Qvib * exp(-hc_k * df.Evibu / Tvib)
+            df["nl_vib"] = df.gvibl / Qvib * exp(-hc_k * df.Evibl / Tvib)
         elif vib_distribution == "treanor":
             df["nu_vib"] = (
-                df.gvibu
-                / df.Qvib
-                * exp(-hc_k * (df.Evibu_h / Tvib + df.Evibu_a / Trot))
+                df.gvibu / Qvib * exp(-hc_k * (df.Evibu_h / Tvib + df.Evibu_a / Trot))
             )
             df["nl_vib"] = (
-                df.gvibl
-                / df.Qvib
-                * exp(-hc_k * (df.Evibl_h / Tvib + df.Evibl_a / Trot))
+                df.gvibl / Qvib * exp(-hc_k * (df.Evibl_h / Tvib + df.Evibl_a / Trot))
             )
         else:
             raise ValueError(
@@ -2501,8 +2448,8 @@ def calc_populations_noneq(
             )
 
         # ... Total
-        df["nu"] = df.nu_vib * df.nu_rot * (df.Qrotu * df.Qvib / df.Q)
-        df["nl"] = df.nl_vib * df.nl_rot * (df.Qrotl * df.Qvib / df.Q)
+        df["nu"] = df.nu_vib * df.nu_rot * (df.Qrotu * Qvib / Q)
+        df["nl"] = df.nl_vib * df.nl_rot * (df.Qrotl * Qvib / Q)
 
         if __debug__:
             assert "nu" in self.df1
@@ -3187,10 +3134,12 @@ def _cutoff_linestrength(self, cutoff=None):
         # ... @dev performance: quite long to select here, but I couldn't find a faster
         # ... alternative
         # TODO: remove useless columns in df1 to save memory
+        # Note @EP : with Vaex; the selection should be updated here
 
         # Ensures abundance, molar mass and partition functions are transfered
         # (needed if they are attributes and not isotopes)
-        transfer_metadata(df, self.df1, [k for k in df_metadata if hasattr(df, k)])
+        transfer_metadata(df, self.df1, [k for k in df_metadata if k in df.attrs])
+        # assert len(self.df1.attrs) > 0
 
         # Store number of lines cut (for information)
         self._Nlines_cutoff = Nlines_cutoff
@@ -3572,19 +3521,62 @@ def get_waverange(
     return wavenum_min, wavenum_max
 
 
+def linestrength_from_Einstein(A, gu, El, Ia, nu, Q, T):
+    """Calculate linestrength at temperature ``T`` from Einstein coefficients.
+
+    Parameters
+    ----------
+    A : float, s-1
+        Einstein emission coefficients
+    gu : int
+        upper state degeneracy
+    El : float, cm-1
+        lower state energy
+    Ia : float
+        isotope abundance
+    nu : cm-1
+        transition wavenumber
+    Q : float
+        partition function at temperature ``T``
+    T : float
+        temperature
+
+    Returns
+    -------
+    S : float
+        linestrength at temperature ``T``.
+
+    References
+    ----------
+
+    .. math::
+        S=\frac{A gu \operatorname{exp}\left(\frac{-hc_k El}{T}\right) \left(1-\operatorname{exp}\left(\frac{-hc_k \nu}{T}\right)\right)}{8\pi c_{CGS} {\nu}^2 Q}
+
+    Combine Eq.(A.5), (A.8) in [Rothman-1998]_
+
+
+    """
+
+    return (
+        A
+        * gu
+        * np.exp(-hc_k * El / T)
+        * (1 - np.exp(-hc_k * nu / T))
+        / (8 * np.pi * c_CGS * nu ** 2 * Q)
+    )
+
+
 def get_abundance(df):
 
     """Returns the abundance
 
     Parameters
     ----------
-
     df: dataframe
 
     Returns
     -------
-
-    The abundance of all the isotopes in the dataframe
+    dict: The abundance of all the isotopes in the dataframe
     """
 
     molpar = MolParams()
diff --git a/radis/lbl/broadening.py b/radis/lbl/broadening.py
index abcb25d42..5017e7eac 100644
--- a/radis/lbl/broadening.py
+++ b/radis/lbl/broadening.py
@@ -69,7 +69,6 @@
 from numpy import pi, sin, sqrt, trapz, zeros_like
 from scipy.signal import oaconvolve
 
-from radis.db.classes import get_molecule_identifier
 from radis.db.molparam import MolParams
 from radis.lbl.base import BaseFactory
 from radis.misc.arrays import add_at, numpy_add_at
@@ -262,8 +261,9 @@ def pressure_broadening_HWHM(
         half-width half max coefficient (HWHM ) for resonant (self) broadening with air
     Tdpair: array like     [length N]
         temperature dependance coefficient for collisional broadening with air
-    Tdpsel: array like     [length N]
-        temperature dependance coefficient for resonant (self) broadening
+    Tdpsel: array like, optional    [length N]
+        temperature dependance coefficient for resonant (self) broadening.
+        If ``None``, use ``Tdpair``.
     pressure_atm: float  [atm]
         pressure in atmosphere (warning, not bar!)
     mole_fraction: float    [0-1]
@@ -276,11 +276,20 @@ def pressure_broadening_HWHM(
 
     Returns
     -------
-    lineshape: pandas Series        [shape N x W]
-        Lorentzian half-width at half-maximum (FWHM) for each line profile
+    pandas Series        [shape N]
+        Lorentzian half-width at half-maximum (HWHM) for each line profile
 
     References
     ----------
+
+    ..math::
+
+        \\gamma_{lb}={\\left(\\frac{T_{ref}}{T_{gas}}\\right)}^{n_{air}} \\gamma_{air} P \\left(1-x\\right)+{\\left(\\frac{T_{ref}}{T_{gas}}\\right)}^{n_{self}} \\gamma_{self} P x
+
+    With :math:`n_{air}, n_{self}` the temperature dependance coefficients
+    ``Tdpair, Tdpsel`` ; :math:`\\gamma_{air}, \\gamma_{self}` the air and resonant
+    broadening ``airbrd, selbrd``, :math:`x` the ``mole_fraction``.
+
     .. [1] `Rothman 1998 (HITRAN 1996) eq (A.14) <https://www.sciencedirect.com/science/article/pii/S0022407398000788>`_
 
     See Also
@@ -802,14 +811,6 @@ def _calc_broadening_HWHM(self):
         # Init variables
         df = self.df1
 
-        isotope_set = list(self.input.isotope)
-        if len(isotope_set) == 1:
-            df.attrs["iso"] = int(isotope_set[0])
-        else:
-            df.attrs["iso"] = isotope_set
-        M = get_molecule_identifier(self.input.molecule)
-        df.attrs["id"] = M
-
         if len(df) == 0:
             return  # no lines
 
@@ -1003,21 +1004,31 @@ def _add_collisional_broadening_HWHM(
         .. [1] `Rothman 1998 (HITRAN 1996) eq (A.12) <https://www.sciencedirect.com/science/article/pii/S0022407398000788>`_
         """
 
-        # Check self broadening is here
+        # Check self broadening temperature-dependance coefficient is here
         if not "Tdpsel" in list(df.keys()):
             self.warn(
-                "Self-broadening temperature coefficient Tdpsel not given in database: used Tdpair instead",
-                "MissingSelfBroadeningWarning",
+                "Self-broadening temperature coefficient `Tdpsel` not given in database: used `Tdpair` instead",
+                "MissingSelfBroadeningTdepWarning",
                 level=2,  # only appear if verbose>=2
             )
-            Tdpsel = None
+            Tdpsel = None  # will be corrected in pressure_broadening_HWHM()
         else:
             Tdpsel = df.Tdpsel
 
+        # Check self broadening is here
+        if not "selbrd" in list(df.keys()):
+            self.warn(
+                "Self-broadening reference width `selbrd` not given in database: used air broadening reference width `airbrd` instead",
+                "MissingSelfBroadeningWarning",
+            )
+            selbrd = df.airbrd
+        else:
+            selbrd = df.selbrd
+
         # Calculate broadening HWHM
         wl = pressure_broadening_HWHM(
             df.airbrd,
-            df.selbrd,
+            selbrd,
             df.Tdpair,
             Tdpsel,
             pressure_atm,
@@ -2840,43 +2851,43 @@ def rough_sum_on_grid():
 
 
 def get_molar_mass(df):
-
-    """Returns molar mass
+    """Returns molar mass.
 
     Parameters
     ----------
-
     df: dataframe
 
     Returns
     -------
-
     The molar mass of all the isotopes in the dataframe
     """
-
     molpar = MolParams()
 
     if "id" in df.columns:
-        id = df.id.unique()
-        id = id[0]
-
-    else:
+        raise NotImplementedError(">1 molecule")
+    elif "id" in df.attrs:
         id = df.attrs["id"]
+    else:
+        # HARDCODED molar mass; for WIP ExoMol implementation, until MolParams
+        # is an attribute and can be updated with definitions from ExoMol.
+        # https://github.com/radis/radis/issues/321
+        HARCODED_MOLAR_MASS = {"SiO": {1: 43.971842}}
+        try:
+            return HARCODED_MOLAR_MASS[df.attrs["molecule"]][df.attrs["iso"]]
+        except KeyError:
+            raise NotImplementedError
 
     if "iso" in df.columns:
         iso_set = df.iso.unique()
-        molar_mass = {}
+        molar_mass_dict = {}
         for iso in iso_set:
-            molar_mass[iso] = molpar.get(id, iso, "mol_mass")
-
-        req = df["iso"].map(molar_mass)
+            molar_mass_dict[iso] = molpar.get(id, iso, "mol_mass")
+        molar_mass = df["iso"].map(molar_mass_dict)
     else:
         iso = df.attrs["iso"]
         molar_mass = molpar.get(id, iso, "mol_mass")
 
-        req = molar_mass
-
-    return req
+    return molar_mass
 
     #
 
diff --git a/radis/lbl/calc.py b/radis/lbl/calc.py
index 48afca7ac..d2bdb483a 100644
--- a/radis/lbl/calc.py
+++ b/radis/lbl/calc.py
@@ -564,20 +564,34 @@ def _is_at_equilibrium():
         export_lines=export_lines,
         **kwargs
     )
-    if databank in [
-        "fetch",
-        "hitran",
-        "hitemp",
-    ]:  # mode to get databank without relying on  Line databases
+    if databank in ["fetch", "hitran", "hitemp", "exomol",] or (
+        isinstance(databank, tuple) and databank[0] == "exomol"
+    ):  # mode to get databank without relying on  Line databases
         # Line database :
         if databank in ["fetch", "hitran"]:
-            conditions = {"source": "hitran"}
+            conditions = {
+                "source": "hitran",
+                "parfuncfmt": "hapi",  # use HAPI (TIPS) partition functions for equilibrium
+            }
         elif databank in ["hitemp"]:
-            conditions = {"source": "hitemp"}
+            conditions = {
+                "source": "hitemp",
+                "parfuncfmt": "hapi",  # use HAPI (TIPS) partition functions for equilibrium}
+            }
+        elif databank in ["exomol"]:
+            conditions = {
+                "source": "exomol",
+                "parfuncfmt": "exomol",  # download & use Exo partition functions for equilibrium}
+            }
+        elif isinstance(databank, tuple) and databank[0] == "exomol":
+            conditions = {
+                "source": "exomol",
+                "exomol_database": databank[1],
+                "parfuncfmt": "exomol",  # download & use Exo partition functions for equilibrium}
+            }
         # Partition functions :
         conditions.update(
             **{
-                "parfuncfmt": "hapi",  # use HAPI (TIPS) partition functions for equilibrium
                 "levelsfmt": None,  # no need to load energies by default
                 "db_use_cached": use_cached,
             }
@@ -657,7 +671,7 @@ def _is_at_equilibrium():
                 path_length=path_length,
                 name=name,
             )
-        else:
+        elif mode == "gpu":
             s = sf.eq_spectrum_gpu(
                 Tgas=Tgas,
                 mole_fraction=mole_fraction,
@@ -665,7 +679,11 @@ def _is_at_equilibrium():
                 path_length=path_length,
                 name=name,
             )
+        else:
+            raise ValueError(mode)
     else:
+        if mode != "cpu":
+            raise NotImplementedError(mode)
         s = sf.non_eq_spectrum(
             Tvib=Tvib,
             Trot=Trot,
diff --git a/radis/lbl/factory.py b/radis/lbl/factory.py
index 483112f25..52f360193 100644
--- a/radis/lbl/factory.py
+++ b/radis/lbl/factory.py
@@ -449,6 +449,7 @@ def __init__(
                     + "\n- under nonequilibrium: {0}".format(
                         MOLECULES_LIST_NONEQUILIBRIUM
                     )
+                    + "\n\nNote that RADIS now has ExoMol support, but not all ExoMol molecules are referenced in RADIS. If a molecule is available in ExoMol but does not appear in RADIS yet, please contact the RADIS team or write on https://github.com/radis/radis/issues/319"
                 )
 
         # Store isotope identifier in str format (list wont work in database queries)
@@ -1773,8 +1774,8 @@ def print_perf_profile(self, number_format="{:.3f}", precision=16):
 
         Examples
         --------
-
         ::
+
             sf.print_perf_profile()
 
             # output >>
@@ -1838,6 +1839,7 @@ def generate_perf_profile(self):
         Examples
         --------
         ::
+
             sf = SpectrumFactory(...)
             sf.eq_spectrum(...)
             sf.generate_perf_profile()
diff --git a/radis/lbl/loader.py b/radis/lbl/loader.py
index edcd28cfb..7da05200f 100644
--- a/radis/lbl/loader.py
+++ b/radis/lbl/loader.py
@@ -61,6 +61,7 @@
 from radis.db.molecules import getMolecule
 from radis.io.cache_files import cache_file_name
 from radis.io.cdsd import cdsd2df
+from radis.io.exomol import fetch_exomol
 from radis.io.hdf5 import hdf2df
 from radis.io.hitemp import fetch_hitemp
 from radis.io.hitran import hit2df, parse_global_quanta, parse_local_quanta
@@ -98,12 +99,14 @@
     "hdf5-radisdb",
 ]
 """list: Known formats for Line Databases:
+
 - ``'hitran'`` : [HITRAN-2016]_ original .par format
 - ``'hitemp'`` : [HITEMP-2010]_ original format (same format as 'hitran')
 - ``'cdsd-hitemp'`` : CDSD-HITEMP original format (CO2 only, same lines as HITEMP-2010)
 - ``'cdsd-4000'`` : [CDSD-4000]_ original format (CO2 only)
 - ``'hitemp-radisdb'`` : HITEMP under RADISDB format (pytables-HDF5 with RADIS column names).
 - ``'hdf5-radisdb'`` : arbitrary HDF5 file with RADIS column names.
+
 To install all databases manually see the :ref:`Configuration file <label_lbl_config_file>`
 and the :ref:`list of databases <label_line_databases>` .
 
@@ -114,6 +117,7 @@
 
 KNOWN_LVLFORMAT = ["radis", "cdsd-pc", "cdsd-pcN", "cdsd-hamil", None]
 """list: Known formats for Energy Level Databases (used in non-equilibrium calculations):
+
 - ``'radis'``: energies calculated with Dunham expansions by
     :class:`~radis.levels.partfunc.PartFunc_Dunham`
 - ``'cdsd-pc'``: energies read from precomputed CDSD energies for CO2, with
@@ -456,7 +460,7 @@ def format_paths(s):
     return s
 
 
-df_metadata = ["Ia", "molar_mass", "Qref", "Qvib", "Q"]
+df_metadata = ["molecule", "iso", "id", "Ia", "molar_mass", "Qref", "Qvib", "Q"]
 """ list: metadata of line DataFrames :py:attr:`~radis.lbl.loader.DatabankLoader.df0`,
 :py:attr:`~radis.lbl.loader.DatabankLoader.df1`.
 @dev: when having only 1 molecule, 1 isotope, these parameters are
@@ -603,8 +607,8 @@ def __init__(self):
 
     def init_databank(self, *args, **kwargs):
         """Method to init databank parameters but only load them when needed.
-        Databank is reloaded by
-        :meth:`~radis.lbl.loader.DatabankLoader._check_line_databank`
+
+        Databank is reloaded by :py:meth:`~radis.lbl.loader.DatabankLoader._check_line_databank`
         Same inputs Parameters as :meth:`~radis.lbl.loader.DatabankLoader.load_databank`:
 
         Parameters
@@ -726,6 +730,7 @@ def init_databank(self, *args, **kwargs):
     def fetch_databank(
         self,
         source="hitran",
+        exomol_database=None,
         parfunc=None,
         parfuncfmt="hapi",
         levels=None,
@@ -742,17 +747,23 @@ def fetch_databank(
 
         Parameters
         ----------
-        source: ``'hitran'``, ``'hitemp'``
-            [Download database lines from the latest HITRAN (see [HITRAN-2016]_)
-            or HITEMP version (see [HITEMP-2010]_  )]
+        source: ``'hitran'``, ``'hitemp'``, ``'exomol'``
+            [Download database lines from the latest HITRAN (see [HITRAN-2016]_),
+            HITEMP (see [HITEMP-2010]_  )] or EXOMOL see [ExoMol-2020]_  ) databases.
+        exomol_database: None
+            if fetching from ''`exomol`'', choose which database to use. Keep
+            ``None`` to use the recommended one. See all available databases
+            with :py:func:`radis.io.exomol.get_exomol_database_list`
+
 
         Other Parameters
         ----------------
-        parfuncfmt: ``'cdsd'``, ``'hapi'``, or any of :data:`~radis.lbl.loader.KNOWN_PARFUNCFORMAT`
+        parfuncfmt: ``'cdsd'``, ``'hapi'``, ``'exomol'``, or any of :data:`~radis.lbl.loader.KNOWN_PARFUNCFORMAT`
             format to read tabulated partition function file. If ``hapi``, then
             HAPI (HITRAN Python interface) [2]_ is used to retrieve them (valid if
             your database is HITRAN data). HAPI is embedded into RADIS. Check the
-            version. If partfuncfmt is None then ``hapi`` is used. Default ``hapi``.
+            version. If partfuncfmt is None then ``hapi`` is used. If ``'exomol'``
+            then partition functions are downloaded from ExoMol. Default ``hapi``.
         parfunc: filename or None
             path to tabulated partition function to use.
             If `parfuncfmt` is `hapi` then `parfunc` should be the link to the
@@ -814,7 +825,7 @@ def fetch_databank(
                 )
             )
             source = "hitran"
-        if source not in ["hitran", "hitemp"]:
+        if source not in ["hitran", "hitemp", "exomol"]:
             raise NotImplementedError("source: {0}".format(source))
         if source == "hitran":
             dbformat = "hitran"
@@ -822,6 +833,14 @@ def fetch_databank(
             dbformat = (
                 "hitemp-radisdb"  # downloaded in RADIS local databases ~/.radisdb
             )
+        elif source == "exomol":
+            dbformat = "exomol-radisdb"  # downloaded in RADIS local databases ~/.radisdb  # Note @EP : still WIP.
+        if exomol_database != None:
+            assert source == "exomol"
+        if [parfuncfmt, source].count("exomol") == 1:
+            raise NotImplementedError(
+                "ExoMol partition functions must be used with ExoMol database - and vice-versa"
+            )
 
         # Get inputs
         molecule = self.input.molecule
@@ -858,7 +877,7 @@ def fetch_databank(
             # Query one isotope at a time
             if isotope == "all":
                 raise ValueError(
-                    "Please define isotope explicitely (cannot use 'all' with fetch_databank)"
+                    "Please define isotope explicitely (cannot use 'all' with fetch_databank('hitran'))"
                 )
             isotope_list = self._get_isotope_list()
 
@@ -882,10 +901,15 @@ def fetch_databank(
                     f"{molecule} has no lines on range "
                     + "{0:.2f}-{1:.2f} cm-1".format(wavenum_min, wavenum_max)
                 )
-            else:
+            if len(frames) > 1:
+                # Note @dev : may be faster/less memory hungry to keep lines separated for each isotope. TODO : test both versions
+                for df in frames:
+                    assert "iso" in df.columns
                 df = pd.concat(frames, ignore_index=True)  # reindex
-
+            else:
+                df = frames[0]
             self.params.dbpath = "fetched from hitran"
+
         elif source == "hitemp":
             # Download, setup local databases, and fetch (use existing if possible)
 
@@ -911,6 +935,57 @@ def fetch_databank(
                     [str(k) for k in self._get_isotope_list(df=df)]
                 )
 
+        elif source == "exomol":
+            # Download, setup local databases, and fetch (use existing if possible)
+
+            if isotope == "all":
+                raise ValueError(
+                    "Please define isotope explicitely (cannot use 'all' with fetch_databank('exomol'))"
+                )
+            isotope_list = self._get_isotope_list()
+
+            local_paths = []
+            frames = []  # lines for all isotopes
+            partition_function_exomol = {
+                molecule: {}
+            }  # partition function tabulators for all isotpes
+            for iso in isotope_list:
+                df, local_path, Z_exomol = fetch_exomol(
+                    molecule,
+                    database=exomol_database,
+                    isotope=iso,
+                    load_wavenum_min=wavenum_min,
+                    load_wavenum_max=wavenum_max,
+                    cache=db_use_cached,
+                    verbose=self.verbose,
+                    return_local_path=True,
+                    return_partition_function=True,
+                )
+                # @dev refactor : have a DatabaseClass from which we load lines and partition functions
+                if len(df) > 0:
+                    frames.append(df)
+                local_paths.append(local_path)
+                partition_function_exomol[molecule][iso] = Z_exomol
+
+            # Merge
+            if frames == []:
+                raise EmptyDatabaseError(
+                    f"{molecule} has no lines on range "
+                    + "{0:.2f}-{1:.2f} cm-1".format(wavenum_min, wavenum_max)
+                )
+            if len(frames) > 1:
+                # Note @dev : may be faster/less memory hungry to keep lines separated for each isotope. TODO : test both versions
+                for df in frames:
+                    assert "iso" in df.columns
+                df = pd.concat(frames, ignore_index=True)  # reindex
+                self.params.dbpath = local_paths
+            else:
+                df = frames[0]
+                self.params.dbpath = local_paths[0]
+
+        else:
+            raise NotImplementedError("source: {0}".format(source))
+
         if len(df) == 0:
             raise EmptyDatabaseError(
                 f"{molecule} has no lines on range "
@@ -939,7 +1014,14 @@ def fetch_databank(
         # %% Init Partition functions (with energies)
         # ------------
 
-        self._init_equilibrium_partition_functions(parfunc, parfuncfmt)
+        if parfuncfmt == "exomol":
+            self._init_equilibrium_partition_functions(
+                parfunc,
+                parfuncfmt,
+                predefined_partition_functions=partition_function_exomol,
+            )
+        else:
+            self._init_equilibrium_partition_functions(parfunc, parfuncfmt)
 
         # If energy levels are given, initialize the partition function calculator
         # (necessary for non-equilibrium). If levelsfmt == 'radis' then energies
@@ -959,21 +1041,7 @@ def fetch_databank(
                     + "in fetch_databank"
                 )
 
-        id_set = df.id.unique()
-
-        if len(id_set) != 1:  # only 1 molecule supported ftm
-            raise NotImplementedError(
-                "Only 1 molecule at a time is currently supported "
-                + "in SpectrumFactory. Use radis.calc_spectrum, which "
-                + "calculates them independently then use MergeSlabs"
-            )
-
-        isotope_set = df.iso.unique()
-
-        if len(isotope_set) == 1:
-            df.drop("iso", axis=1, inplace=True)
-            df_metadata.append("iso")
-            df.attrs["iso"] = isotope_set[0]
+        self._remove_unecessary_columns(df)
 
         return
 
@@ -1427,12 +1495,14 @@ def init_database(
     # _reload_databank
     # _check_line_databank
     # _retrieve_from_database
-    # _get_partition_function_interpolator
-    # _get_partition_function_calculator
+    # get_partition_function_interpolator
+    # get_partition_function_calculator
     #
     # =========================================================================
 
-    def _init_equilibrium_partition_functions(self, parfunc, parfuncfmt):
+    def _init_equilibrium_partition_functions(
+        self, parfunc, parfuncfmt, predefined_partition_functions={}
+    ):
         """Initializes equilibrium partition functions in ``self.parsum_tab``
 
         Parameters
@@ -1446,6 +1516,12 @@ def _init_equilibrium_partition_functions(self, parfunc, parfuncfmt):
             path to tabulated partition function to use.
             If ``parfuncfmt`` is ``hapi`` then ``parfunc`` should be the link to the
             hapi.py file. If not given, then the hapi.py embedded in RADIS is used (check version)
+
+        Other Parameters
+        ----------------
+        predefined_partition_functions: dict
+            ::
+                {molecule: {isotope: PartitionFunctionTabulator object}}
         """
 
         # Let's get the tabulated partition function (to calculate eq spectra)
@@ -1455,7 +1531,11 @@ def _init_equilibrium_partition_functions(self, parfunc, parfuncfmt):
         for iso in self._get_isotope_list():
             self.parsum_tab[molecule][iso] = {}
             ParsumTab = self._build_partition_function_interpolator(
-                parfunc, parfuncfmt, self.input.molecule, isotope=iso
+                parfunc,
+                parfuncfmt,
+                self.input.molecule,
+                isotope=iso,
+                predefined_partition_functions=predefined_partition_functions,
             )
             self.parsum_tab[molecule][iso][state] = ParsumTab
 
@@ -1848,27 +1928,43 @@ def load_and_concat(files):
                 )
             )
 
-        id_set = df.id.unique()
+        self._remove_unecessary_columns(df)
 
-        if len(id_set) != 1:  # only 1 molecule supported ftm
-            raise NotImplementedError(
-                "Only 1 molecule at a time is currently supported "
-                + "in SpectrumFactory. Use radis.calc_spectrum, which "
-                + "calculates them independently then use MergeSlabs"
-            )
+        return df
+
+    def _remove_unecessary_columns(self, df):
+        """Remove unecessary columns and add values as attributes
+
+        Returns
+        -------
+        None: DataFrame updated inplace
+        """
 
-        df.drop("id", axis=1, inplace=True)
-        df_metadata.append("id")
-        df.attrs["id"] = id_set[0]
+        # Discard molecule column if unique
+        if "id" in df.columns:
+            id_set = df.id.unique()
+            if len(id_set) != 1:  # only 1 molecule supported ftm
+                raise NotImplementedError(
+                    "Only 1 molecule at a time is currently supported "
+                    + "in SpectrumFactory. Use radis.calc_spectrum, which "
+                    + "calculates them independently then use MergeSlabs"
+                )
 
-        isotope_set = df.iso.unique()
+            df.drop("id", axis=1, inplace=True)
+            df_metadata.append("id")
+            df.attrs["id"] = id_set[0]
+        else:
+            assert "id" in df.attrs or "molecule" in df.attrs
 
-        if len(isotope_set) == 1:
-            df.drop("iso", axis=1, inplace=True)
-            df_metadata.append("iso")
-            df.attrs["iso"] = isotope_set[0]
+        if "iso" in df.columns:
+            isotope_set = df.iso.unique()
 
-        return df
+            if len(isotope_set) == 1:
+                df.drop("iso", axis=1, inplace=True)
+                df_metadata.append("iso")
+                df.attrs["iso"] = isotope_set[0]
+        else:
+            assert "iso" in df.attrs
 
     def _get_isotope_list(self, molecule=None, df=None):
         """Returns list of isotopes for given molecule Parse the Input
@@ -1979,13 +2075,19 @@ def get_best_match():
                 return None
 
     def _build_partition_function_interpolator(
-        self, parfunc, parfuncfmt, molecule, isotope
+        self, parfunc, parfuncfmt, molecule, isotope, predefined_partition_functions={}
     ):
         """Returns an universal partition function object ``parsum`` with the
         following methods defined::
             parsum.at(T)
 
         Partition functions are interpolated from tabulated values
+
+        Other Parameters
+        ----------------
+        predefined_partition_functions: dict
+            ::
+                {molecule: {isotope: PartitionFunctionTabulator object}}
         """
 
         if __debug__:
@@ -1996,16 +2098,22 @@ def _build_partition_function_interpolator(
 
         isotope = int(isotope)
 
-        # Use HAPI (HITRAN Python interface, integrated in RADIS)
-        # no tabulated partition functions defined. Only non-eq spectra can
-        # be calculated if energies are also given
         if parfuncfmt == "hapi" or parfuncfmt is None:
+            assert len(predefined_partition_functions) == 0
+            # Use TIPS-2017 through HAPI (HITRAN Python interface, integrated in RADIS)
+            # no tabulated partition functions defined. Only non-eq spectra can
+            # be calculated if energies are also given
             parsum = PartFuncHAPI(
                 M=molecule, I=isotope, path=parfunc, verbose=self.verbose
             )
         elif parfuncfmt == "cdsd":  # Use tabulated CDSD partition functions
+            assert len(predefined_partition_functions) == 0
             assert molecule == "CO2"
             parsum = PartFuncCO2_CDSDtab(isotope, parfunc)
+        elif parfuncfmt == "exomol":
+            # Just read dictionary of predefined partition function
+            assert len(predefined_partition_functions) > 0
+            parsum = predefined_partition_functions[molecule][isotope]
         else:
             raise ValueError(
                 "Unknown format for partition function: {0}".format(parfuncfmt)
diff --git a/radis/levels/partfunc.py b/radis/levels/partfunc.py
index 3381ec7e8..31ffc015a 100644
--- a/radis/levels/partfunc.py
+++ b/radis/levels/partfunc.py
@@ -110,39 +110,8 @@ class RovibPartitionFunction(object):
     :class:`~radis.levels.partfunc.RovibParFuncCalculator`
     """
 
-    def __init__(self, electronic_state):
-
-        #        # Check inputs
-        #        if type(molecule) is int:
-        #            molecule = get_molecule(molecule)
-        #
-        #        # Get molecule from radis.db
-        #        mol_list = import_from_module('radis.db.molecules', molecule)
-        #        try:
-        #            mol = mol_list[isotope]
-        #        except KeyError:
-        #            raise KeyError('Isotope {0} not defined for molecule {1}'.format(int(isotope), molecule))
-
-        ElecState = electronic_state
-
-        try:
-            ElecState.Erovib
-        except AttributeError:
-            raise AttributeError(
-                "{0} has no energy function defined in RADIS".format(
-                    ElecState.get_fullname()
-                )
-            )
-
-        # Store
-        self.ElecState = ElecState  # electronic state object
-        self.molecule = ElecState.name  # molecule name
-        self.isotope = ElecState.iso
-
-        # Line database
-        # pandas Dataframe that holds all the lines
-        self.df = pd.DataFrame({})
-        # updated on inherited classes initialization
+    def __init__(self):
+        pass
 
 
 # %% Subclasses :
@@ -151,6 +120,9 @@ def __init__(self, electronic_state):
 
 
 class RovibParFuncTabulator(RovibPartitionFunction):
+    def __init__(self):
+        super(RovibParFuncTabulator, self).__init__()
+
     def at(self, T, **kwargs):
         """Get partition function at temperature T under equilibrium
         conditions, from tabulated data.
@@ -220,7 +192,28 @@ class RovibParFuncCalculator(RovibPartitionFunction):
 
     def __init__(self, electronic_state):
 
-        super(RovibParFuncCalculator, self).__init__(electronic_state=electronic_state)
+        super(RovibParFuncCalculator, self).__init__()
+
+        ElecState = electronic_state
+
+        try:
+            ElecState.Erovib
+        except AttributeError:
+            raise AttributeError(
+                "{0} has no energy function defined in RADIS".format(
+                    ElecState.get_fullname()
+                )
+            )
+
+        # Store
+        self.ElecState = ElecState  # electronic state object
+        self.molecule = ElecState.name  # molecule name
+        self.isotope = ElecState.iso
+
+        # Line database
+        # pandas Dataframe that holds all the lines
+        self.df = pd.DataFrame({})
+        # updated on inherited classes initialization
 
     def at(self, T, update_populations=False):
         """Get partition function at temperature T under equilibrium
@@ -690,6 +683,31 @@ def _get_rovib_populations(self):
 # %% Variants of Tabulated partition functions (interpolate)
 
 
+class PartFuncExoMol(RovibParFuncTabulator):
+    """Return partition function using interpolation of tabulated values.
+
+    Parameters
+    ----------
+    name: str
+        exomol isotope full name
+
+    See Also
+    --------
+    :py:func:`~radis.io.exomol.fetch_exomol`
+    """
+
+    def __init__(self, name, T_range, Q_range):
+
+        super(PartFuncExoMol, self).__init__()
+
+        self.name = name
+        self.T_range = T_range
+        self.Q_range = Q_range
+
+    def _at(self, T):
+        return np.interp(T, self.T_range, self.Q_range)
+
+
 class PartFuncHAPI(RovibParFuncTabulator):
     """Return partition function using interpolation of tabulated values.
 
@@ -727,6 +745,8 @@ class PartFuncHAPI(RovibParFuncTabulator):
 
     def __init__(self, M, I, path=None, verbose=True):
 
+        super(PartFuncHAPI, self).__init__()
+
         self.verbose = verbose
 
         if path is not None:
diff --git a/radis/los/slabs.py b/radis/los/slabs.py
index b7b01e02d..1ebacb259 100644
--- a/radis/los/slabs.py
+++ b/radis/los/slabs.py
@@ -40,7 +40,7 @@ def intersect(a, b):
     return c
 
 
-def SerialSlabs(*slabs, **kwargs):
+def SerialSlabs(*slabs, **kwargs) -> Spectrum:
     # type: (*Spectrum, **dict) -> Spectrum
     r"""Adds several slabs along the line-of-sight.
     If adding two slabs only, you can also use::
@@ -478,7 +478,7 @@ def same_wavespace(wl):
     return slabs
 
 
-def MergeSlabs(*slabs, **kwargs):
+def MergeSlabs(*slabs, **kwargs) -> Spectrum:
     # type: (*Spectrum, **dict) -> Spectrum
     r"""Combines several slabs into one. Useful to calculate multi-gas slabs.
     Linear absorption coefficient is calculated as the sum of all linear absorption
diff --git a/radis/misc/basics.py b/radis/misc/basics.py
index 9ff9570da..6f6595377 100644
--- a/radis/misc/basics.py
+++ b/radis/misc/basics.py
@@ -346,27 +346,32 @@ def transfer_metadata(df1, df2, metadata):
     """Transfer metadata between a DataFrame df1 and df2.
 
     For some reason metadata are sometimes not copied when a DataFrame is
-    sliced or copied, even if they explicitely figure in the df._metadata
-    attribute. Here we copy them back
+    sliced or copied, even if they explicitely appear in the ``df.attrs``
+    attribute. See https://github.com/pandas-dev/pandas/issues/28283
+
+    Here we copy them back. Attributes can be :
+        - keys of the ``df1.attrs`` dictionary
+        - simple attributes of ``df1``, i.e., ``df1.X``
 
     Parameters
     ----------
-
     df1: pandas DataFrame
         copy from df1
-
     df2: pandas DataFrame
         copy to df2
     """
 
     for k in metadata:
-        if __debug__ and k not in df1._metadata:
-            from radis.misc.debug import printdbg
-
-            printdbg("WARNING. {0} not in _metadata: {1}".format(k, df1._metadata))
         if not hasattr(df2, k):
             assert k not in df1.columns  # point is not to copy columns!
-            setattr(df2, k, getattr(df1, k))
+            if k in df1.attrs:  # Keys of attribute dictionary
+                df2.attrs[k] = df1.attrs[k]
+            else:  # Direct attributes of df1
+                setattr(df2, k, getattr(df1, k))
+
+    # @dev: refactor : we're updating the database to properly store values
+    # either in columns either in the .attrs dict, but so they can always
+    # be accessed with df.X
 
 
 def expand_metadata(df, metadata):
diff --git a/radis/misc/warning.py b/radis/misc/warning.py
index c6c57b65c..0ee6f807b 100644
--- a/radis/misc/warning.py
+++ b/radis/misc/warning.py
@@ -119,15 +119,33 @@ class NegativeEnergiesWarning(UserWarning):
     pass
 
 
+class MissingSelfBroadeningTdepWarning(UserWarning):
+    """Self broadening temperature dependance coefficient is missing in Line Database.
+
+    Usually, use Air broadening temperature dependance coefficient instead. See
+    :py:meth:`~radis.lbl.broadening.BroadenFactory._add_collisional_broadening_HWHM`
+    """
+
+    pass
+
+
 class MissingSelfBroadeningWarning(UserWarning):
-    """Self broadening is missing in Line Database.
+    """Self broadening tabulated width is missing in Line Database.
 
-    Usually, use Air broadening instead
+    Usually, use Air broadening tabulated width instead. See
+    :py:meth:`~radis.lbl.broadening.BroadenFactory._add_collisional_broadening_HWHM`
     """
 
     pass
 
 
+class MissingPressureShiftWarning(UserWarning):
+    """Pressure-shift coefficient is missing in Line Database."""
+
+    # TODO : add docstring link to references of line database columns.
+    pass
+
+
 class LinestrengthCutoffWarning(AccuracyWarning):
     """Warning triggered when the cumulated linestrength after intensity cutoff
     has changed too much."""
@@ -183,7 +201,9 @@ class DatabaseNotFoundError(FileNotFoundError):
     "OutOfRangeLinesWarning": OutOfRangeLinesWarning,
     "HighTemperatureWarning": HighTemperatureWarning,
     "NegativeEnergiesWarning": NegativeEnergiesWarning,
+    "MissingSelfBroadeningTdepWarning": MissingSelfBroadeningTdepWarning,
     "MissingSelfBroadeningWarning": MissingSelfBroadeningWarning,
+    "MissingPressureShiftWarning": MissingPressureShiftWarning,
     "LinestrengthCutoffWarning": LinestrengthCutoffWarning,
     "InputConditionsWarning": InputConditionsWarning,
     "DeprecatedFileWarning": DeprecatedFileWarning,
@@ -227,7 +247,9 @@ class DatabaseNotFoundError(FileNotFoundError):
     "HighTemperatureWarning": "warn",
     "NegativeEnergiesWarning": "warn",  # warning if negative energies in database
     # warning if self-broadening abs coefficnet missing (Air is used instead)
+    "MissingSelfBroadeningTdepWarning": "warn",
     "MissingSelfBroadeningWarning": "warn",
+    "MissingPressureShiftWarning": "warn",
     "InputConditionsWarning": "warn",
     "DeprecatedFileWarning": "warn",
     "IrrelevantFileWarning": "warn",
diff --git a/radis/spectrum/models.py b/radis/spectrum/models.py
index b6f4de257..bc0dd6889 100644
--- a/radis/spectrum/models.py
+++ b/radis/spectrum/models.py
@@ -46,7 +46,7 @@ def calculated_spectrum(
     cond_units=None,
     populations=None,
     name=None,
-):  # -> Spectrum:
+) -> Spectrum:
     """Convert ``(w, I)`` into a :py:class:`~radis.spectrum.spectrum.Spectrum`
     object that has unit conversion, plotting and slit convolution
     capabilities.
@@ -121,7 +121,7 @@ def calculated_spectrum(
 
 def transmittance_spectrum(
     w, T, wunit="nm", Tunit="", conditions=None, cond_units=None, name=None
-):  # -> Spectrum:
+) -> Spectrum:
     """Convert ``(w, I)`` into a :py:class:`~radis.spectrum.spectrum.Spectrum`
     object that has unit conversion, plotting and slit convolution
     capabilities.
@@ -191,7 +191,7 @@ def transmittance_spectrum(
 
 def experimental_spectrum(
     w, I, wunit="nm", Iunit="counts", conditions={}, cond_units=None, name=None
-):  # -> Spectrum:
+) -> Spectrum:
     """Convert ``(w, I)`` into a :py:class:`~radis.spectrum.spectrum.Spectrum`
     object that has unit conversion and plotting capabilities. Convolution is
     not available as the spectrum is assumed to have be measured experimentally
diff --git a/radis/spectrum/spectrum.py b/radis/spectrum/spectrum.py
index 8e2665527..58ebca889 100644
--- a/radis/spectrum/spectrum.py
+++ b/radis/spectrum/spectrum.py
@@ -3583,6 +3583,7 @@ def print_perf_profile(self, number_format="{:.3f}", precision=16):
         Example
         -------
         ::
+
             Spectrum.print_perf_profile()
 
             # output >>
@@ -3649,6 +3650,7 @@ def generate_perf_profile(self):
         Examples
         --------
         ::
+
             s = calc_spectrum(...)
             s.generate_perf_profile()
 
diff --git a/radis/spectrum/utils.py b/radis/spectrum/utils.py
index 1cb600552..fe47c24f5 100644
--- a/radis/spectrum/utils.py
+++ b/radis/spectrum/utils.py
@@ -461,6 +461,7 @@ def print_perf_profile(
     Example
     -------
     ::
+
         Spectrum.print_perf_profile()
 
         # output >>
diff --git a/setup.py b/setup.py
index 3b6986014..f9ce431ba 100644
--- a/setup.py
+++ b/setup.py
@@ -243,6 +243,9 @@ def run_setup(with_binary):
         packages=find_packages(),
         install_requires=[
             "hitran-api",
+            "beautifulsoup4",  # parse ExoMol website
+            "lxml",  # parser used for ExoMol website
+            "pyarrow",  # for the feather format (temporarily needed for ExoMol)
             "numpy",
             "scipy>=1.4.0",
             "matplotlib",  # ">=3.4.0" to suppress the Ruler warning, but only available for Python >= 3.7


Test patch:
diff --git a/radis/test/db/test_molparams.py b/radis/test/db/test_molparams.py
new file mode 100644
index 000000000..e4a17e880
--- /dev/null
+++ b/radis/test/db/test_molparams.py
@@ -0,0 +1,19 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul 21 00:47:30 2021
+
+@author: erwan
+"""
+
+from radis.db.molparam import MolParams
+
+
+def test_molparams(verbose=True, *args, **kwargs):
+
+    molpar = MolParams()
+
+    assert molpar.get(2, 2, "abundance") == 0.0110574
+
+
+if __name__ == "__main__":
+    test_molparams()
diff --git a/radis/test/io/test_exomol.py b/radis/test/io/test_exomol.py
new file mode 100644
index 000000000..c59579218
--- /dev/null
+++ b/radis/test/io/test_exomol.py
@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Jul 19 22:44:39 2021
+
+@author: erwan
+"""
+
+import pytest
+
+from radis.io.exomol import get_exomol_database_list, get_exomol_full_isotope_name
+
+
+def test_exomol_parsing_functions(verbose=True, *args, **kwargs):
+    """Test functions used to parse ExoMol website"""
+
+    assert get_exomol_full_isotope_name("H2O", 1) == "1H2-16O"
+    assert get_exomol_full_isotope_name("H2O", 4) == "1H-2H-16O"
+    assert get_exomol_full_isotope_name("CO2", 3) == "16O-12C-18O"
+    assert get_exomol_full_isotope_name("H2O", 1) == "1H2-16O"
+    assert get_exomol_full_isotope_name("NH3", 1) == "14N-1H3"
+    assert get_exomol_full_isotope_name("H2S", 1) == "1H2-32S"
+    assert get_exomol_full_isotope_name("FeH", 1) == "56Fe-1H"
+    assert get_exomol_full_isotope_name("CH4", 1) == "12C-1H4"
+    assert get_exomol_full_isotope_name("C2H4", 1) == "12C2-1H4"
+    assert get_exomol_full_isotope_name("C2H2", 1) == "12C2-1H2"
+
+    #  Note : this may change if new databases are added ; test would have to
+    # be updated in that case.
+    databases, recommended = get_exomol_database_list("CH4", "12C-1H4")
+    assert databases == ["xsec-YT10to10", "YT10to10", "YT34to10"]
+    assert recommended == "YT34to10"
+
+    # Test that these databases are found
+    KNOWN_EXOMOL_DATABASE_NAMES = {
+        "H2O": {"1H2-16O": ["POKAZATEL"]},
+        "NH3": {"14N-1H3": ["CoYuTe"]},
+        "H2S": {"1H2-32S": ["AYT2"]},
+        "FeH": {"56Fe-1H": ["MoLLIST"]},
+    }
+    for molecule, val in KNOWN_EXOMOL_DATABASE_NAMES.items():
+        for isotope_name, lookup_databases in val.items():
+            known_databases, _ = get_exomol_database_list(molecule, isotope_name)
+            for dat in lookup_databases:
+                assert dat in known_databases
+
+
+@pytest.mark.needs_connection
+def test_calc_exomol_spectrum(verbose=True, plot=True, *args, **kwargs):
+    """Auto-fetch and calculate a SiO spectrum from the ExoMol database
+
+    https://github.com/radis/radis/pull/320#issuecomment-884508206
+    """
+    from radis import calc_spectrum
+
+    s = calc_spectrum(
+        1080,
+        1320,  # cm-1
+        molecule="SiO",
+        isotope="1",
+        pressure=1.01325,  # bar
+        Tgas=1000,  # K
+        mole_fraction=0.1,
+        path_length=1,  # cm
+        databank="exomol",  # or use ('exomol', "EBJT") for a specific database ("EBJT")
+        verbose=verbose,
+    )
+    s.apply_slit(1, "cm-1")  # simulate an experimental slit
+    if plot:
+        s.plot("radiance")
+
+
+if __name__ == "__main__":
+    test_exomol_parsing_functions()
+    test_calc_exomol_spectrum()
diff --git a/radis/test/lbl/test_base.py b/radis/test/lbl/test_base.py
index e8af220bc..f3c0abf9e 100644
--- a/radis/test/lbl/test_base.py
+++ b/radis/test/lbl/test_base.py
@@ -71,8 +71,8 @@ def test_populations(plot=True, verbose=True, warnings=True, *args, **kwargs):
         plt.close()  # no need to keep it open, we just tested the function
 
     # Test calculated quantities are there
-    assert hasattr(sf.df1, "Qref")
-    assert hasattr(sf.df1, "Qvib")
+    assert "Qref" in sf.df1.attrs
+    assert "Qvib" in sf.df1.attrs
     assert hasattr(sf.df1, "Qrotu")
     assert hasattr(sf.df1, "Qrotl")
 
diff --git a/radis/test/validation/test_validation_vs_specair_noneqCO.py b/radis/test/validation/test_validation_vs_specair_noneqCO.py
index 2f67723cb..2c7ac5b96 100644
--- a/radis/test/validation/test_validation_vs_specair_noneqCO.py
+++ b/radis/test/validation/test_validation_vs_specair_noneqCO.py
@@ -122,8 +122,8 @@ def test_validation_vs_specair(rtol=1e-2, verbose=True, plot=False, *args, **kwa
     )
 
     # Compare partition functions to hardcoded values
-    b2 = np.isclose(s_2000_300.lines.Q, 139, atol=1)  # Specair: 139
-    b2 *= np.isclose(s_300_2000.lines.Q, 727, atol=1)  # Specair: 727
+    b2 = np.isclose(s_2000_300.lines.attrs["Q"], 139, atol=1)  # Specair: 139
+    b2 *= np.isclose(s_300_2000.lines.attrs["Q"], 727, atol=1)  # Specair: 727
 
     if verbose:
         printm(