Repository: radis/radis (language: python)

Issue:
Exomol api update
- avoid duplicate creation of *trans.bz2.hdf5 (without Sij0), *trans.hdf5 (with Sij0) files in ExoMol
- fail with more explicit error if missing nu_lines  (see #434)
- prepare refactor of ExoMolDb
- improve ExoMol molarmass parsing 
- fix reported errors in ExoMol definition file for C2H2. Fixes #434

Also : 
- add vaex/jax options to output a DataFrame loaded from HITRAN/HITEMP   (prepare Exojax/RADIS common API)
Bug Report: fetch_exomol() does not work for H2O, CO2, CH4, NH3, OH, NO, HCl
### 🐛 Describe the bug

``radis.io.exomol.fetch_exomol()`` does not work for H2O, CO2, CH4, NH3, OH, NO, HCl (and maybe for other species also, but I didn't test all species).
NB: ``radis.io.exomol.fetch_exomol()`` works perfectly for CO (and maybe for other species also, but I didn't test all species).

I tried the following script with molecule being ``"H2O"``, ``"CO2"``, ``"CH4"``, ``"NH3"``, ``"OH"``, ``"NO"``, or ``"HCl"``.
```
import radis
from radis.io.exomol import fetch_exomol
molecule="H2O" # "H2O", "CO2", "CH4", "NH3", "OH", "NO", "HCl"
fetch_exomol(
    molecule=molecule,
    database=None,
    isotope="1",
    cache='regen',
    verbose=True,
    clean_cache_files=True,
    return_local_path=False,
    engine="default",
)
```
The script was run on an HPC cluster (node: 1, cpu: 1, mem: 175GB) so there should be enough memory to load large databases.

For ``H2O``, I got the following error:
```
Using ExoMol database POKAZATEL for 1H2-16O (recommended by the ExoMol team). All available databases are ['BT2', 'VT', 'WAT-UV296', 'xsec-POKAZATEL', 'CKYKKY', 'xsec-VUV-DTU', 'YYZYK-HF', 'POKAZATEL']. 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))` 

Background atmosphere:  H2
Reading transition file
Note: Caching line transition data to the HDF5 format with vaex. After the second time, it will become much faster.
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 fetch_exomol(
      2         molecule="H2O",
      3         database=None,
      4         # local_databases=None,
      5         # databank_name="EXOMOL-{molecule}",
      6         isotope="1",
      7         # load_wavenum_min=0.,
      8         # load_wavenum_max=100000.,
      9         # columns=None,
     10         cache='regen',
     11         verbose=True,
     12         clean_cache_files=True,
     13         return_local_path=False,
     14         engine="default",
     15     )

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:370, in fetch_exomol(molecule, database, local_databases, databank_name, isotope, load_wavenum_min, load_wavenum_max, columns, cache, verbose, clean_cache_files, return_local_path, return_partition_function, engine)
    368 if load_wavenum_max is None:
    369     load_wavenum_max = np.inf
--> 370 mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max], engine=engine)
    372 # Attributes of the DataFrame
    373 from radis.db.classes import HITRAN_MOLECULES

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:769, in MdbExomol.__init__(self, path, nurange, margin, crit, bkgdatm, broadf, engine, verbose)
    767 b = trans.nu_lines.isnan()  # TODO: check if can be made faster?
    768 if b.sum() > 0:
--> 769     raise NotImplementedError(
    770         "transition wavenumber not given. Has to be computed from upper/lower level energy look up"
    771     )
    773 # TODO Load extra columns (precomputed Sij0, etc.)
    774 # trans_file_extra = self.trans_file.with_name(
    775 #    self.trans_file.stem + "_extra" + self.trans_file.suffix
   (...)
    782 
    783 # mask needs to be applied
    784 mask_needed = True

NotImplementedError: transition wavenumber not given. Has to be computed from upper/lower level energy look up

```

For ``CO2``:
```
Using ExoMol database UCL-4000 for 12C-16O2 (recommended by the ExoMol team). All available databases are ['xsec-UCL-4000', 'xsec-VUV-Venot-2018', 'AMES', 'Ames-2016', 'CDSD-4000', 'Venot', 'HITEMP', 'Zak-16-12-16', 'xsec-VUV-DTU', 'UCL-4000']. 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))` 

Background atmosphere:  H2
Downloading http://www.exomol.com/db/CO2/12C-16O2/12C-16O2__H2.broad
Error: Couldn't download .broad file and save.
Downloading http://www.exomol.com/db/CO2/12C-16O2/12C-16O2__He.broad
Error: Couldn't download .broad file and save.
Downloading http://www.exomol.com/db/CO2/12C-16O2/12C-16O2__air.broad
[05/17/22 09:47:04] ERROR    error opening '/gpfs/workdir/mariottop/.radisdb/exomol/CO __init__.py:259
                             2/12C-16O2/UCL-4000/12C-16O2__UCL-4000.states.bz2.hdf5'                  
                             Traceback (most recent call last):                                       
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/__init__.py", line 232, in open                
                                 ds = vaex.dataset.open(path, fs_options=fs_options,                  
                             fs=fs, **kwargs)                                                         
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/dataset.py", line 73, in open                  
                                 return opener.open(path, fs_options=fs_options,                      
                             fs=fs, *args, **kwargs)                                                  
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/dataset.py", line 1438, in open                
                                 return cls(path, *args, **kwargs)                                    
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 71, in                  
                             __init__                                                                 
                                 self._load()                                                         
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 202, in                 
                             _load                                                                    
                                 self._load_columns(self.h5file["/table"])                            
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 374, in                 
                             _load_columns                                                            
                                 elif hasattr(column["data"], "dtype"):                               
                               File "h5py/_objects.pyx", line 54, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/_objects.pyx", line 55, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/h5py/_hl/group.py", line 305, in                    
                             __getitem__                                                              
                                 oid = h5o.open(self.id, self._e(name),                               
                             lapl=self._lapl)                                                         
                               File "h5py/_objects.pyx", line 54, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/_objects.pyx", line 55, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/h5o.pyx", line 190, in h5py.h5o.open                        
                             KeyError: "Unable to open object (object 'data' doesn't                  
                             exist)"                                                                  
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 fetch_exomol(
      2         molecule="CO2",
      3         database=None,
      4         # local_databases=None,
      5         # databank_name="EXOMOL-{molecule}",
      6         isotope="1",
      7         # load_wavenum_min=0.,
      8         # load_wavenum_max=100000.,
      9         # columns=None,
     10         cache='regen',
     11         verbose=True,
     12         clean_cache_files=True,
     13         return_local_path=False,
     14         engine="default",
     15     )

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:370, in fetch_exomol(molecule, database, local_databases, databank_name, isotope, load_wavenum_min, load_wavenum_max, columns, cache, verbose, clean_cache_files, return_local_path, return_partition_function, engine)
    368 if load_wavenum_max is None:
    369     load_wavenum_max = np.inf
--> 370 mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max], engine=engine)
    372 # Attributes of the DataFrame
    373 from radis.db.classes import HITRAN_MOLECULES

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:577, in MdbExomol.__init__(self, path, nurange, margin, crit, bkgdatm, broadf, engine, verbose)
    575 if self.states_file.with_suffix(".bz2.hdf5").exists():
    576     states_file_cache = self.states_file.with_suffix(".bz2.hdf5")
--> 577     states = vaex.open(states_file_cache)
    578 elif self.states_file.with_suffix(".hdf5").exists():
    579     states_file_cache = self.states_file.with_suffix(".hdf5")

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/__init__.py:232, in open(path, convert, progress, shuffle, fs_options, fs, *args, **kwargs)
    230     ds = vaex.dataset.open(path_output, fs_options=fs_options, fs=fs, **kwargs)
    231 else:
--> 232     ds = vaex.dataset.open(path, fs_options=fs_options, fs=fs, **kwargs)
    233 df = vaex.from_dataset(ds)
    234 if df is None:

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/dataset.py:73, in open(path, fs_options, fs, *args, **kwargs)
     71     if opener.quick_test(path, fs_options=fs_options, fs=fs):
     72         if opener.can_open(path, fs_options=fs_options, fs=fs, *args, **kwargs):
---> 73             return opener.open(path, fs_options=fs_options, fs=fs, *args, **kwargs)
     75 # otherwise try all openers
     76 for opener in opener_classes:

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/dataset.py:1438, in DatasetFile.open(cls, path, *args, **kwargs)
   1436 @classmethod
   1437 def open(cls, path, *args, **kwargs):
-> 1438     return cls(path, *args, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:71, in Hdf5MemoryMapped.__init__(self, path, write, fs_options, fs, nommap, group, _fingerprint)
     69 self.group = group
     70 self._version = 1
---> 71 self._load()
     72 if not write:  # in write mode, call freeze yourself, so the hashes are computed
     73     self._freeze()

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:202, in Hdf5MemoryMapped._load(self)
    200 if "table" in self.h5file:
    201     self._version = 2
--> 202     self._load_columns(self.h5file["/table"])
    203     self.group = "/table"
    204 root_datasets = [dataset for name, dataset in self.h5file.items() if isinstance(dataset, h5py.Dataset)]

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:374, in Hdf5MemoryMapped._load_columns(self, h5data, first)
    372     column = self._map_hdf5_array(column)
    373     self.add_column(column_name, column)
--> 374 elif hasattr(column["data"], "dtype"):
    375     column = self._map_column(column)
    376     self.add_column(column_name, column)

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/h5py/_hl/group.py:305, in Group.__getitem__(self, name)
    303         raise ValueError("Invalid HDF5 object reference")
    304 elif isinstance(name, (bytes, str)):
--> 305     oid = h5o.open(self.id, self._e(name), lapl=self._lapl)
    306 else:
    307     raise TypeError("Accessing a group is done with bytes or str, "
    308                     " not {}".format(type(name)))

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File h5py/h5o.pyx:190, in h5py.h5o.open()

KeyError: "Unable to open object (object 'data' doesn't exist)"
```

For ``CH4``:
```
Using ExoMol database YT34to10 for 12C-1H4 (recommended by the ExoMol team). All available databases are ['xsec-YT10to10', 'YT10to10', 'YT34to10']. 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))` 

Background atmosphere:  H2
Reading transition file
Note: Caching line transition data to the HDF5 format with vaex. After the second time, it will become much faster.
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 fetch_exomol(
      2         molecule="CH4",
      3         database=None,
      4         # local_databases=None,
      5         # databank_name="EXOMOL-{molecule}",
      6         isotope="1",
      7         # load_wavenum_min=0.,
      8         # load_wavenum_max=100000.,
      9         # columns=None,
     10         cache='regen',
     11         verbose=True,
     12         clean_cache_files=True,
     13         return_local_path=False,
     14         engine="default",
     15     )

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:370, in fetch_exomol(molecule, database, local_databases, databank_name, isotope, load_wavenum_min, load_wavenum_max, columns, cache, verbose, clean_cache_files, return_local_path, return_partition_function, engine)
    368 if load_wavenum_max is None:
    369     load_wavenum_max = np.inf
--> 370 mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max], engine=engine)
    372 # Attributes of the DataFrame
    373 from radis.db.classes import HITRAN_MOLECULES

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:769, in MdbExomol.__init__(self, path, nurange, margin, crit, bkgdatm, broadf, engine, verbose)
    767 b = trans.nu_lines.isnan()  # TODO: check if can be made faster?
    768 if b.sum() > 0:
--> 769     raise NotImplementedError(
    770         "transition wavenumber not given. Has to be computed from upper/lower level energy look up"
    771     )
    773 # TODO Load extra columns (precomputed Sij0, etc.)
    774 # trans_file_extra = self.trans_file.with_name(
    775 #    self.trans_file.stem + "_extra" + self.trans_file.suffix
   (...)
    782 
    783 # mask needs to be applied
    784 mask_needed = True

NotImplementedError: transition wavenumber not given. Has to be computed from upper/lower level energy look up
```

For ``NH3``: 
```
Using ExoMol database CoYuTe for 14N-1H3 (recommended by the ExoMol team). All available databases are ['BYTe', 'C2018', 'xsec-VUV-DTU', 'xsec-BYTe', 'CoYuTe']. 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))` 

Background atmosphere:  H2
Note: Caching states data to the hdf5 format with vaex. After the second time, it will become much faster.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File /gpfs/workdir/mariottop/radis/radis/io/exomolapi.py:254, in read_states(statesf, dic_def, engine)
    253 try:
--> 254     dat = vaex.from_csv(
    255         statesf,
    256         compression="bz2",
    257         sep=r"\s+",
    258         usecols=usecol,
    259         names=names,
    260         convert=True,
    261     )
    262 except:

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/__init__.py:567, in from_csv(filename_or_buffer, copy_index, chunk_size, convert, fs_options, fs, progress, **kwargs)
    566 path_output = convert if isinstance(convert, str) else vaex.convert._convert_name(filename_or_buffer)
--> 567 vaex.convert.convert_csv(
    568     path_input=filename_or_buffer, fs_options_input=fs_options, fs_input=fs,
    569     path_output=path_output, fs_options_output=fs_options, fs_output=fs,
    570     chunk_size=chunk_size,
    571     copy_index=copy_index,
    572     progress=progress,
    573     **kwargs
    574 )
    575 return open(path_output, fs_options=fs_options, fs=fs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/convert.py:52, in convert_csv(path_input, fs_options_input, fs_input, path_output, fs_options_output, fs_output, progress, *args, **kwargs)
     51     _from_csv_convert_and_read(path_input, path_output=path_output, fs_options=fs_options_input, fs=fs_input, progress=progress, **kwargs)
---> 52 cached_output(*args, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/cache.py:422, in output_file.<locals>.wrapper1.<locals>.wrapper2.<locals>.call(callable, *args, **kwargs)
    421 log.info('file %s does not exist yet, running conversion %s → %s', path_output_meta, path_input, path_output)
--> 422 value = callable(*args, **kwargs)
    423 write_fingerprint()

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/convert.py:51, in convert_csv.<locals>.cached_output(*args, **kwargs)
     50     kwargs['chunk_size'] = 5_000_000
---> 51 _from_csv_convert_and_read(path_input, path_output=path_output, fs_options=fs_options_input, fs=fs_input, progress=progress, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/convert.py:70, in _from_csv_convert_and_read(filename_or_buffer, path_output, chunk_size, fs_options, fs, copy_index, progress, **kwargs)
     69 with vaex.file.open(filename_or_buffer, fs_options=fs_options, fs=fs, for_arrow=True) as f:
---> 70     csv_reader = pd.read_csv(filename_or_buffer, chunksize=chunk_size, **kwargs)
     71     for i, df_pandas in enumerate(csv_reader):

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/pandas/util/_decorators.py:311, in deprecate_nonkeyword_arguments.<locals>.decorate.<locals>.wrapper(*args, **kwargs)
    306     warnings.warn(
    307         msg.format(arguments=arguments),
    308         FutureWarning,
    309         stacklevel=stacklevel,
    310     )
--> 311 return func(*args, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/pandas/io/parsers/readers.py:680, in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, error_bad_lines, warn_bad_lines, on_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options)
    678 kwds.update(kwds_defaults)
--> 680 return _read(filepath_or_buffer, kwds)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/pandas/io/parsers/readers.py:572, in _read(filepath_or_buffer, kwds)
    571 # Check for duplicates in names.
--> 572 _validate_names(kwds.get("names", None))
    574 # Create the parser.

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/pandas/io/parsers/readers.py:532, in _validate_names(names)
    531 if len(names) != len(set(names)):
--> 532     raise ValueError("Duplicate names are not allowed.")
    533 if not (
    534     is_list_like(names, allow_sets=False) or isinstance(names, abc.KeysView)
    535 ):

ValueError: Duplicate names are not allowed.

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 fetch_exomol(
      2         molecule="NH3",
      3         database=None,
      4         # local_databases=None,
      5         # databank_name="EXOMOL-{molecule}",
      6         isotope="1",
      7         # load_wavenum_min=0.,
      8         # load_wavenum_max=100000.,
      9         # columns=None,
     10         cache='regen',
     11         verbose=True,
     12         clean_cache_files=True,
     13         return_local_path=False,
     14         engine="default",
     15     )

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:370, in fetch_exomol(molecule, database, local_databases, databank_name, isotope, load_wavenum_min, load_wavenum_max, columns, cache, verbose, clean_cache_files, return_local_path, return_partition_function, engine)
    368 if load_wavenum_max is None:
    369     load_wavenum_max = np.inf
--> 370 mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max], engine=engine)
    372 # Attributes of the DataFrame
    373 from radis.db.classes import HITRAN_MOLECULES

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:585, in MdbExomol.__init__(self, path, nurange, margin, crit, bkgdatm, broadf, engine, verbose)
    581     else:
    582         print(
    583             "Note: Caching states data to the hdf5 format with vaex. After the second time, it will become much faster."
    584         )
--> 585         states = exomolapi.read_states(self.states_file, dic_def, engine="vaex")
    586     ndstates = vaex.array_types.to_numpy(states)
    588 # load pf

File /gpfs/workdir/mariottop/radis/radis/io/exomolapi.py:263, in read_states(statesf, dic_def, engine)
    254         dat = vaex.from_csv(
    255             statesf,
    256             compression="bz2",
   (...)
    260             convert=True,
    261         )
    262     except:
--> 263         dat = vaex.read_csv(
    264             statesf, sep=r"\s+", usecols=usecol, names=names, convert=True
    265         )
    266 elif engine == "csv":
    267     try:

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/__init__.py:594, in read_csv(filepath_or_buffer, **kwargs)
    592 def read_csv(filepath_or_buffer, **kwargs):
    593     '''Alias to from_csv.'''
--> 594     return from_csv(filepath_or_buffer, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/__init__.py:567, in from_csv(filename_or_buffer, copy_index, chunk_size, convert, fs_options, fs, progress, **kwargs)
    565 import vaex.convert
    566 path_output = convert if isinstance(convert, str) else vaex.convert._convert_name(filename_or_buffer)
--> 567 vaex.convert.convert_csv(
    568     path_input=filename_or_buffer, fs_options_input=fs_options, fs_input=fs,
    569     path_output=path_output, fs_options_output=fs_options, fs_output=fs,
    570     chunk_size=chunk_size,
    571     copy_index=copy_index,
    572     progress=progress,
    573     **kwargs
    574 )
    575 return open(path_output, fs_options=fs_options, fs=fs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/convert.py:52, in convert_csv(path_input, fs_options_input, fs_input, path_output, fs_options_output, fs_output, progress, *args, **kwargs)
     50         kwargs['chunk_size'] = 5_000_000
     51     _from_csv_convert_and_read(path_input, path_output=path_output, fs_options=fs_options_input, fs=fs_input, progress=progress, **kwargs)
---> 52 cached_output(*args, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/cache.py:422, in output_file.<locals>.wrapper1.<locals>.wrapper2.<locals>.call(callable, *args, **kwargs)
    420 if not vaex.file.exists(path_output, fs_options=fs_options_output_, fs=fs_output):
    421     log.info('file %s does not exist yet, running conversion %s → %s', path_output_meta, path_input, path_output)
--> 422     value = callable(*args, **kwargs)
    423     write_fingerprint()
    424     return value

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/convert.py:51, in convert_csv.<locals>.cached_output(*args, **kwargs)
     48 if 'chunk_size' not in kwargs:
     49     # make it memory efficient by default
     50     kwargs['chunk_size'] = 5_000_000
---> 51 _from_csv_convert_and_read(path_input, path_output=path_output, fs_options=fs_options_input, fs=fs_input, progress=progress, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/convert.py:70, in _from_csv_convert_and_read(filename_or_buffer, path_output, chunk_size, fs_options, fs, copy_index, progress, **kwargs)
     68     print("Converting csv to chunk files")
     69 with vaex.file.open(filename_or_buffer, fs_options=fs_options, fs=fs, for_arrow=True) as f:
---> 70     csv_reader = pd.read_csv(filename_or_buffer, chunksize=chunk_size, **kwargs)
     71     for i, df_pandas in enumerate(csv_reader):
     72         df = vaex.from_pandas(df_pandas, copy_index=copy_index)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/pandas/util/_decorators.py:311, in deprecate_nonkeyword_arguments.<locals>.decorate.<locals>.wrapper(*args, **kwargs)
    305 if len(args) > num_allow_args:
    306     warnings.warn(
    307         msg.format(arguments=arguments),
    308         FutureWarning,
    309         stacklevel=stacklevel,
    310     )
--> 311 return func(*args, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/pandas/io/parsers/readers.py:680, in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, error_bad_lines, warn_bad_lines, on_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options)
    665 kwds_defaults = _refine_defaults_read(
    666     dialect,
    667     delimiter,
   (...)
    676     defaults={"delimiter": ","},
    677 )
    678 kwds.update(kwds_defaults)
--> 680 return _read(filepath_or_buffer, kwds)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/pandas/io/parsers/readers.py:572, in _read(filepath_or_buffer, kwds)
    569 nrows = kwds.get("nrows", None)
    571 # Check for duplicates in names.
--> 572 _validate_names(kwds.get("names", None))
    574 # Create the parser.
    575 parser = TextFileReader(filepath_or_buffer, **kwds)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/pandas/io/parsers/readers.py:532, in _validate_names(names)
    530 if names is not None:
    531     if len(names) != len(set(names)):
--> 532         raise ValueError("Duplicate names are not allowed.")
    533     if not (
    534         is_list_like(names, allow_sets=False) or isinstance(names, abc.KeysView)
    535     ):
    536         raise ValueError("Names should be an ordered collection.")

ValueError: Duplicate names are not allowed.
```
For ``OH``:
```
Using ExoMol database MoLLIST for 16O-1H (recommended by the ExoMol team). All available databases are ['HITEMP', 'xsec-HITEMP', 'MoLLIST']. 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))` 

Background atmosphere:  H2
Downloading http://www.exomol.com/db/OH/16O-1H/16O-1H__H2.broad
Error: Couldn't download .broad file and save.
Downloading http://www.exomol.com/db/OH/16O-1H/16O-1H__He.broad
Error: Couldn't download .broad file and save.
Downloading http://www.exomol.com/db/OH/16O-1H/16O-1H__air.broad
Error: Couldn't download .broad file and save.
/gpfs/workdir/mariottop/radis/radis/io/exomol.py:542: Warning: 
                    No default broadening exponent in def file. Assigned n = 0.5
                    
  warnings.warn(
/gpfs/workdir/mariottop/radis/radis/io/exomol.py:552: Warning: 
                    No default broadening in def file. Assigned alpha_ref = 0.07
                    
  warnings.warn(
[05/17/22 09:51:32] ERROR    error opening '/gpfs/workdir/mariottop/.radisdb/exomol/OH __init__.py:259
                             /16O-1H/MoLLIST/16O-1H__MoLLIST.states.bz2.hdf5'                         
                             Traceback (most recent call last):                                       
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/__init__.py", line 232, in open                
                                 ds = vaex.dataset.open(path, fs_options=fs_options,                  
                             fs=fs, **kwargs)                                                         
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/dataset.py", line 73, in open                  
                                 return opener.open(path, fs_options=fs_options,                      
                             fs=fs, *args, **kwargs)                                                  
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/dataset.py", line 1438, in open                
                                 return cls(path, *args, **kwargs)                                    
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 71, in                  
                             __init__                                                                 
                                 self._load()                                                         
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 202, in                 
                             _load                                                                    
                                 self._load_columns(self.h5file["/table"])                            
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 374, in                 
                             _load_columns                                                            
                                 elif hasattr(column["data"], "dtype"):                               
                               File "h5py/_objects.pyx", line 54, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/_objects.pyx", line 55, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/h5py/_hl/group.py", line 305, in                    
                             __getitem__                                                              
                                 oid = h5o.open(self.id, self._e(name),                               
                             lapl=self._lapl)                                                         
                               File "h5py/_objects.pyx", line 54, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/_objects.pyx", line 55, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/h5o.pyx", line 190, in h5py.h5o.open                        
                             KeyError: "Unable to open object (object 'data' doesn't                  
                             exist)"                                                                  
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Input In [10], in <cell line: 1>()
----> 1 fetch_exomol(
      2         molecule="OH",
      3         database=None,
      4         # local_databases=None,
      5         # databank_name="EXOMOL-{molecule}",
      6         isotope="1",
      7         # load_wavenum_min=0.,
      8         # load_wavenum_max=100000.,
      9         # columns=None,
     10         cache='regen',
     11         verbose=True,
     12         clean_cache_files=True,
     13         return_local_path=False,
     14         engine="default",
     15     )

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:370, in fetch_exomol(molecule, database, local_databases, databank_name, isotope, load_wavenum_min, load_wavenum_max, columns, cache, verbose, clean_cache_files, return_local_path, return_partition_function, engine)
    368 if load_wavenum_max is None:
    369     load_wavenum_max = np.inf
--> 370 mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max], engine=engine)
    372 # Attributes of the DataFrame
    373 from radis.db.classes import HITRAN_MOLECULES

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:577, in MdbExomol.__init__(self, path, nurange, margin, crit, bkgdatm, broadf, engine, verbose)
    575 if self.states_file.with_suffix(".bz2.hdf5").exists():
    576     states_file_cache = self.states_file.with_suffix(".bz2.hdf5")
--> 577     states = vaex.open(states_file_cache)
    578 elif self.states_file.with_suffix(".hdf5").exists():
    579     states_file_cache = self.states_file.with_suffix(".hdf5")

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/__init__.py:232, in open(path, convert, progress, shuffle, fs_options, fs, *args, **kwargs)
    230     ds = vaex.dataset.open(path_output, fs_options=fs_options, fs=fs, **kwargs)
    231 else:
--> 232     ds = vaex.dataset.open(path, fs_options=fs_options, fs=fs, **kwargs)
    233 df = vaex.from_dataset(ds)
    234 if df is None:

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/dataset.py:73, in open(path, fs_options, fs, *args, **kwargs)
     71     if opener.quick_test(path, fs_options=fs_options, fs=fs):
     72         if opener.can_open(path, fs_options=fs_options, fs=fs, *args, **kwargs):
---> 73             return opener.open(path, fs_options=fs_options, fs=fs, *args, **kwargs)
     75 # otherwise try all openers
     76 for opener in opener_classes:

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/dataset.py:1438, in DatasetFile.open(cls, path, *args, **kwargs)
   1436 @classmethod
   1437 def open(cls, path, *args, **kwargs):
-> 1438     return cls(path, *args, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:71, in Hdf5MemoryMapped.__init__(self, path, write, fs_options, fs, nommap, group, _fingerprint)
     69 self.group = group
     70 self._version = 1
---> 71 self._load()
     72 if not write:  # in write mode, call freeze yourself, so the hashes are computed
     73     self._freeze()

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:202, in Hdf5MemoryMapped._load(self)
    200 if "table" in self.h5file:
    201     self._version = 2
--> 202     self._load_columns(self.h5file["/table"])
    203     self.group = "/table"
    204 root_datasets = [dataset for name, dataset in self.h5file.items() if isinstance(dataset, h5py.Dataset)]

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:374, in Hdf5MemoryMapped._load_columns(self, h5data, first)
    372     column = self._map_hdf5_array(column)
    373     self.add_column(column_name, column)
--> 374 elif hasattr(column["data"], "dtype"):
    375     column = self._map_column(column)
    376     self.add_column(column_name, column)

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/h5py/_hl/group.py:305, in Group.__getitem__(self, name)
    303         raise ValueError("Invalid HDF5 object reference")
    304 elif isinstance(name, (bytes, str)):
--> 305     oid = h5o.open(self.id, self._e(name), lapl=self._lapl)
    306 else:
    307     raise TypeError("Accessing a group is done with bytes or str, "
    308                     " not {}".format(type(name)))

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File h5py/h5o.pyx:190, in h5py.h5o.open()

KeyError: "Unable to open object (object 'data' doesn't exist)"
```

For ``NO``:
```
Using ExoMol database XABC for 14N-16O (recommended by the ExoMol team). All available databases are ['HITEMP', 'xsec-XABC', 'NOname', 'XABC']. 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))` 

Background atmosphere:  H2
Downloading http://www.exomol.com/db/NO/14N-16O/14N-16O__H2.broad
Error: Couldn't download .broad file and save.
Downloading http://www.exomol.com/db/NO/14N-16O/14N-16O__He.broad
Error: Couldn't download .broad file and save.
Downloading http://www.exomol.com/db/NO/14N-16O/14N-16O__air.broad
[05/17/22 09:52:20] ERROR    error opening '/gpfs/workdir/mariottop/.radisdb/exomol/NO __init__.py:259
                             /14N-16O/XABC/14N-16O__XABC.states.bz2.hdf5'                             
                             Traceback (most recent call last):                                       
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/__init__.py", line 232, in open                
                                 ds = vaex.dataset.open(path, fs_options=fs_options,                  
                             fs=fs, **kwargs)                                                         
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/dataset.py", line 73, in open                  
                                 return opener.open(path, fs_options=fs_options,                      
                             fs=fs, *args, **kwargs)                                                  
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/dataset.py", line 1438, in open                
                                 return cls(path, *args, **kwargs)                                    
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 71, in                  
                             __init__                                                                 
                                 self._load()                                                         
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 202, in                 
                             _load                                                                    
                                 self._load_columns(self.h5file["/table"])                            
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/vaex/hdf5/dataset.py", line 374, in                 
                             _load_columns                                                            
                                 elif hasattr(column["data"], "dtype"):                               
                               File "h5py/_objects.pyx", line 54, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/_objects.pyx", line 55, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "/gpfs/users/mariottop/.conda/envs/pynabam/lib/pyt                
                             hon3.8/site-packages/h5py/_hl/group.py", line 305, in                    
                             __getitem__                                                              
                                 oid = h5o.open(self.id, self._e(name),                               
                             lapl=self._lapl)                                                         
                               File "h5py/_objects.pyx", line 54, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/_objects.pyx", line 55, in                                  
                             h5py._objects.with_phil.wrapper                                          
                               File "h5py/h5o.pyx", line 190, in h5py.h5o.open                        
                             KeyError: "Unable to open object (object 'data' doesn't                  
                             exist)"                                                                  
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 fetch_exomol(
      2         molecule="NO",
      3         database=None,
      4         # local_databases=None,
      5         # databank_name="EXOMOL-{molecule}",
      6         isotope="1",
      7         # load_wavenum_min=0.,
      8         # load_wavenum_max=100000.,
      9         # columns=None,
     10         cache='regen',
     11         verbose=True,
     12         clean_cache_files=True,
     13         return_local_path=False,
     14         engine="default",
     15     )

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:370, in fetch_exomol(molecule, database, local_databases, databank_name, isotope, load_wavenum_min, load_wavenum_max, columns, cache, verbose, clean_cache_files, return_local_path, return_partition_function, engine)
    368 if load_wavenum_max is None:
    369     load_wavenum_max = np.inf
--> 370 mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max], engine=engine)
    372 # Attributes of the DataFrame
    373 from radis.db.classes import HITRAN_MOLECULES

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:577, in MdbExomol.__init__(self, path, nurange, margin, crit, bkgdatm, broadf, engine, verbose)
    575 if self.states_file.with_suffix(".bz2.hdf5").exists():
    576     states_file_cache = self.states_file.with_suffix(".bz2.hdf5")
--> 577     states = vaex.open(states_file_cache)
    578 elif self.states_file.with_suffix(".hdf5").exists():
    579     states_file_cache = self.states_file.with_suffix(".hdf5")

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/__init__.py:232, in open(path, convert, progress, shuffle, fs_options, fs, *args, **kwargs)
    230     ds = vaex.dataset.open(path_output, fs_options=fs_options, fs=fs, **kwargs)
    231 else:
--> 232     ds = vaex.dataset.open(path, fs_options=fs_options, fs=fs, **kwargs)
    233 df = vaex.from_dataset(ds)
    234 if df is None:

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/dataset.py:73, in open(path, fs_options, fs, *args, **kwargs)
     71     if opener.quick_test(path, fs_options=fs_options, fs=fs):
     72         if opener.can_open(path, fs_options=fs_options, fs=fs, *args, **kwargs):
---> 73             return opener.open(path, fs_options=fs_options, fs=fs, *args, **kwargs)
     75 # otherwise try all openers
     76 for opener in opener_classes:

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/dataset.py:1438, in DatasetFile.open(cls, path, *args, **kwargs)
   1436 @classmethod
   1437 def open(cls, path, *args, **kwargs):
-> 1438     return cls(path, *args, **kwargs)

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:71, in Hdf5MemoryMapped.__init__(self, path, write, fs_options, fs, nommap, group, _fingerprint)
     69 self.group = group
     70 self._version = 1
---> 71 self._load()
     72 if not write:  # in write mode, call freeze yourself, so the hashes are computed
     73     self._freeze()

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:202, in Hdf5MemoryMapped._load(self)
    200 if "table" in self.h5file:
    201     self._version = 2
--> 202     self._load_columns(self.h5file["/table"])
    203     self.group = "/table"
    204 root_datasets = [dataset for name, dataset in self.h5file.items() if isinstance(dataset, h5py.Dataset)]

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/vaex/hdf5/dataset.py:374, in Hdf5MemoryMapped._load_columns(self, h5data, first)
    372     column = self._map_hdf5_array(column)
    373     self.add_column(column_name, column)
--> 374 elif hasattr(column["data"], "dtype"):
    375     column = self._map_column(column)
    376     self.add_column(column_name, column)

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File ~/.conda/envs/pynabam/lib/python3.8/site-packages/h5py/_hl/group.py:305, in Group.__getitem__(self, name)
    303         raise ValueError("Invalid HDF5 object reference")
    304 elif isinstance(name, (bytes, str)):
--> 305     oid = h5o.open(self.id, self._e(name), lapl=self._lapl)
    306 else:
    307     raise TypeError("Accessing a group is done with bytes or str, "
    308                     " not {}".format(type(name)))

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File h5py/h5o.pyx:190, in h5py.h5o.open()

KeyError: "Unable to open object (object 'data' doesn't exist)"
```

For ``HCl``:
```
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 fetch_exomol(
      2         molecule="HCl",
      3         database=None,
      4         # local_databases=None,
      5         # databank_name="EXOMOL-{molecule}",
      6         isotope="1",
      7         # load_wavenum_min=0.,
      8         # load_wavenum_max=100000.,
      9         # columns=None,
     10         cache='regen',
     11         verbose=True,
     12         clean_cache_files=True,
     13         return_local_path=False,
     14         engine="default",
     15     )

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:328, in fetch_exomol(molecule, database, local_databases, databank_name, isotope, load_wavenum_min, load_wavenum_max, columns, cache, verbose, clean_cache_files, return_local_path, return_partition_function, engine)
    322     raise ValueError(
    323         f"In fetch_exomol, ``isotope`` must be an integer. Got `{isotope}` "
    324         + "Only one isotope can be queried at a time. "
    325     )
    327 full_molecule_name = get_exomol_full_isotope_name(molecule, isotope)
--> 328 known_exomol_databases, recommended_database = get_exomol_database_list(
    329     molecule, full_molecule_name
    330 )
    332 _exomol_use_hint = "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"
    333 if database is None or database == "default":

File /gpfs/workdir/mariottop/radis/radis/io/exomol.py:137, in get_exomol_database_list(molecule, isotope_full_name)
    131         print(
    132             f"Multiple recommended databases found for {molecule} in ExoMol : {databases_recommended}. This is unexpected. Using the first"
    133         )
    135 databases = databases + databases_recommended
--> 137 return databases, databases_recommended[0]

IndexError: list index out of range
```



### 💡 Possible solutions

_No response_

### 🎲 Radis version

0.13.0

### 💻 Operating system

Linux


Discussion / hints:
# [Codecov](https://codecov.io/gh/radis/radis/pull/464?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) Report
> Merging [#464](https://codecov.io/gh/radis/radis/pull/464?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) (344214d) into [develop](https://codecov.io/gh/radis/radis/commit/a0c9bf3d9444b05ad4196959398eddf1879f8d40?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) (a0c9bf3) will **decrease** coverage by `0.13%`.
> The diff coverage is `51.85%`.


```diff
@@             Coverage Diff             @@
##           develop     #464      +/-   ##
===========================================
- Coverage    77.13%   76.99%   -0.14%     
===========================================
  Files          161      161              
  Lines        18687    18713      +26     
===========================================
- Hits         14414    14409       -5     
- Misses        4273     4304      +31     
```



Fix patch:
diff --git a/.travis.yml b/.travis.yml
index 955bf3183..bc9b01c22 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -46,9 +46,12 @@ jobs:
         # Your test script goes here
         - echo ">>> Run tests"
         # use XVFB to have headless display port, and still run the Matplotlib tests.
-        - xvfb-run pytest --cov=./
+        - xvfb-run -a pytest --cov=./
         # lookup 'addopts' in setup.cfg>[tools:pytest] for default tests
 
+      services:
+        - xvfb
+
       after_success:
         # codecoverage
         - pip install codecov
diff --git a/examples/plot_compare_CO_exomol_hitemp.py b/examples/plot_compare_CO_exomol_hitemp.py
index 40ee2ffee..b9c692452 100644
--- a/examples/plot_compare_CO_exomol_hitemp.py
+++ b/examples/plot_compare_CO_exomol_hitemp.py
@@ -27,7 +27,7 @@
     "wmin": 2002 / u.cm,
     "wmax": 2300 / u.cm,
     "molecule": "CO",
-    "isotope": "2",
+    "isotope": "1",
     "pressure": 1.01325,  # bar
     "Tgas": 1000,  # K
     "mole_fraction": 0.1,
@@ -37,7 +37,7 @@
 }
 
 s_exomol = calc_spectrum(
-    **conditions, databank="exomol", name="ExoMol's HITEMP (H2 broadened)"
+    **conditions, databank="exomol", name="ExoMol's HITEMP (default broadening)"
 )
 s_hitemp = calc_spectrum(
     **conditions,
diff --git a/radis/db/classes.py b/radis/db/classes.py
index e99fdba16..85c3ee003 100644
--- a/radis/db/classes.py
+++ b/radis/db/classes.py
@@ -248,7 +248,7 @@ def get_molecule_identifier(molecule_name):
     except KeyError:
         raise NotImplementedError(
             "Molecule '{0}' not supported. Choose one of {1}".format(
-                molecule_name, list(trans2.keys())
+                molecule_name, sorted(list(trans2.keys()))
             )
         )
 
diff --git a/radis/io/cache_files.py b/radis/io/cache_files.py
index f745de919..efac7845f 100644
--- a/radis/io/cache_files.py
+++ b/radis/io/cache_files.py
@@ -29,7 +29,7 @@
 # Note: don't import unicode_literals because it breaks the df.to_hdf of
 # save_to_hdf because of a stupid unicode/str error in Python 2.7
 import os
-from os.path import exists, splitext
+from os.path import exists
 from warnings import warn
 
 from packaging.version import parse
@@ -37,9 +37,9 @@
 import radis
 
 try:
-    from .hdf5 import HDF5Manager
+    from .hdf5 import DataFileManager
 except ImportError:  # if file is ran as a module
-    from radis.io.hdf5 import HDF5Manager
+    from radis.io.hdf5 import DataFileManager
 from radis.misc.basics import compare_dict, is_float
 from radis.misc.printer import printm, printr
 from radis.misc.warning import DeprecatedFileWarning, IrrelevantFileWarning
@@ -179,8 +179,8 @@ def load_h5_cache_file(
         printm("Reading cache file ({0})".format(cachefile))
     try:
         # Load file :
-        manager = HDF5Manager(engine)
-        df = manager.load(cachefile, columns=columns, key="df")
+        manager = DataFileManager(engine)
+        df = manager.read(cachefile, columns=columns, key="df")
 
     except KeyError as err:  # An error happened during file reading.
         # Fail safe by deleting cache file (unless we explicitely wanted it
@@ -222,8 +222,8 @@ def get_cache_file(fcache, engine="pytables", verbose=True):
     # TODO @dev: refactor : we probably don't need this function (only used in astroquery)
 
     # Load file
-    manager = HDF5Manager(engine)
-    df = manager.load(fcache, key="df")
+    manager = DataFileManager(engine)
+    df = manager.read(fcache, key="df")
 
     # Check file
     # ... 'object' columns slow everything down (not fixed format strings!)
@@ -368,10 +368,10 @@ def check_not_deprecated(
         which HDF5 library to use. If ``'guess'``, try to guess.
     """
     if engine == "guess":
-        engine = HDF5Manager.guess_engine(file)
+        engine = DataFileManager.guess_engine(file)
 
     # Get metadata :
-    manager = HDF5Manager(engine)
+    manager = DataFileManager(engine)
 
     try:
         file_metadata = manager.read_metadata(file)
@@ -526,10 +526,10 @@ def check_relevancy(
 
     """
     if engine == "guess":
-        engine = HDF5Manager.guess_engine(file)
+        engine = DataFileManager.guess_engine(file)
 
     # Get metadata :
-    manager = HDF5Manager(engine)
+    manager = DataFileManager(engine)
     file_metadata = manager.read_metadata(file, key=key)
 
     for k, v in relevant_if_metadata_above.items():
@@ -607,7 +607,7 @@ def save_to_hdf(
      ``None`` values are not stored
     """
     # Check file
-    assert fname.endswith(".h5") or fname.endswith(".hdf5")
+    assert str(fname).endswith(".h5") or str(fname).endswith(".hdf5")
     assert "version" not in metadata
     # ... 'object' columns slow everything down (not fixed format strings!)
     if verbose >= 2:
@@ -621,7 +621,7 @@ def save_to_hdf(
         raise ValueError("File exist: {0}".format(fname))
 
     # start by exporting dataframe
-    manager = HDF5Manager(engine)
+    manager = DataFileManager(engine)
     manager.write(fname, df, append=False, key=key)
 
     # Add metadata
@@ -686,19 +686,7 @@ def some_function(*args):
 
 
 def cache_file_name(fname, engine="pytables"):
-    """Return the corresponding cache file name for fname.
-
-    Other Parameters
-    ----------------
-    engine: ``'h5py'``, ``'pytables'``, ``'vaex'``
-       which HDF5 library to use. Default ``pytables``
-    """
-    if engine in ["pytables", "pytables-fixed"]:
-        return splitext(fname)[0] + ".h5"
-    elif engine in ["h5py", "vaex"]:
-        return splitext(fname)[0] + ".hdf5"
-    else:
-        raise ValueError(engine)
+    raise DeprecationWarning("Use DataFileManager.cache_file() instead")
 
 
 if __name__ == "__main__":
diff --git a/radis/io/cdsd.py b/radis/io/cdsd.py
index d3518b50a..13cac495c 100644
--- a/radis/io/cdsd.py
+++ b/radis/io/cdsd.py
@@ -23,7 +23,8 @@
 from os.path import exists, getmtime
 
 import radis
-from radis.io.cache_files import cache_file_name, load_h5_cache_file, save_to_hdf
+from radis.io.cache_files import load_h5_cache_file, save_to_hdf
+from radis.io.hdf5 import DataFileManager
 from radis.io.tools import (
     drop_object_format_columns,
     parse_hitran_file,
@@ -241,7 +242,7 @@ def cdsd2df(
         raise ValueError("Unknown CDSD version: {0}".format(version))
 
     # Use cache file if possible
-    fcache = cache_file_name(fname, engine=engine)
+    fcache = str(DataFileManager(engine).cache_file(fname))
     if cache and exists(fcache):
         relevant_if_metadata_above = (
             {"wavenum_max": load_wavenum_min} if load_wavenum_min else {}
@@ -312,6 +313,6 @@ def cdsd2df(
 
 
 if __name__ == "__main__":
-    from radis.test.test_io import test_hitemp
+    from radis.test.io.test_hitran_cdsd import _run_testcases
 
-    print("Testing cdsd: ", test_hitemp())
+    print("Testing cdsd: ", _run_testcases())
diff --git a/radis/io/dbmanager.py b/radis/io/dbmanager.py
index c61e12398..232e2f1ac 100644
--- a/radis/io/dbmanager.py
+++ b/radis/io/dbmanager.py
@@ -16,9 +16,9 @@
 
 try:
     from .cache_files import check_not_deprecated
-    from .hdf5 import HDF5Manager, hdf2df
+    from .hdf5 import DataFileManager
 except ImportError:
-    from radis.io.hdf5 import hdf2df, HDF5Manager
+    from radis.io.hdf5 import DataFileManager
     from radis.io.cache_files import check_not_deprecated
 
 from datetime import date
@@ -324,8 +324,10 @@ def register(self, dict_entries):
             self.verbose,
         )
 
-    def get_hdf5_manager(self):
-        return HDF5Manager(engine=self.engine)
+    def get_datafile_manager(self, engine=None):
+        if engine is None:
+            engine = self.engine
+        return DataFileManager(engine=engine)
 
     def download_and_parse(self, urlnames, local_files):
         all_local_files, _ = self.get_filenames()
@@ -425,10 +427,10 @@ def clean_download_files(self):
     def load(
         self,
         local_files,
-        isotope,
-        columns,
-        load_wavenum_min,
-        load_wavenum_max,
+        columns=None,
+        lower_bound=[],
+        upper_bound=[],
+        within=[],
         output="pandas",
     ):
         """
@@ -439,20 +441,31 @@ def load(
         output: 'pandas', 'vaex', 'jax'
             format of the output DataFrame. If ``'jax'``, returns a dictionary of
             jax arrays.
+        lower_bound: list of tuples [(column, lower_bound), etc.]
+            ::
+
+                lower_bound =[("wav", load_wavenum_min)]
+        upper_bound_bound: list of tuples [(column, upper_bound), etc.]
+            ::
+
+                upper_bound=[("wav", load_wavenum_max)]
+        within: list of tuples [(column, within_list), etc.]
+            ::
+
+                within=[("iso", isotope.split(","))]
         """
         engine = self.engine
-        if engine == "pytables":
+        mgr = self.get_datafile_manager()
+        if engine in ["pytables", "feather"]:
             df_all = []
             for local_file in local_files:
                 df_all.append(
-                    hdf2df(
+                    mgr.load(
                         local_file,
                         columns=columns,
-                        isotope=isotope,
-                        load_wavenum_min=load_wavenum_min,
-                        load_wavenum_max=load_wavenum_max,
-                        verbose=self.verbose,
-                        engine=engine,
+                        lower_bound=lower_bound,
+                        upper_bound=upper_bound,
+                        within=within,
                         output=output,
                     )
                 )
@@ -460,14 +473,12 @@ def load(
 
         elif engine == "vaex":
             # vaex can open several files at the same time:
-            return hdf2df(
+            return mgr.load(
                 local_files,
                 columns=columns,
-                isotope=isotope,
-                load_wavenum_min=load_wavenum_min,
-                load_wavenum_max=load_wavenum_max,
-                verbose=self.verbose,
-                engine=engine,
+                lower_bound=lower_bound,
+                upper_bound=upper_bound,
+                within=within,
                 output=output,
             )
         else:
@@ -503,6 +514,37 @@ def get_nrows(self, local_file):
             raise ValueError(engine)
         return nrows
 
+    def add_column(self, df, key, value):
+        """Create column ``key`` in DataFrame or dictionary ``df`` with value ``value``"""
+        from radis.misc.basics import is_number
+
+        if is_number(self.alpha_ref):
+            import vaex
+
+            if isinstance(df, vaex.dataframe.DataFrameLocal):
+                # see https://github.com/vaexio/vaex/pull/1570
+                df[key] = vaex.vconstant(float(value), length=len(df))
+            else:
+                df[key] = value
+        else:
+            df[key] = value
+
+    def rename_columns(self, df, rename_dict):
+        """Example::
+
+        mdb.rename_columns(df, {"nu_lines":"wav"})
+        """
+        import vaex
+
+        if isinstance(df, vaex.dataframe.DataFrameLocal):
+            for k, v in rename_dict.items():
+                df.rename(k, v)
+        elif isinstance(df, pd.DataFrame):
+            df.rename(columns=rename_dict, inplace=True)
+        elif isinstance(df, dict):
+            for k, v in rename_dict.items():
+                df[v] = df.pop(k)
+
 
 def register_database(databank_name, dict_entries, verbose):
     """Add to registered databases in RADIS config file.
diff --git a/radis/io/exomol.py b/radis/io/exomol.py
index 37d367d02..78228959c 100644
--- a/radis/io/exomol.py
+++ b/radis/io/exomol.py
@@ -6,23 +6,26 @@
 code (which you should also have a look at !), by @HajimeKawahara, under MIT License.
 
 """
+import os
 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 .dbmanager import DatabaseManager
     from .exomol_utils import e2s
 except ImportError:  # if local import
     from radis.io import exomolapi
+    from radis.io.dbmanager import DatabaseManager
     from radis.io.exomol_utils import e2s
 
-from radis.db.classes import get_molecule_identifier
+from radis.db.classes import (
+    EXOMOL_MOLECULES,
+    EXOMOL_ONLY_ISOTOPES_NAMES,
+    get_molecule_identifier,
+)
 
 EXOMOL_URL = "http://www.exomol.com/db/"
 
@@ -51,7 +54,10 @@ def get_exomol_full_isotope_name(molecule, isotope):
     --------
     :py:func:`~radis.io.exomol.get_exomol_database_list`"""
 
-    from radis.db.classes import EXOMOL_ONLY_ISOTOPES_NAMES
+    if molecule not in EXOMOL_MOLECULES:
+        raise ValueError(
+            f"Molecule {molecule} not available in ExoMol. Change database, or choose one of ExoMol available molecules: {sorted(EXOMOL_MOLECULES)}"
+        )
 
     if (molecule, isotope) in EXOMOL_ONLY_ISOTOPES_NAMES:
         return EXOMOL_ONLY_ISOTOPES_NAMES[(molecule, isotope)]
@@ -134,7 +140,12 @@ def get_exomol_database_list(molecule, isotope_full_name):
 
     databases = databases + databases_recommended
 
-    return databases, databases_recommended[0]
+    if len(databases_recommended) > 0:
+        recommended_database = databases_recommended[0]
+    else:
+        recommended_database = False
+
+    return databases, recommended_database
 
 
 # def fetch_exomol_molecule_list():
@@ -226,6 +237,8 @@ def fetch_exomol(
     return_local_path=False,
     return_partition_function=False,
     engine="default",
+    output="pandas",
+    skip_optional_data=True,
 ):
     """Stream ExoMol file from EXOMOL website. Unzip and build a HDF5 file directly.
 
@@ -270,9 +283,9 @@ def fetch_exomol(
 
     Other Parameters
     ----------------
-    cache: bool, or ``'regen'``
+    cache: bool, or ``'regen'`` or ``'force'``
         if ``True``, use existing HDF5 file. If ``False`` or ``'regen'``, rebuild it.
-        Default ``True``.
+        If ``'force'``, crash if not cache file found. Default ``True``.
     verbose: bool
     clean_cache_files: bool
         if ``True`` clean downloaded cache files after HDF5 are created.
@@ -282,6 +295,25 @@ def fetch_exomol(
         if ``True``, also returns a :py:class:`~radis.levels.partfunc.PartFuncExoMol` object.
     engine: 'vaex', 'feather'
         which memory-mapping library to use. If 'default' use the value from ~/radis.json
+    output: 'pandas', 'vaex', 'jax'
+        format of the output DataFrame. If ``'jax'``, returns a dictionary of
+        jax arrays.
+    skip_optional_data : bool
+        If False, fetch all fields which are marked as available in the ExoMol definition
+        file. If True, load only the first 4 columns of the states file
+        ("i", "E", "g", "J"). The structure of the columns above 5 depend on the
+        the definitions file (*.def) and the Exomol version.
+        If ``skip_optional_data=False``, two errors may occur:
+
+            - a field is marked as present/absent in the *.def field but is
+              absent/present in the *.states file (ie both files are inconsistent).
+            - in the updated version of Exomol, new fields have been added in the
+              states file of some species. But it has not been done for all species,
+              so both structures exist. For instance, the states file of
+              https://exomol.com/data/molecules/HCl/1H-35Cl/HITRAN-HCl/ follows the
+              structure described in [1]_, unlike the states file of
+              https://exomol.com/data/molecules/NO/14N-16O/XABC/ which follows the
+              structure described in [2]_.
 
     Returns
     -------
@@ -306,6 +338,12 @@ def fetch_exomol(
     the expected wavenumber range & isotopes are returned. The .HFD5 parsing uses
     :py:func:`~radis.io.hdf5.hdf2df`
 
+    References
+    ----------
+
+    .. [1] Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., Barton, E. J., Chubb, K. L., Coles, P. A., … Zak, E. (2016). The ExoMol database: molecular line lists for exoplanet and other hot atmospheres. https://doi.org/10.1016/j.jms.2016.05.002
+    .. [2] Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., Clark, V. H. J., Chubb, K. L., Conway, E. K., … Yurchenko, O. P. (2020). The 2020 release of the ExoMol database: Molecular line lists for exoplanet and other hot atmospheres. Journal of Quantitative Spectroscopy and Radiative Transfer, 255, 107228. https://doi.org/10.1016/j.jqsrt.2020.107228
+
     See Also
     --------
     :py:func:`~radis.io.hitran.fetch_hitran`, :py:func:`~radis.io.hitemp.fetch_hitemp`
@@ -363,26 +401,80 @@ def fetch_exomol(
         / database
     )
 
-    if load_wavenum_min is None:
-        load_wavenum_min = -np.inf
-    if load_wavenum_max is None:
-        load_wavenum_max = np.inf
-    mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max], engine=engine)
+    # TODO: add deprecation if missing columns in cache file
+
+    # Init database, download files if needed.
+    mdb = MdbExomol(
+        local_path,
+        molecule=molecule,
+        name=databank_name,
+        local_databases=local_databases,
+        nurange=[
+            load_wavenum_min if load_wavenum_min is not None else 0.0,
+            load_wavenum_max if load_wavenum_max is not None else np.inf,
+        ],
+        engine=engine,
+        cache=cache,
+        skip_optional_data=skip_optional_data,
+    )
+
+    # Get local files
+    local_files = mdb.trans_file
+    if not isinstance(local_files, list):
+        local_files = [local_files]
+    mgr = mdb.get_datafile_manager()
+    local_files = [mgr.cache_file(f) for f in local_files]
+
+    # Specific for RADIS : rename columns
+    radis2exomol_columns = {
+        "wav": "nu_lines",
+        "airbrd": "alpha_ref",
+        "El": "elower",
+        "ju": "jupper",
+        "jl": "jlower",
+        "gp": "gupper",
+        "gpp": "glower",
+        "Tdpair": "n_Texp",
+    }
+    # get column name converting to exomol/exojax format if possible, else use the same
+    if columns is not None:
+        columns_exomol = [radis2exomol_columns.get(c, c) for c in columns] + [
+            "Sij0",
+            "jlower",
+            "jupper",
+        ]  # needed for broadening
+    else:
+        columns_exomol = None
+
+    df = mdb.load(
+        local_files,
+        columns=columns_exomol,
+        lower_bound=([("nu_lines", load_wavenum_min)] if load_wavenum_min else [])
+        + ([("Sij0", mdb.crit)] if not np.isneginf(mdb.crit) else []),
+        upper_bound=([("nu_lines", load_wavenum_max)] if load_wavenum_max else []),
+        output=output,
+    )
+
+    if "jlower" not in df:
+        raise KeyError(
+            f"jlower not found. Maybe try to delete cache file {local_files} and restart?"
+        )
 
-    # Attributes of the DataFrame
-    from radis.db.classes import HITRAN_MOLECULES
+    # Add broadening
+    mdb.set_broadening(df, output=output)
 
-    attrs = {}
-    if molecule in HITRAN_MOLECULES:
-        attrs["id"] = get_molecule_identifier(
-            molecule
-        )  # HITRAN id-number (if available)
-    attrs["molecule"] = molecule
-    attrs["iso"] = isotope
+    # Specific for RADIS :
+    # ... Get RADIS column names:
+    exomol2radis_columns = {v: k for k, v in radis2exomol_columns.items()}
+    mdb.rename_columns(df, {k: v for k, v in exomol2radis_columns.items() if k in df})
 
-    df = mdb.to_df(attrs=attrs)
+    assert "wav" in df
 
-    # Replace Linestrength with Line Intensity taking into account Terrestrial isotopic abundance
+    # ... include isotopic abundance in linestrength :
+    # Note : ExoMol treats isotopes as independant molecules ; linestrength is not
+    # corrected by isotopic abundance.
+    # Below, replace Linestrength with Line Intensity taking into account
+    # Terrestrial isotopic abundance (to be compatible with HITRAN/HITEMP/etc. )
     from radis.db.molparam import MolParams
 
     try:
@@ -392,12 +484,35 @@ def fetch_exomol(
 
         Ia = config["molparams"]["abundances"][molecule][str(isotope)]
 
-    df["Sij0"] *= Ia
-    df.rename(columns={"Sij0": "int"}, inplace=True)
+    if output == "jax":
+        try:
+            import jax.numpy as jnp
+        except:
+            import numpy as jnp
+        df["logsij0"] += jnp.log(Ia)
+    else:
+        df["Sij0"] *= Ia
+        mdb.rename_columns(df, {"Sij0": "int"})
 
-    # Return:
-    out = [df]
+    # Add Attributes of the DataFrame
+    if output == "pandas":  # no attribtes in "Jax" or "Vaex" mode
+        from radis.db.classes import HITRAN_MOLECULES
+
+        attrs = {}
+        if molecule in HITRAN_MOLECULES:
+            attrs["id"] = get_molecule_identifier(
+                molecule
+            )  # HITRAN id-number (if available)
+        attrs["molecule"] = molecule
+        attrs["iso"] = isotope
 
+        for k, v in attrs.items():
+            df.attrs[k] = v
+
+    # Return:
+    out = df
+    if return_local_path or return_partition_function:
+        out = [out]
     if return_local_path:
         out.append(str(mdb.path))
     if return_partition_function:
@@ -407,30 +522,7 @@ def fetch_exomol(
     return out
 
 
-class MdbExomol(object):
-    """molecular database of ExoMol
-
-    MdbExomol is a class for ExoMol.
-
-    Parameters
-    ----------
-    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
-
-    .. minigallery:: radis.fetch_exomol
-
-    """
+class MdbExomol(DatabaseManager):
 
     # TODO : inherit from DatabaseManager or similar
 
@@ -441,15 +533,23 @@ class MdbExomol(object):
     def __init__(
         self,
         path,
-        nurange=[-np.inf, np.inf],
+        molecule,
+        database=None,
+        local_databases=None,
+        name="EXOMOL-{molecule}",
+        nurange=[0.0, np.inf],
         margin=1.0,
         crit=-np.inf,
-        bkgdatm="H2",  # TODO: use Air whenever possible, to be consistent with HITRAN/HITEMP
+        bkgdatm="Air",  # TODO: use Air whenever possible, to be consistent with HITRAN/HITEMP
         broadf=True,
         engine="vaex",
         verbose=True,
+        cache=True,
+        skip_optional_data=True,
     ):
-        """Molecular database for Exomol form
+        """molecular database of ExoMol
+
+        MdbExomol is a class for ExoMol.
 
         Parameters
         ----------
@@ -460,6 +560,27 @@ def __init__(
         bkgdatm: background atmosphere for broadening. e.g. H2, He,
         broadf: if False, the default broadening parameters in .def file is used
 
+        Other Parameters
+        ----------------
+        engine : str
+            which memory mapping engine to use : 'vaex', 'pytables' (HDF5), 'feather'
+        skip_optional_data : bool
+            If False, fetch all fields which are marked as available in the ExoMol definition
+            file. If True, load only the first 4 columns of the states file
+            ("i", "E", "g", "J"). The structure of the columns above 5 depend on the
+            the definitions file (*.def) and the Exomol version.
+            If ``skip_optional_data=False``, two errors may occur:
+
+                - a field is marked as present/absent in the *.def field but is
+                  absent/present in the *.states file (ie both files are inconsistent).
+                - in the updated version of Exomol, new fields have been added in the
+                  states file of some species. But it has not been done for all species,
+                  so both structures exist. For instance, the states file of
+                  https://exomol.com/data/molecules/HCl/1H-35Cl/HITRAN-HCl/ follows the
+                  structure described in [1]_, unlike the states file of
+                  https://exomol.com/data/molecules/NO/14N-16O/XABC/ which follows the
+                  structure described in [2]_.
+
         Notes
         -----
 
@@ -467,7 +588,70 @@ def __init__(
         we convert it to the feather or hdf5-format. After the second-time,
         we use the feather/hdf5 format instead.
 
+        Examples
+        --------
+        ::
+
+            # Init database, download files if needed.
+            mdb = MdbExomol(
+                local_path,
+                molecule=molecule,
+                name=databank_name,
+                local_databases=local_databases,
+                # nurange=[load_wavenum_min, load_wavenum_max],
+                engine="vaex",
+            )
+
+            # Get cache files to load :
+            mgr = mdb.get_datafile_manager()
+            local_files = [mgr.cache_file(f) for f in mdb.trans_file]
+
+            # Load files
+            df = mdb.load(
+                local_files,
+                columns=columns_exomol,
+                lower_bound=([('nu_lines', load_wavenum_min)] if load_wavenum_min else []) + ([("Sij0", mdb.crit)] if not np.isneginf(mdb.crit) else []),
+                upper_bound=([('nu_lines', load_wavenum_max)] if load_wavenum_max else []),
+                output="jax", # or "pytables", "vaex"
+            )
+
+        .. minigallery:: radis.fetch_exomol
+
+
+        DataFrame columns
+        -----------------
+
+        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
+
+        References
+        ----------
+
+        .. [1] Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., Barton, E. J., Chubb, K. L., Coles, P. A., … Zak, E. (2016). The ExoMol database: molecular line lists for exoplanet and other hot atmospheres. https://doi.org/10.1016/j.jms.2016.05.002
+        .. [2] Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., Clark, V. H. J., Chubb, K. L., Conway, E. K., … Yurchenko, O. P. (2020). The 2020 release of the ExoMol database: Molecular line lists for exoplanet and other hot atmospheres. Journal of Quantitative Spectroscopy and Radiative Transfer, 255, 107228. https://doi.org/10.1016/j.jqsrt.2020.107228
+
         """
+        super().__init__(
+            name,
+            molecule,
+            local_databases,
+            engine,
+            verbose=verbose,
+        )
+
+        assert cache  # cache only used for cache='regen' or cache='force' modes, cache=False is not expected
+
         from os import environ
 
         if engine == "default":
@@ -487,18 +671,11 @@ def __init__(
                 else:
                     engine = "vaex"
 
-        if engine == "vaex":
-            import vaex
-        elif engine == "feather":
-            pass
-            # vaex will be the future default, but it still fails on some systems
-            # Ex : on Spyder : See https://github.com/spyder-ide/spyder/issues/16183
-            # We keep "feather" as backup)
-        else:
-            raise NotImplementedError(f"{engine} is not implemented")
-        self.engine = engine  # which memory mapping engine to use : 'vaex', 'pytables' (HDF5), 'feather'
-
         self.path = pathlib.Path(path)
+
+        if local_databases is not None:
+            self.path = pathlib.Path(local_databases).expanduser() / self.path
+
         t0 = self.path.parents[0].stem
         molec = t0 + "__" + str(self.path.stem)
         self.bkgdatm = bkgdatm
@@ -508,6 +685,7 @@ def __init__(
         self.crit = crit
         self.margin = margin
         self.nurange = [np.min(nurange), np.max(nurange)]
+        self.wmin, self.wmax = self.nurange
         self.broadf = broadf
         # Where exomol files are
         self.states_file = self.path / pathlib.Path(molec + ".states.bz2")
@@ -528,62 +706,57 @@ def __init__(
         tag = molec.split("__")
         self.isotope_fullname = tag[0]
         self.molecule = e2s(tag[0])
-        self.isotope = 1  # Placeholder. TODO : impement parsing of other isotopes.
+        # 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:
+            self.n_Texp_def = 0.5
             warnings.warn(
                 Warning(
-                    """
-                    No default broadening exponent in def file. Assigned n = 0.5
+                    f"""
+                    No default broadening exponent in def file. Assigned n = {self.n_Texp_def}
                     """
                 )
             )
-            self.n_Texp_def = 0.5
         #  default alpha_ref value if not given
         if self.alpha_ref_def is None:
+            self.alpha_ref_def = 0.07
             warnings.warn(
                 Warning(
-                    """
-                    No default broadening in def file. Assigned alpha_ref = 0.07
+                    f"""
+                    No default broadening in def file. Assigned alpha_ref = {self.alpha_ref_def}
                     """
                 )
             )
-            self.alpha_ref_def = 0.07
 
         # load states
-        if engine == "feather":
-            if self.states_file.with_suffix(".feather").exists():
-                states = pd.read_feather(self.states_file.with_suffix(".feather"))
-            else:
-                print(
-                    "Note: Caching states data to the feather format. After the second time, it will become much faster."
-                )
-                states = exomolapi.read_states(self.states_file, dic_def, engine="csv")
-                states.to_feather(self.states_file.with_suffix(".feather"))
-            ndstates = states.to_numpy()[
-                :, :4
-            ]  # the i, E, g, J are in the 4 first columns
-        elif engine == "vaex":
-            if self.states_file.with_suffix(".bz2.hdf5").exists():
-                states_file_cache = self.states_file.with_suffix(".bz2.hdf5")
-                states = vaex.open(states_file_cache)
-            elif self.states_file.with_suffix(".hdf5").exists():
-                states_file_cache = self.states_file.with_suffix(".hdf5")
-                states = vaex.open(states_file_cache)
-            else:
-                print(
-                    "Note: Caching states data to the hdf5 format with vaex. After the second time, it will become much faster."
+        mgr = self.get_datafile_manager()
+        if cache == "regen" and mgr.cache_file(self.states_file).exists():
+            if verbose:
+                print("Removing existing file ", mgr.cache_file(self.states_file))
+            os.remove(mgr.cache_file(self.states_file))
+        if mgr.cache_file(self.states_file).exists():
+            states = mgr.read(mgr.cache_file(self.states_file))
+        else:
+            if cache == "force":
+                raise ValueError(
+                    f"Cache file {str(mgr.cache_file(self.states_file))} does not exist"
                 )
-                states = exomolapi.read_states(self.states_file, dic_def, engine="vaex")
-            ndstates = vaex.array_types.to_numpy(states)
+            print(
+                f"Note: Caching states data to the {engine} format. After the second time, it will become much faster."
+            )
+            states = exomolapi.read_states(
+                self.states_file,
+                dic_def,
+                engine="vaex" if engine == "vaex" else "csv",
+                skip_optional_data=skip_optional_data,
+            )
+            mgr.write(mgr.cache_file(self.states_file), states)
 
         # load pf
         pf = exomolapi.read_pf(self.pf_file)
@@ -597,354 +770,109 @@ def __init__(
         self.Tref = 296.0
         self.QTref = np.array(self.QT_interp(self.Tref))
 
+        # Download files
+
+        # Generate list of files
+        # ---------------------
         # Case1 : transitions are stored in a single 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 engine == "vaex":
-                if self.trans_file.with_suffix(".hdf5").exists():
-                    trans = vaex.open(self.trans_file.with_suffix(".hdf5"))
-                    # Apply filter
-                    cdt = 1
-                    if not np.isneginf(self.nurange[0]):
-                        cdt *= trans.nu_lines > self.nurange[0] - self.margin
-                    if not np.isinf(self.nurange[1]):
-                        cdt *= trans.nu_lines < self.nurange[1] - self.margin
-                    if not np.isneginf(self.crit):
-                        cdt *= trans.Sij0 > self.crit
-                    if cdt != 1:
-                        trans = trans[cdt]
-
-                    # mask has been alraedy applied
-                    mask_needed = False
-
-                else:
-                    print(
-                        "Note: Caching line transition data to the HDF5 format with vaex. After the second time, it will become much faster."
-                    )
-                    trans = exomolapi.read_trans(self.trans_file, engine="vaex")
-
-                    # Check validity
-                    b = trans.nu_lines.isnan()  # TODO: check if can be made faster?
-                    if b.sum() > 0:
-                        raise NotImplementedError(
-                            "transition wavenumber not given. Has to be computed from upper/lower level energy look up"
-                        )
-                    # mask needs to be applied
-                    mask_needed = True
-
-                # TODO Load extra columns (precomputed Sij0, etc.)
-                # trans_file_extra = self.trans_file.with_name(
-                #    self.trans_file.stem + "_extra" + self.trans_file.suffix
-                # )
-                # if trans_file_extra.exists():
-                #    trans_extra = vaex.open(trans_file_extra)
-                #    # trans = trans.join(trans_extra)  # merge without column name will merge on rows basis
-
-                ndtrans = vaex.array_types.to_numpy(trans)
-
-            elif engine == "feather":
-                if self.trans_file.with_suffix(".feather").exists():
-                    trans = pd.read_feather(self.trans_file.with_suffix(".feather"))
-                else:
-                    print(
-                        "Note: Caching line transition data to the feather format. After the second time, it will become much faster."
-                    )
-                    trans = exomolapi.read_trans(self.trans_file, engine="csv")
-                    trans.to_feather(self.trans_file.with_suffix(".feather"))
-                ndtrans = trans.to_numpy()
-                # mask needs to be applied   (in feather mode we don't sleect wavneumbers)
-                mask_needed = True
-
-            # compute extra parameters (Sij0, gup and elower, etc.)
-            (
-                self._A,
-                self.nu_lines,
-                self._elower,
-                self._gpp,
-                self._jlower,
-                self._jupper,
-                mask_zeronu,
-                self._quantumNumbers,
-            ) = exomolapi.pickup_gE(ndstates, ndtrans, self.trans_file, dic_def)
-
-            if engine == "vaex" and self.trans_file.with_suffix(".hdf5").exists():
-                # if engine == 'feather' we recompute all the time
-                # Todo : we should get the column name instead ?
-                self.Sij0 = ndtrans[:, 4]  # TODO: Refactor remove array index lookups.
-                # TODO Refactor. Add extra here:   trans = trans.join(trans_extra)
-            else:
-                ##Line strength:
-                from radis.lbl.base import (  # TODO: move elsewhere
-                    linestrength_from_Einstein,
-                )
-
-                self.Sij0 = linestrength_from_Einstein(
-                    A=self._A,
-                    gu=self._gpp,
-                    El=self._elower,
-                    Ia=1,  #  Sij0 is a linestrength calculated without taking into account isotopic abundance (unlike line intensity parameter of HITRAN. In RADIS this is corrected for in fetch_exomol()  )
-                    nu=self.nu_lines,
-                    Q=self.QTref,
-                    T=self.Tref,
-                )
-
-                trans["nu_positive"] = mask_zeronu
-                if engine == "vaex":
-                    # exclude the lines whose nu_lines evaluated inside exomolapi.pickup_gE (thus sometimes different from the "nu_lines" column in trans) is not positive
-
-                    trans = trans[trans.nu_positive].extract()
-                    trans.drop("nu_positive", inplace=True)
-                else:
-                    if False in mask_zeronu:
-                        raise NotImplementedError(
-                            "some wavenumber is not defined;  masking not impleemtend so far in 'feather' engine"
-                        )
-
-                trans["nu_lines"] = self.nu_lines
-                trans["Sij0"] = self.Sij0
-
-                if engine == "vaex":
-                    trans.export(self.trans_file.with_suffix(".hdf5"))
-                #  TODO : implement masking in 'feather' mode
+            self.trans_file = [self.path / pathlib.Path(molec + ".trans.bz2")]
+            self.num_tag = [None]
 
         # Case2 : Transitions are stored in multiple files:
         else:  # dic_def["numinf"] is not None
-            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
+            # dic_def["numinf"] contains the limit points ``[w(0), w(1), ..., w(n)]``
+            # (n+1 elements) defining the spectral ranges appearing in the list
+            # dic_def["numtag"] which looks like ``["w(0)-w(1)", "w(1)-w(2)", ..., w(n-1)-w(n)]``
+            # (n elements)
+            # imin :
+            # index i of the "w(i)-w(i+1)" element in dic_def["numtag"] such
+            # that nurange[0]<=w(i)
+            imin = np.searchsorted(dic_def["numinf"], nurange[0], side="right") - 1
+            # imax :
+            # index i of the "w(i)-w(i+1)" element in dic_def["numtag"] such
+            # that w(i+1)<=nurange[1]
+            imax = np.searchsorted(dic_def["numinf"], nurange[1], side="right") - 2
             self.trans_file = []
+            self.num_tag = []
             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 engine == "feather":
-                    trans_feather_path = trans_file.with_suffix(".bz2.feather")
-                    if trans_feather_path.exists():
-                        trans = pd.read_feather(trans_feather_path)
-                    else:
-                        print(
-                            "Note: Caching line transition data to the feather format. After the second time, it will become much faster."
-                        )
-                        trans = exomolapi.read_trans(trans_file, engine="csv")
-                        trans.to_feather(trans_feather_path)
-                        #!!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)]
-                    ndtrans = trans.to_numpy()
-                    # mask needs to be applied   (in feather mode we don't sleect wavneumbers)
-                    mask_needed = True
-                elif engine == "vaex":
-                    if trans_file.with_suffix(".hdf5").exists():
-                        trans = vaex.open(trans_file.with_suffix(".hdf5"))
-                        cdt = 1
-                        if not np.isneginf(self.nurange[0]):
-                            cdt *= trans.nu_lines > self.nurange[0] - self.margin
-                        if not np.isinf(self.nurange[1]):
-                            cdt *= trans.nu_lines < self.nurange[1] - self.margin
-                        if not np.isneginf(self.crit):
-                            cdt *= trans.Sij0 > self.crit
-                        if cdt != 1:
-                            trans = trans[cdt]
-
-                        # mask has been already applied
-                        mask_needed = False
-                    else:
-                        print(
-                            "Note: Caching line transition data to the HDF5 format with vaex. After the second time, it will become much faster."
-                        )
-                        trans = exomolapi.read_trans(trans_file, engine="vaex")
-
-                        # Check validity
-                        b = trans.nu_lines.isnan()  # TODO: check if can be made faster?
-                        if b.sum() > 0:
-                            raise NotImplementedError(
-                                "transition wavenumber not given. Has to be computed from upper/lower level energy look up"
-                            )
-
-                        # TODO Load extra columns (precomputed Sij0, etc.)
-                        # trans_file_extra = self.trans_file.with_name(
-                        #    self.trans_file.stem + "_extra" + self.trans_file.suffix
-                        # )
-                        # if trans_file_extra.exists():
-                        #    trans_extra = vaex.open(trans_file_extra)
-                        #    trans = trans.join(
-                        #        trans_extra
-                        #    )  # merge without column name will merge on rows basis
-
-                        # mask needs to be applied
-                        mask_needed = True
-
-                    ndtrans = vaex.array_types.to_numpy(trans)
                 self.trans_file.append(trans_file)
+                self.num_tag.append(dic_def["numtag"][i])
 
-                # Complete transition data with lookup on upper & lower state :
-                # In particular, compute gup and elower
-                if k == 0:
-                    (
-                        self._A,
-                        self.nu_lines,
-                        self._elower,
-                        self._gpp,
-                        self._jlower,
-                        self._jupper,
-                        mask_zeronu,
-                        self._quantumNumbers,
-                    ) = exomolapi.pickup_gE(ndstates, ndtrans, trans_file, dic_def)
-
-                    if engine == "vaex" and trans_file.with_suffix(".hdf5").exists():
-                        # Sij0x already stored in cache file
-                        self.Sij0 = ndtrans[:, 4]
-                    else:
-                        ##Recompute Line strength:
-                        from radis.lbl.base import (  # TODO: move elsewhere
-                            linestrength_from_Einstein,
-                        )
+        # Look-up missing parameters and write file
+        # -----------------------------------------
+        for trans_file, num_tag in zip(self.trans_file, self.num_tag):
 
-                        self.Sij0 = linestrength_from_Einstein(
-                            A=self._A,
-                            gu=self._gpp,
-                            El=self._elower,
-                            Ia=1,  #  Sij0 is a linestrength calculated without taking into account isotopic abundance (unlike line intensity parameter of HITRAN. In RADIS this is corrected for in fetch_exomol()  )
-                            nu=self.nu_lines,
-                            Q=self.QTref,
-                            T=self.Tref,
-                        )
+            if cache == "regen" and mgr.cache_file(trans_file).exists():
+                if verbose:
+                    print("Removing existing file ", mgr.cache_file(trans_file))
+                os.remove(mgr.cache_file(trans_file))
 
-                        # exclude the lines whose nu_lines evaluated inside exomolapi.pickup_gE (thus sometimes different from the "nu_lines" column in trans) is not positive
-                        if engine == "vaex":
-                            trans["nu_positive"] = mask_zeronu
-                            trans = trans[trans.nu_positive].extract()
-                            trans.drop("nu_positive", inplace=True)
-                        else:
-                            if False in mask_zeronu:
-                                raise NotImplementedError(
-                                    "some wavenumber is not defined;  masking not impleemtend so far in 'feather' engine"
-                                )
-
-                        trans["nu_lines"] = self.nu_lines
-                        trans["Sij0"] = self.Sij0
-                else:  # k!=0
-
-                    (
-                        Ax,
-                        nulx,
-                        elowerx,
-                        gppx,
-                        jlowerx,
-                        jupperx,
-                        mask_zeronu,
-                        quantumNumbersx,
-                    ) = exomolapi.pickup_gE(ndstates, ndtrans, trans_file, dic_def)
-                    if engine == "vaex" and trans_file.with_suffix(".hdf5").exists():
-                        # Sij0x already stored in cache file
-                        Sij0x = ndtrans[:, 4]
-                    else:
-                        ##Recompute Line strength:
-                        from radis.lbl.base import (  # TODO: move elsewhere
-                            linestrength_from_Einstein,
-                        )
+            if not mgr.cache_file(trans_file).exists():
 
-                        Sij0x = linestrength_from_Einstein(
-                            A=Ax,
-                            gu=gppx,
-                            El=elowerx,
-                            Ia=1,  #  Sij0 is a linestrength calculated without taking into account isotopic abundance (unlike line intensity parameter of HITRAN. In RADIS this is corrected for in fetch_exomol()  )
-                            nu=nulx,
-                            Q=self.QTref,
-                            T=self.Tref,
-                        )
+                if cache == "force":
+                    raise ValueError(
+                        f"Cache file {str(mgr.cache_file(trans_file))} does not exist"
+                    )
 
-                        # exclude the lines whose nu_lines evaluated inside exomolapi.pickup_gE (thus sometimes different from the "nu_lines" column in trans) is not positive
-                        if engine == "vaex":
-                            trans["nu_positive"] = mask_zeronu
-                            trans = trans[trans.nu_positive].extract()
-                            trans.drop("nu_positive", inplace=True)
-                        else:
-                            if False in mask_zeronu:
-                                raise NotImplementedError(
-                                    "some wavenumber is not defined;  masking not impleemtend so far in 'feather' engine"
-                                )
-
-                        trans["nu_lines"] = nulx
-                        trans["Sij0"] = Sij0x
-
-                    self._A = np.hstack([self._A, Ax])
-                    self.nu_lines = np.hstack([self.nu_lines, nulx])
-                    self.Sij0 = np.hstack([self.Sij0, Sij0x])
-                    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])
-                    for key in self._quantumNumbers.keys():
-                        self._quantumNumbers[key] = np.hstack(
-                            [self._quantumNumbers[key], quantumNumbersx[key]]
-                        )
+                if not trans_file.exists():
+                    self.download(molec, extension=[".trans.bz2"], numtag=num_tag)
+                    # TODO: add option to delete file at the end
 
-                    #  TODO refactpr : store into "extra" file ; create column if needed ;if it did exist already make sure the value was computed for ALL lines
-                    if (
-                        engine == "vaex"
-                        and not trans_file.with_suffix(".hdf5").exists()
-                    ):
-                        trans.export(trans_file.with_suffix(".hdf5"))
-
-        ### MASKING ###
-        mask = (
-            (self.nu_lines > self.nurange[0] - self.margin)
-            * (self.nu_lines < self.nurange[1] + self.margin)
-            * (self.Sij0 > self.crit)
-        )
+                print(
+                    f"Note: Caching line transition data to the {engine} format. After the second time, it will become much faster."
+                )
+                trans = exomolapi.read_trans(
+                    trans_file, engine="vaex" if engine == "vaex" else "csv"
+                )
 
-        self.masking(mask, mask_needed)
+                # Complete transition data with lookup on upper & lower state :
+                # In particular, compute gup and elower
 
-    def masking(self, mask, mask_needed=True):
-        """applying mask and (re)generate jnp.arrays
+                exomolapi.pickup_gE(
+                    states,
+                    trans,
+                    trans_file,
+                    dic_def,
+                    skip_optional_data=skip_optional_data,
+                )
 
-        Args:
-           mask: mask to be applied. self.mask is updated.
-           mask_needed: whether mask needs to be applied or not
+                ##Recompute Line strength:
+                from radis.lbl.base import (  # TODO: move elsewhere
+                    linestrength_from_Einstein,
+                )
 
-        """
-        # TODO : with Vaex-HDF5, selection should happen on disk.
-
-        # numpy float 64 Do not convert them jnp array
-        if mask_needed:
-            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]
-            for key, array in self._quantumNumbers.items():
-                self._quantumNumbers[key] = array[mask]
-
-        # jnp arrays
-        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()
-
-    def set_broadening(self, alpha_ref_def=None, n_Texp_def=None):
+                self.Sij0 = linestrength_from_Einstein(
+                    A=trans["A"],
+                    gu=trans["gup"],
+                    El=trans["elower"],
+                    Ia=1,  #  Sij0 is a linestrength calculated without taking into account isotopic abundance (unlike line intensity parameter of HITRAN. In RADIS this is corrected for in fetch_exomol()  )
+                    nu=trans["nu_lines"],
+                    Q=self.QTref,
+                    T=self.Tref,
+                )
+
+                trans["Sij0"] = self.Sij0
+
+                mgr.write(mgr.cache_file(trans_file), trans)
+
+        # Database ready to be loaded.
+        # Proceed with mdb.load()
+
+    def set_broadening(self, df, alpha_ref_def=None, n_Texp_def=None, output=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
+        Parameters
+        ----------
+        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
+
+        Returns
+        -------
+        None. Updates ``df`` with ``alpha_ref``, ``n_Texp_def``
         """
         if alpha_ref_def:
             self.alpha_ref_def = alpha_ref_def
@@ -962,38 +890,42 @@ def set_broadening(self, alpha_ref_def=None, n_Texp_def=None):
                         bdat,
                         alpha_ref_default=self.alpha_ref_def,
                         n_Texp_default=self.n_Texp_def,
-                        jlower_max=np.max(self._jlower),
+                        jlower_max=np.max(df["jlower"]),
                     )
-                    self.alpha_ref = np.array(j2alpha_ref[self._jlower])
-                    self.n_Texp = np.array(j2n_Texp[self._jlower])
+                    self.alpha_ref = np.array(j2alpha_ref[df["jlower"]])
+                    self.n_Texp = np.array(j2n_Texp[df["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),
+                        jlower_max=np.max(df["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),
+                        jupper_max=np.max(df["jupper"]),
                     )
-                    self.alpha_ref = np.array(jj2alpha_ref[self._jlower, self._jupper])
-                    self.n_Texp = np.array(jj2n_Texp[self._jlower, self._jupper])
-            except:
+                    self.alpha_ref = np.array(jj2alpha_ref[df["jlower"], df["jupper"]])
+                    self.n_Texp = np.array(jj2n_Texp[df["jlower"], df["jupper"]])
+            except FileNotFoundError:
                 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.alpha_ref_def * np.ones_like(df["jlower"])
                 )
-                self.n_Texp = np.array(self.n_Texp_def * np.ones_like(self._jlower))
+                self.n_Texp = np.array(self.n_Texp_def * np.ones_like(df["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))
+            self.alpha_ref = np.array(self.alpha_ref_def * np.ones_like(df["jlower"]))
+            self.n_Texp = np.array(self.n_Texp_def * np.ones_like(df["jlower"]))
+
+        # Add values
+        self.add_column(df, "alpha_ref", self.alpha_ref)
+        self.add_column(df, "n_Texp", self.n_Texp)
 
     def QT_interp(self, T):
         """interpolated partition function
@@ -1064,64 +996,15 @@ def download(self, molec, extension, numtag=None):
                 print("Downloading " + pfpath)
                 try:
                     urllib.request.urlretrieve(pfpath, str(self.path / pfname))
-                except:
-                    print("Error: Couldn't download " + ext + " file and save.")
+                except HTTPError:
+                    print(f"Error: Couldn't download {ext} file at {pfpath} and save.")
 
     def to_partition_function_tabulator(self):
+        """Generate a :py:class:`~radis.levels.partfunc.PartFuncExoMol` object"""
         from radis.levels.partfunc import PartFuncExoMol
 
         return PartFuncExoMol(self.isotope_fullname, self.T_gQT, self.gQT)
 
-    def to_df(self, attrs={}):
-        """Export Line Database to a RADIS-friendly Pandas DataFrame
-
-        Parameters
-        ----------
-        attrs: dict
-            add these as attributes of the 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,
-                "Sij0": self.Sij0,  #  linestrength (not corrected for isotopic abundance)
-                "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
-            }
-        )
-
-        # Check format :
-        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")
-
-        for k, v in attrs.items():
-            df.attrs[k] = v
-
-        assert self.Tref == 296  # default in RADIS
-
-        return df
-
-    def to_jax(self):
-        """Generate Jax arrays for use in ExoJax"""
-
-        pass  # TODO
-
 
 if __name__ == "__main__":
     # mdb=MdbExomol("/home/kawahara/exojax/data/CO/12C-16O/Li2015/")
@@ -1146,6 +1029,172 @@ def to_jax(self):
     # # s_hit = sf.eq_spectrum(500, name='HITRAN')
 
     #%% Test by direct caclulation
-    import pytest
+    # import pytest
+
+    # print("Testing factory:", pytest.main(["../test/io/test_exomol.py"]))
+
+    #%% RADIS-like Example
+    # uses fetch_exomol() internally
+
+    from radis import calc_spectrum
+
+    s = calc_spectrum(
+        wavelength_min=1.630e4,
+        wavelength_max=1.6305e4,
+        molecule="CH4",
+        isotope="1",
+        pressure=1.01325,  # bar
+        Tgas=1000,  # K
+        mole_fraction=0.1,
+        path_length=1,  # cm
+        # broadening_method="fft",  # @ dev: Doesn't work with 'voigt'
+        databank=(
+            "exomol",
+            "YT10to10",
+        ),  # Simply use 'exomol' for the recommended database
+    )
+    # s.apply_slit(1, "cm-1")  # simulate an experimental slit
+    s.plot("xsection")
+
+    # %% Exojax like Example
+
+    class mdbExoMol:
+
+        # hardcode attribute names, to prevent typos and the declaration of unwanted parameters
+        __slots__ = [
+            "Sij0",
+            "logsij0",
+            "nu_lines",
+            "A",
+            "elower",
+            "eupper",
+            "gupper",
+            "jlower",
+            "jupper",
+        ]
+
+        def __init__(
+            self,
+            molecule,
+            path,
+            nurange=[-np.inf, np.inf],
+            crit=-np.inf,
+            local_databases="~/exojax",
+        ):
+            """
+            Parameters
+            ----------
+            molecule: molecule name
+            path : local path, mirror of ExoMol path
+            nurange : TYPE, optional
+                DESCRIPTION. The default is [-np.inf, np.inf].
+            crit : TYPE, optional
+                DESCRIPTION. The default is -np.inf.
+
+            Returns
+            -------
+            DataFrame
+
+            Examples
+            --------
+            ::
+
+                mdbCH4 = mdbExoMol("CH4", '.database/CH4/12C-1H4/YT10to10/', nus, crit=1.e-30)
+                print(len(mdbCH4.nu_lines), "lines")
+                mdbCH4.elower
+
+            Available columns::
+
+                [
+                    "Sij0",
+                    "logsij0",
+                    "nu_lines",
+                    "A",
+                    "elower",
+                    "eupper",
+                    "gupper",
+                    "jlower",
+                    "jupper",
+                ]
+
+            """
+
+            wavenum_min, wavenum_max = np.min(nurange), np.max(nurange)
+            if wavenum_min == -np.inf:
+                wavenum_min = None
+            if wavenum_max == np.inf:
+                wavenum_max = None
+
+            # Set-up database, download files and set-up cache files if needed
+            mdb = MdbExomol(
+                path,
+                molecule=molecule,
+                local_databases=local_databases,
+                nurange=[wavenum_min, wavenum_max],
+            )
+
+            # Get cache files to load :
+            mgr = mdb.get_datafile_manager()
+            local_files = [mgr.cache_file(f) for f in mdb.trans_file]
+
+            # Load them:
+            jdict = mdb.load(
+                local_files,
+                columns=[k for k in self.__slots__ if k not in ["logsij0"]],
+                lower_bound=([("nu_lines", wavenum_min)] if wavenum_min else [])
+                + ([("Sij0", mdb.crit)] if not np.isneginf(mdb.crit) else []),
+                upper_bound=([("nu_lines", wavenum_max)] if wavenum_max else []),
+                output="jax",
+            )
+
+            # set attributes, accessible as e.g:  mdb.nu_lines
+            for k in jdict.keys():
+                setattr(self, k, jdict[k])
+
+    nus = np.linspace(1e7 / 1.630e4, 1e7 / 1.6305e4)
+
+    # Download new ExoMol repo (in ~/exomol)
+    mdbCH4 = mdbExoMol(
+        "CH4",
+        ".database/CH4/12C-1H4/YT10to10/",
+        nus,
+        crit=1.0e-30,
+        local_databases=".",  # use local folder
+    )
 
-    print("Testing factory:", pytest.main(["../test/io/test_exomol.py"]))
+    print(len(mdbCH4.nu_lines), "lines")
+    mdbCH4.elower
+
+    # Or use RADIS's folder  (# by default ~/.radisdb/exomol)
+    import radis
+
+    mdbCH4_2 = mdbExoMol(
+        "CH4",
+        "CH4/12C-1H4/YT10to10/",
+        nus,
+        crit=1.0e-30,
+        local_databases=pathlib.Path(radis.config["DEFAULT_DOWNLOAD_PATH"]) / "exomol",
+    )
+    # ... ready to run Jax calculations
+
+    #%%
+
+    # """ExoMol lines can be downloaded and accessed separately using
+    # :py:func:`~radis.io.exomol.fetch_exomol`
+    # """
+
+    # # See line data:
+    # from radis.io.exomol import fetch_exomol
+
+    # df = fetch_exomol("SiO", database="EBJT", isotope="1", load_wavenum_max=5000)
+    # print(df)
+
+    # #%%
+    # # See the list of recommended databases for the 1st isotope of SiO :
+    # from radis.io.exomol import get_exomol_database_list, get_exomol_full_isotope_name
+
+    # databases, recommended = get_exomol_database_list(
+    #     "SiO", get_exomol_full_isotope_name("SiO", 1)
+    # )
+    # print("Databases for SiO: ", databases)
+    # print("Database recommended by ExoMol: ", recommended)
diff --git a/radis/io/exomolapi.py b/radis/io/exomolapi.py
index 3c5fa5235..5456324b4 100644
--- a/radis/io/exomolapi.py
+++ b/radis/io/exomolapi.py
@@ -4,9 +4,19 @@
 code (which you should also have a look at !), by @HajimeKawahara, under MIT License.
 
 """
+import bz2
+import re
+
 import numpy as np
 import pandas as pd
 
+try:
+    from .hdf5 import vaexsafe_colname
+except:
+    from radis.io.hdf5 import vaexsafe_colname
+
+from radis.misc.warning import InconsistentDatabaseError
+
 
 def read_def(deff):
     """Exomol IO for a definition file
@@ -17,16 +27,22 @@ def read_def(deff):
 
     Returns
     -------
-    temperature exponent n_Texp
-    broadening parameter alpha_ref
-    molecular mass
-    numinf: nu minimum for trans
-    numtag: tag for wavelength range
+    n_Texp : float
+        temperature exponent
+    alpha_ref : float
+        broadening parameter
+    molmass : float
+        molecular mass
+    numinf : List[float]
+        limit points (``[w(0), w(1), ..., w(n)]``, n+1 elements) defining
+        the spectral ranges appearing in the name of *.trans.bz2 files
+        (``["w(0)-w(1)", "w(1)-w(2)", ..., w(n-1)-w(n)]``, n elements)
+    numtag : List[str]
+        tag for wavelength ranges.
 
     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"))
@@ -50,10 +66,9 @@ def read_def(deff):
             molmass = float(dat["VAL"][i].split()[0])  # in Da (atomic unit)
 
         elif "Lifetime availability" in com:
-            lifetime = dat["VAL"][i] == 1
+            lifetime = int(dat["VAL"][i]) == 1
         elif "Lande g-factor availability" in com:
-            lande = dat["VAL"][i] == 1
-
+            lande = int(dat["VAL"][i]) == 1
         elif "Quantum label" in com:
             quantum_labels.append(dat["VAL"][i].strip(" "))
 
@@ -177,12 +192,17 @@ def read_trans(transf, engine="vaex"):
                 convert=False,  #  file is created by MdbMol
             )
         except:
-            dat = vaex.read_csv(
-                transf,
-                sep=r"\s+",
-                names=("i_upper", "i_lower", "A", "nu_lines"),
-                convert=False,  #  file is created by MdbMol
-            )
+            try:
+                dat = vaex.read_csv(
+                    transf,
+                    sep=r"\s+",
+                    names=("i_upper", "i_lower", "A", "nu_lines"),
+                    convert=False,  #  file is created by MdbMol
+                )
+            except Exception as err:
+                raise Exception(
+                    f"Error reading {transf}. Maybe the file was corrupted during download ? You can try to delete it."
+                ) from err
     elif engine == "csv":
         try:  # bz2 compression
             dat = pd.read_csv(
@@ -192,14 +212,21 @@ def read_trans(transf, engine="vaex"):
                 names=("i_upper", "i_lower", "A", "nu_lines"),
             )
         except:
-            dat = pd.read_csv(
-                transf, sep=r"\s+", names=("i_upper", "i_lower", "A", "nu_lines")
-            )
+            try:
+                dat = pd.read_csv(
+                    transf, sep=r"\s+", names=("i_upper", "i_lower", "A", "nu_lines")
+                )
+            except Exception as err:
+                raise Exception(
+                    f"Error reading {transf}. Maybe the file was corrupted during download ? You can try to delete it."
+                ) from err
+    else:
+        raise NotImplementedError(engine)
 
     return dat
 
 
-def read_states(statesf, dic_def, engine="vaex"):
+def read_states(statesf, dic_def, engine="vaex", skip_optional_data=True):
     """Exomol IO for a state file
 
     Notes
@@ -216,15 +243,42 @@ def read_states(statesf, dic_def, engine="vaex"):
 
     Parameters
     ----------
-    statesf: state file
-    dic_def: Info from def file to read extra quantum numbers
-    engine: parsing engine to use ('vaex', 'csv')
+    statesf: str
+        state file
+    dic_def: dict
+        Info from def file to read extra quantum numbers
+    engine: str
+        parsing engine to use ('vaex', 'csv')
+    skip_optional_data: bool
+        If False, fetch all fields which are marked as available in the ExoMol definition
+        file. If True, load only the first 4 columns of the states file
+        ("i", "E", "g", "J"). The structure of the columns above 5 depend on the
+        the definitions file (*.def) and the Exomol version.
+        If ``skip_optional_data=False``, two errors may occur:
+
+            - a field is marked as present/absent in the *.def field but is
+              absent/present in the *.states file (ie both files are inconsistent).
+            - in the updated version of Exomol, new fields have been added in the
+              states file of some species. But it has not been done for all species,
+              so both structures exist. For instance, the states file of
+              https://exomol.com/data/molecules/HCl/1H-35Cl/HITRAN-HCl/ follows the
+              structure described in [1]_, unlike the states file of
+              https://exomol.com/data/molecules/NO/14N-16O/XABC/ which follows the
+              structure described in [2]_.
 
     Returns
     -------
     states data in pandas DataFrame
 
     If ``'vaex'``, also writes a local hdf5 file `statesf.with_suffix('.hdf5')`
+
+
+    References
+    ----------
+
+    .. [1] Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., Barton, E. J., Chubb, K. L., Coles, P. A., … Zak, E. (2016). The ExoMol database: molecular line lists for exoplanet and other hot atmospheres. https://doi.org/10.1016/j.jms.2016.05.002
+    .. [2] Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., Clark, V. H. J., Chubb, K. L., Conway, E. K., … Yurchenko, O. P. (2020). The 2020 release of the ExoMol database: Molecular line lists for exoplanet and other hot atmospheres. Journal of Quantitative Spectroscopy and Radiative Transfer, 255, 107228. https://doi.org/10.1016/j.jqsrt.2020.107228
+
     """
     # we read first 4 columns for ("i", "E", "g", "J"),
     # skip lifetime, skip Landé g-factor,
@@ -232,24 +286,45 @@ def read_states(statesf, dic_def, engine="vaex"):
     N = np.size(dic_def["quantum_labels"])
     quantum_labels = dic_def["quantum_labels"]
 
-    usecol = np.arange(4)
-    names = ("i", "E", "g", "J")
-    if dic_def["Landé"] and dic_def["lifetime"]:
-        usecol = np.concatenate((usecol, 5 + 2 + np.arange(N)))
-        names = names + ("Lande", "lifetime") + tuple(quantum_labels)
-    elif dic_def["Landé"]:
-        usecol = np.concatenate((usecol, 5 + 1 + np.arange(N)))
-        names = names + ("Lande") + tuple(quantum_labels)
-    elif dic_def["lifetime"]:
-        usecol = np.concatenate((usecol, 5 + 1 + np.arange(N)))
-        names = names + ("lifetime") + tuple(quantum_labels)
-    else:  # no lifetime, nor landé according to the def file
-        usecol = np.concatenate((usecol, 5 + np.arange(N)))
-        names = names + tuple(quantum_labels)
+    mandatory_usecol = np.arange(4)
+    mandatory_fields = ("i", "E", "g", "J")
+    if skip_optional_data:
+        usecol = mandatory_usecol
+        names = mandatory_fields
+    else:
+        if dic_def["Landé"] and dic_def["lifetime"]:
+            usecol = np.concatenate((mandatory_usecol, 5 + 2 + np.arange(N)))
+            names = mandatory_fields + ("Lande", "lifetime") + tuple(quantum_labels)
+        elif dic_def["Landé"]:
+            usecol = np.concatenate((mandatory_usecol, 5 + 1 + np.arange(N)))
+            names = mandatory_fields + ("Lande",) + tuple(quantum_labels)
+        elif dic_def["lifetime"]:
+            usecol = np.concatenate((mandatory_usecol, 5 + 1 + np.arange(N)))
+            names = mandatory_fields + ("lifetime",) + tuple(quantum_labels)
+        else:  # no lifetime, nor landé according to the def file
+            usecol = np.concatenate((mandatory_usecol, 5 + np.arange(N)))
+            names = mandatory_fields + tuple(quantum_labels)
+        # The definitions file (*.def) specifies which fields are available
+        # in the states fiel (*.states). Check the number of available fields
+        # (see *.def) match the numbers of columns in the states file (*.states).
+        # Otherwise, both files are inconsistent.
+        with bz2.open(statesf, "rb") as f:
+            firstline = f.readline().decode("utf-8")  # read 1dt line
+            splitline = [x for x in re.split(r"\s+", firstline) if x != ""]
+            # splitline = re.split(r"\s+", firstline)
+            f.close()
+        if len(splitline) != len(names):
+            raise InconsistentDatabaseError(
+                "The EXOMOL definitions and states files are inconsistent.\n"
+                + "Some data are specified as available by the *.def file, but are absent in the *.bz2 file.\n"
+                + "Set `skip_optional_data=False` in `fetch_exomol()` to load only the required data for a LTE computation.\n"
+                + "The problematic optional data/columns will be ignored."
+            )
 
     if engine == "vaex":
         import vaex
 
+        # TODO Refactor: move in DataFileManager
         try:
             dat = vaex.from_csv(
                 statesf,
@@ -257,12 +332,22 @@ def read_states(statesf, dic_def, engine="vaex"):
                 sep=r"\s+",
                 usecols=usecol,
                 names=names,
-                convert=True,
+                convert=False,  # written in MolDB
             )
         except:
-            dat = vaex.read_csv(
-                statesf, sep=r"\s+", usecols=usecol, names=names, convert=True
-            )
+            try:
+                dat = vaex.read_csv(
+                    statesf,
+                    sep=r"\s+",
+                    usecols=usecol,
+                    names=names,
+                    convert=False,  # written in MolDB
+                )
+
+            except Exception as err:
+                raise Exception(
+                    f"Error reading {statesf}. Maybe the file was corrupted during download ? You can try to delete it."
+                ) from err
     elif engine == "csv":
         try:
             dat = pd.read_csv(
@@ -276,17 +361,20 @@ def read_states(statesf, dic_def, engine="vaex"):
     return dat
 
 
-def pickup_gE(ndstates, ndtrans, trans_file, dic_def, trans_lines=False):
+def pickup_gE(
+    states, trans, trans_file, dic_def, trans_lines=False, skip_optional_data=True
+):
     """extract g_upper (gup), E_lower (elower), and J_lower and J_upper from states
     DataFrame and insert them into the transition DataFrame.
 
     Parameters
     ----------
-    ndstates: states numpy array    - the i, E, g, J are in the 4 first columns
-    ndtrans: transition numpy array
+    states: states DataFrame  - the i, E, g, J are in the 4 first columns
+    trans: transition numpy array
     trans_file: name of the transition file
     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
+    add_quantum_labels: bool . If True fetch all quantum labels in dic_def['quantum_labels'] from states into transitions
 
     Returns
     -------
@@ -296,89 +384,69 @@ def pickup_gE(ndstates, ndtrans, trans_file, dic_def, trans_lines=False):
     -----
     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.
 
+    States by default has columns ::
+
+        #       i      E             g    J    v
+        0       1      0.0           1    0    0
+        1       2      1.448467      3    1    0
+        2       3      4.345384      5    2    0
+        3       4      8.690712      7    3    0
+
 
     """
     ### Step 1. Essential quantum number for spectra
+    # ----------------------------------------------
+
+    def map_add(col, new_col, trans_key, states_key="i"):
+        """Lookup `key` in states and add it in trans, using the level ``i``
+        for upper and lower state
+
+        Examples
+        --------
+        ::
+
+            map_add("E", "E_lower", "i_lower")
+        """
+        try:  # pytable
+            trans[new_col] = trans[trans_key].map(dict(states[col]))
+        except:  # a priori, vaex version  (TODO : replace with dict() approach in vaex too)
+            col = vaexsafe_colname(col)
+            new_col = vaexsafe_colname(new_col)
+
+            # WIP. First implementation with join(). Use Map() will probably be faster.
+            trans.join(
+                states[states_key, col],
+                left_on=trans_key,
+                right_on=states_key,
+                inplace=True,
+            )
+            trans.drop(states_key, inplace=True)
+            trans.rename(col, new_col)
 
-    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:]
+    map_add("g", "gup", "i_upper")
+    map_add("J", "jlower", "i_lower")
+    map_add("J", "jupper", "i_upper")
 
-    i_upper = np.array(ndtrans[:, 0], dtype=int)
-    i_lower = np.array(ndtrans[:, 1], dtype=int)
+    map_add("E", "elower", "i_lower")
 
-    ### 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 = ndstates.to_numpy()[
-    ndstates_extra = ndstates[:, 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]
-
-    # 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)
-    del i_lower
-    jupper = np.array(newstates[i_upper, 2], dtype=int)
-    del i_upper
-
-    A = ndtrans[:, 2]
-
-    if trans_lines:
-        nu_lines = ndtrans[:, 3]
-    else:
-        nu_lines = eupper - elower
-    del eupper
-
-    # ### MASKING ###
-    mask = nu_lines > 0.0
-    if False in mask:
-        len_org = len(nu_lines)
-
-        A = A[mask]
-        nu_lines = nu_lines[mask]
-        elower = elower[mask]
-        gup = gup[mask]
-        jlower = jlower[mask]
-        jupper = jupper[mask]
-        for key, array in dico_quantumNumbers.items():
-            dico_quantumNumbers[key] = array[mask]
-        print(
-            "WARNING: {0:,} transitions with the wavenumber=zero in {1} have been ignored.".format(
-                len_org - len(nu_lines), trans_file
-            )
-        )
-        if trans_lines:
-            print(
-                "This is because the value for the wavenumber column in the transition file is zero for those transitions."
-            )
-        else:
-            print(
-                "This is because the upper and lower state IDs in the transition file indicate the same energy level when referring to the states file for those transitions."
-            )
+    def has_nan(column):
+        try:  # Vaex
+            return column.countnan() > 0
+        except AttributeError:  # Pandas
+            return column.hasnans
 
-    # See Issue exojax#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()
+    if not "nu_lines" in trans or has_nan(trans["nu_lines"]):
+        map_add("E", "eupper", "i_upper")
+        trans["nu_lines"] = trans["eupper"] - trans["elower"]
+
+    ### Step 2. Extra quantum numbers (e/f parity, vib and rot numbers)
+    # -----------------------------------------------------------------
+    if not skip_optional_data:
+        for q in dic_def["quantum_labels"]:
+            map_add(q, f"{q}_l", "i_lower")
+            # map_add(q, f"{q}_u", "i_upper")
 
-    return A, nu_lines, elower, gup, jlower, jupper, mask, dico_quantumNumbers
+    return trans
 
 
 # def pickup_gEslow(states, trans):
@@ -559,69 +627,3 @@ def make_jj2b(bdat, j2alpha_ref_def, j2n_Texp_def, jupper_max=None):
 
 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/geisa.py b/radis/io/geisa.py
index 35b8120da..eaa33cb56 100644
--- a/radis/io/geisa.py
+++ b/radis/io/geisa.py
@@ -19,10 +19,17 @@
 import numpy as np
 
 import radis
-from radis.io.cache_files import cache_file_name, load_h5_cache_file, save_to_hdf
-from radis.io.dbmanager import DatabaseManager
-from radis.io.tools import drop_object_format_columns, parse_hitran_file
 
+try:
+    from .cache_files import load_h5_cache_file, save_to_hdf
+    from .dbmanager import DatabaseManager
+    from .hdf5 import DataFileManager
+    from .tools import drop_object_format_columns, parse_hitran_file
+except ImportError:
+    from radis.io.cache_files import load_h5_cache_file, save_to_hdf
+    from radis.io.dbmanager import DatabaseManager
+    from radis.io.hdf5 import DataFileManager
+    from radis.io.tools import drop_object_format_columns, parse_hitran_file
 # from typing import Union
 
 
@@ -229,7 +236,7 @@ def gei2df(
         print("Last Modification time: {0}".format(metadata["last_modification"]))
 
     # Attempt to use cache file
-    fcache = cache_file_name(fname, engine=engine)
+    fcache = DataFileManager(engine).cache_file(fname)
     if cache and exists(fcache):
         relevant_if_metadata_above = (
             {"wavenum_max": load_wavenum_min} if load_wavenum_min else {}
@@ -388,7 +395,7 @@ def parse_to_local_file(
         # molecule = self.molecule
         # Ntotal_lines_expected = GEISA_MOLECULES_Nlines[molecule.upper()]
 
-        writer = self.get_hdf5_manager()
+        writer = self.get_datafile_manager()
 
         with opener.open(urlname) as gfile:  # locally downloaded file
 
diff --git a/radis/io/hdf5.py b/radis/io/hdf5.py
index 7824d0129..61b5972b7 100644
--- a/radis/io/hdf5.py
+++ b/radis/io/hdf5.py
@@ -6,6 +6,7 @@
 """
 
 import os
+import pathlib
 import sys
 from os.path import exists, expanduser, splitext
 from time import time
@@ -15,6 +16,14 @@
 from tables.exceptions import NoSuchNodeError
 
 
+def vaexsafe_colname(name):
+    """replace '/' (forbidden in HDF5 vaex column names with '_'
+    https://github.com/radis/radis/issues/473
+    https://github.com/vaexio/vaex/issues/1255
+    """
+    return name.replace("/", "_")
+
+
 def update_pytables_to_vaex(fname, remove_initial=False, verbose=True, key="df"):
     """Convert a HDF5 file generated from PyTables to a
     Vaex-friendly HDF5 format, preserving metadata"""
@@ -30,7 +39,7 @@ def update_pytables_to_vaex(fname, remove_initial=False, verbose=True, key="df")
     df = pd.read_hdf(fname)
     df = vaex.from_pandas(df)
 
-    pytables_manager = HDF5Manager(engine="pytables")
+    pytables_manager = DataFileManager(engine="pytables")
 
     # Read metadata
     file_metadata = pytables_manager.read_metadata(fname)
@@ -40,7 +49,7 @@ def update_pytables_to_vaex(fname, remove_initial=False, verbose=True, key="df")
     df.close()  # try to fix file not closed()  TODO: remove?
     del df  # same TODO
 
-    vaex_manager = HDF5Manager(engine="vaex")
+    vaex_manager = DataFileManager(engine="vaex")
     vaex_manager.add_metadata(fname_vaex, file_metadata)
 
     if verbose:
@@ -56,8 +65,13 @@ def update_pytables_to_vaex(fname, remove_initial=False, verbose=True, key="df")
 
 
 class HDF5Manager(object):
+    def __init__(*args, **kwargs):
+        raise DeprecationWarning("HDF5Manager replaced with DataFileManager")
+
+
+class DataFileManager(object):
     def __init__(self, engine=None):
-        """Class to handle all memory-mapping-librairies with one common API
+        """Class to handle all DataFrame-librairies with one common API
 
         All functions may not be fully implemetned, will raise a NotImplementedError
         if that's not the case.
@@ -67,13 +81,14 @@ def __init__(self, engine=None):
             'vaex'     > HDF5,  column-based
             'pytables' > Pandas's HDF5,  row-based
             'h5py'     > HDF5
+            'feather'  > feather
 
         Functions ::
 
             add_metadata
             read_metadata
             write
-            load
+            read
             guess_engine
 
         Examples
@@ -107,13 +122,17 @@ def write(
         self,
         file,
         df,
-        append=True,
+        append=False,
         key="default",
         format="table",
-        data_columns=["iso", "wav"],
+        data_columns=["iso", "wav", "nu_lines"],
     ):
         """Write dataframe ``df`` to ``file``
 
+        Parameters
+        ----------
+        df: DataFrame
+
         Other Parameters
         ----------------
         key: str
@@ -141,9 +160,16 @@ def write(
             # export dataframe
             df.to_hdf(file, key, format="fixed", mode="w", complevel=9, complib="blosc")
         elif self.engine == "vaex":
+            if isinstance(df, pd.DataFrame):
+                import vaex
+
+                df = vaex.from_pandas(df)
+
+            for c in df.columns:  # remove "/" in columns (forbidden)
+                df.rename(c, vaexsafe_colname(c))
+
             if key == "default":
                 key = r"/table"
-            import vaex
 
             if append == True:
                 # In vaex we cannot append. Here we write lots of small files then combine them.
@@ -163,10 +189,9 @@ def write(
                     os.remove(file)
                 self._temp_batch_files.append(file)
             # Write:
-            try:
-                df.export_hdf5(file, group=key, mode="w")
-            except AttributeError:  # case where df is not a Vaex dataFrame but (likely) a Pandas Dataframe
-                vaex.from_pandas(df).export_hdf5(file, group=key, mode="w")
+            df.export_hdf5(file, group=key, mode="w")
+        elif self.engine == "feather":
+            df.to_feather(file)
         else:
             raise NotImplementedError(self.engine)
             # h5py is not designed to write Pandas DataFrames
@@ -213,6 +238,104 @@ def __del__(self):
             self._close_temp_batch_files()
 
     def load(
+        self,
+        fname,
+        columns=None,
+        lower_bound=[],
+        upper_bound=[],
+        within=[],
+        output="pandas",
+        **store_kwargs,
+    ):
+        """
+        Other Parameters
+        ----------------
+        columns: list of str
+            list of columns to load. If ``None``, returns all columns in the file.
+        output: 'pandas', 'vaex', 'jax'
+            format of the output DataFrame. If ``'jax'``, returns a dictionary of
+            jax arrays.
+        lower_bound: list of tuples [(column, lower_bound), etc.]
+            ::
+
+                lower_bound =[("wav", load_wavenum_min)]
+        upper_bound_bound: list of tuples [(column, upper_bound), etc.]
+            ::
+
+                upper_bound=[("wav", load_wavenum_max)]
+        within: list of tuples [(column, within_list), etc.]
+            ::
+
+                within=[("iso", isotope.split(","))]
+        """
+        if not lower_bound and not upper_bound and not within:
+            df = self.read(
+                fname,
+                columns,
+                **store_kwargs,
+            )
+        else:
+            df = self.read_filter(
+                fname,
+                columns,
+                lower_bound=lower_bound,
+                upper_bound=upper_bound,
+                within=within,
+                **store_kwargs,
+            )
+
+        # Load
+        # Convert to output format
+        engine = self.engine
+
+        if engine == output:
+            df = df
+        elif engine == "vaex":
+            # in vaex, column selection has to happen now
+            if columns:  # load only these columns (if they exist)
+                columns = [c for c in columns if c in df.columns]
+            if output == "pandas":
+                df = df.to_pandas_df(column_names=columns)
+            elif output == "jax":
+                if columns == None:
+                    columns = list(df.columns)
+                out = {}
+                try:
+                    import jax.numpy as jnp
+                except ImportError:
+                    print("Jax not found. Using Numpy.")
+                    import numpy as jnp
+                for c in columns:
+                    if (
+                        c == "Sij0"
+                    ):  # (special case for Sij0 : we store logSij0 instead, to store them as jax np.32 arrays)
+                        import numpy as np
+
+                        out["logsij0"] = jnp.array(np.log(df[c].values))
+                    else:
+                        out[c] = jnp.array(df[c].values)
+                df = out
+            else:
+                raise NotImplementedError(f"output {output} for engine {engine}")
+        elif engine == "pytables":
+            if output == "pandas":
+                df = df
+            else:
+                raise NotImplementedError(f"output {output} for engine {engine}")
+        elif engine == "feather":
+            if output == "pandas":
+                df = df
+            else:
+                raise NotImplementedError(f"output {output} for engine {engine}")
+        else:
+            raise NotImplementedError(output)
+
+        if output == "vaex":
+            df = df.extract()  # return DataFrame containing only the filtered rows
+
+        return df
+
+    def read(
         self,
         fname,
         columns=None,
@@ -224,6 +347,7 @@ def load(
         """
         Parameters
         ----------
+        fname: str
         columns: list of str
             list of columns to load. If ``None``, returns all columns in the file.
         where: list of str
@@ -290,7 +414,7 @@ def load(
                 except (FileNotFoundError, OSError) as err:
                     # error message with suggestion on how to convert from existing file
                     for f in fname_list:
-                        if exists(f.replace(".hdf5", ".h5")):
+                        if f.endswith(".hdf5") and exists(f.replace(".hdf5", ".h5")):
                             raise FileNotFoundError(
                                 f"`{f}` not found but `{f.replace('.hdf5', '.h5')}` exists (probably a row-based pytables HDF5 file). Try (1) using engine='pytables' in the calling function (`hdf2df`, `fetch_hitemp`, etc.)  ; (2) delete the file to re-download and re-parse it (this may take a lot of time !) ;  or (3, recommended) set `import radis; radis.config['AUTO_UPDATE_DATABASE']= True` in your script to auto-update to Vaex HDF5 file"
                             ) from err
@@ -322,11 +446,107 @@ def load(
                     out[k] = f[k][()]
             return pd.DataFrame(out)
 
+        elif self.engine == "feather":
+            assert where is None
+            fname = expanduser(fname)
+            return pd.read_feather(fname)
+
+        else:
+            raise NotImplementedError(self.engine)
+
+        return df
+
+    def read_filter(
+        self,
+        fname,
+        columns=None,
+        lower_bound=[],
+        upper_bound=[],
+        within=[],
+        **store_kwargs,
+    ):
+        """
+        Parameters
+        ----------
+        fname: str
+        columns: list of str
+            list of columns to load. If ``None``, returns all columns in the file.
+        lower_bound: list of tuples [(column, lower_bound), etc.]
+            ::
+
+                lower_bound =[("wav", load_wavenum_min)]
+        upper_bound_bound: list of tuples [(column, upper_bound), etc.]
+            ::
+
+                upper_bound=[("wav", load_wavenum_max)]
+        within: list of tuples [(column, within_list), etc.]
+            ::
+
+                within=[("iso", isotope.split(","))]
+
+        """
+
+        # Selection
+        if self.engine == "pytables":
+            # Selection
+            where = []
+            for (column, lbound) in lower_bound:
+                where.append(f"{column} > {lbound}")
+            for (column, ubound) in upper_bound:
+                where.append(f"{column} < {ubound}")
+            for (column, withinv) in within:
+                where.append(f"{column} in {withinv.split(',')}")
+
+        elif self.engine in ["vaex", "feather"]:
+            # Selection is done after opening the file time in vaex
+            # see end of this function
+            where = None
         else:
             raise NotImplementedError(self.engine)
 
+        # Load :
+        df = self.read(fname, columns=columns, where=where, **store_kwargs)
+
+        #  Selection in vaex
+        if self.engine in ["vaex", "feather"]:
+            # (note that in Vaex, the selection happens on disk whereas Feather
+            # is already loaded as a Pandas DataFrame in RAM)
+
+            # Selection
+            b = True
+            for (column, lbound) in lower_bound:
+                b *= df[column] > lbound
+            for (column, ubound) in upper_bound:
+                b *= df[column] < ubound
+            for (column, withinv) in within:
+                b2 = False
+                for val in withinv.split(","):
+                    b2 += df[column] == float(val)
+                b *= b2
+            if b is not True and False in b:
+                df = df[
+                    b
+                ]  # note in Vaex mode, this is a vaex Expression, not the DataFrame yet
+
         return df
 
+    def cache_file(self, fname):
+        """Return the corresponding cache file name for fname.
+
+        Other Parameters
+        ----------------
+        engine: ``'h5py'``, ``'pytables'``, ``'vaex'``
+           which HDF5 library to use. Default ``pytables``
+        """
+        if self.engine in ["pytables", "pytables-fixed"]:
+            return pathlib.Path(fname).with_suffix(".h5")
+        elif self.engine in ["h5py", "vaex"]:
+            return pathlib.Path(fname).with_suffix(".hdf5")
+        elif self.engine == "feather":
+            return pathlib.Path(fname).with_suffix(".feather")
+        else:
+            raise ValueError(self.engine)
+
     def add_metadata(
         self, fname: str, metadata: dict, key="default", create_empty_dataset=False
     ):
@@ -424,6 +644,14 @@ def read_metadata(self, fname: str, key="default") -> dict:
                 except KeyError as err:
                     print(f"Error reading metadata from {fname}")
                     raise err
+                except AttributeError as err:
+                    if "Attribute 'metadata' does not exist" in str(err):
+                        metadata = {}
+                    else:
+                        raise err
+
+        elif self.engine == "feather":
+            return {}  # no metadata
 
         elif self.engine == "h5py":
             fname = expanduser(fname)
@@ -477,6 +705,18 @@ def read_metadata(self, fname: str, key="default") -> dict:
 
         return metadata
 
+    def to_numpy(self, df):
+        """Convert DataFrame to numpy"""
+
+        if self.engine == "vaex":
+            import vaex
+
+            return vaex.array_types.to_numpy(df)
+        elif self.engine == "feather":
+            return df.to_numpy()
+        else:
+            raise NotImplementedError(self.engine)
+
     @classmethod
     def guess_engine(self, file, verbose=True):
         """Guess which HDF5 library ``file`` is compatible with
@@ -496,26 +736,38 @@ def guess_engine(self, file, verbose=True):
             mgr.read_metadata(file)
 
         """
-        # See if it looks like PyTables
-        import tables
-
-        if tables.is_pytables_file(file):
-            engine = "pytables"
+        if file.endswith(".feather"):
+            engine = "feather"
         else:
-            # Try Vaex
-            file = expanduser(file)
-            with h5py.File(file, mode="r") as hf:
-                try:
-                    hf[r"/table"]
-                except KeyError:
-                    engine = "h5py"
-                else:
-                    engine = "vaex"
+            # See if it looks like PyTables
+            import tables
+
+            if tables.is_pytables_file(file):
+                engine = "pytables"
+            else:
+                # Try Vaex
+                file = expanduser(file)
+                with h5py.File(file, mode="r") as hf:
+                    try:
+                        hf[r"/table"]
+                    except KeyError:
+                        engine = "h5py"
+                    else:
+                        engine = "vaex"
         if verbose:
-            print(f"Guessed that {file} was compatible with `{engine}` hdf5 engine")
+            print(f"Guessed that {file} was compatible with `{engine}` engine")
         # raise
         return engine
 
+    def has_nan(self, column):
+        if self.engine == "vaex":
+            b = column.isnan()  # TODO: check if can be made faster?
+            return b.sum() > 0
+        elif self.engine in ["pytables", "feather"]:
+            return column.hasnans
+        else:
+            raise NotImplementedError(self.engine)
+
 
 def hdf2df(
     fname,
@@ -557,7 +809,7 @@ def hdf2df(
 
     Returns
     -------
-    df: pandas Dataframe
+    df: pandas Dataframe, or vaex DataFrameLocal, or dictionary of Jax arrays
         dataframe containing all lines or energy levels
 
     Examples
@@ -588,88 +840,26 @@ def hdf2df(
     """
 
     if engine == "guess":
-        engine = HDF5Manager.guess_engine(fname)
+        engine = DataFileManager.guess_engine(fname)
 
     t0 = time()
 
-    # Selection
-    selection = False
-    if engine == "pytables":
-        # Selection
-        selection = True
-        where = []
-        if load_wavenum_min is not None:
-            where.append(f"wav > {load_wavenum_min}")
-        if load_wavenum_max is not None:
-            where.append(f"wav < {load_wavenum_max}")
-        if isotope:
-            where.append(f'iso in {isotope.split(",")}')
-
-    elif engine == "vaex":
-        # Selection is done after opening the file time in vaex
-        # see end of this function
-        where = None
-    else:
-        raise NotImplementedError(engine)
-
-    # Load :
-    manager = HDF5Manager(engine)
-    df = manager.load(fname, columns=columns, where=where, **store_kwargs)
-
-    #  Selection in vaex
-    if engine == "vaex":
-
-        # Selection
-        selection = True
-        b = True
-        if load_wavenum_min is not None:
-            b *= df.wav > load_wavenum_min
-        if load_wavenum_max is not None:
-            b *= df.wav < load_wavenum_max
-        if isotope is not None:
-            from radis.misc.basics import is_float
-
-            if is_float(isotope):
-                b *= df.iso == int(isotope)
-            else:
-                b2 = False
-                for iso in isotope.split(","):
-                    b2 += df.iso == int(iso)
-                b *= b2
-        if b != True and False in b:
-            df = df[b]  # note that this is a vaex Expression, not the DataFrame yet
+    manager = DataFileManager(engine)
 
-        # Load
-        if columns:  # load only these columns (if they exist)
-            columns = [c for c in columns if c in df.columns]
-
-    # Convert to output format
-    if output == "pandas" and engine == "vaex":
-        df = df.to_pandas_df(column_names=columns)
-    elif output == "vaex" and engine == "pytables":
-        raise NotImplementedError
-    elif output == "jax":
-        if columns == None:
-            columns = ["wav", "int", "A", "El", "gpp"]
-        out = {}
-        try:
-            import jax as jnp
-        except ImportError:
-            import numpy as jnp
-        for c in columns:
-            out[c] = jnp.array(df[c].values)
-        df = out
-    elif output == engine or (
-        engine == "pytables" and output == "pandas"
-    ):  # pandas > pandas; vaex -> vaex
-        df = df
-    else:
-        raise NotImplementedError(output)
+    df = manager.load(
+        fname,
+        columns,
+        lower_bound=[("wav", load_wavenum_min)] if load_wavenum_min is not None else [],
+        upper_bound=[("wav", load_wavenum_max)] if load_wavenum_max is not None else [],
+        within=[("iso", isotope)] if isotope is not None else [],
+        output=output,
+    )
 
     # Read and add metadata in the DataFrame
     metadata = manager.read_metadata(fname)
 
     # Sanity Checks if loading the full file
+    selection = isotope or load_wavenum_min or load_wavenum_max
     if not selection:
         if "total_lines" in metadata:
             assert len(df) == metadata["total_lines"]
@@ -679,7 +869,7 @@ def hdf2df(
             assert df["wav"].max() == metadata["wavenumber_max"]
 
     # Update metadata
-    if output != "jax":
+    if output != "jax":  # metadata not given in Jax output
         if isinstance(metadata, list):
             metadata_dict = {}
             for k, v in metadata[0].items():
diff --git a/radis/io/hitemp.py b/radis/io/hitemp.py
index c6dddc379..323e413ec 100644
--- a/radis/io/hitemp.py
+++ b/radis/io/hitemp.py
@@ -312,7 +312,7 @@ def parse_to_local_file(
         wmin = np.inf
         wmax = 0
 
-        writer = self.get_hdf5_manager()
+        writer = self.get_datafile_manager()
 
         with opener.open(urlname) as gfile:  # locally downloaded file
 
@@ -622,9 +622,10 @@ def fetch_hitemp(
     df = ldb.load(
         files_loaded,  # filter other files,
         columns=columns,
-        isotope=isotope,
-        load_wavenum_min=load_wavenum_min,  # for relevant files, get only the right range
-        load_wavenum_max=load_wavenum_max,
+        within=[("iso", isotope)] if isotope is not None else [],
+        # for relevant files, get only the right range :
+        lower_bound=[("wav", load_wavenum_min)] if load_wavenum_min is not None else [],
+        upper_bound=[("wav", load_wavenum_max)] if load_wavenum_max is not None else [],
         output=output,
     )
 
diff --git a/radis/io/hitran.py b/radis/io/hitran.py
index 9726770a0..87e22ad29 100644
--- a/radis/io/hitran.py
+++ b/radis/io/hitran.py
@@ -50,13 +50,26 @@
     HITRAN_GROUP6,
     get_molecule,
 )
-from radis.io.cache_files import cache_file_name, load_h5_cache_file, save_to_hdf
-from radis.io.dbmanager import DatabaseManager
-from radis.io.tools import (
-    drop_object_format_columns,
-    parse_hitran_file,
-    replace_PQR_with_m101,
-)
+
+try:
+    from .cache_files import load_h5_cache_file, save_to_hdf
+    from .dbmanager import DatabaseManager
+    from .hdf5 import DataFileManager
+    from .tools import (
+        drop_object_format_columns,
+        parse_hitran_file,
+        replace_PQR_with_m101,
+    )
+except ImportError:
+    from radis.io.cache_files import load_h5_cache_file, save_to_hdf
+    from radis.io.dbmanager import DatabaseManager
+    from radis.io.hdf5 import DataFileManager
+    from radis.io.tools import (
+        drop_object_format_columns,
+        parse_hitran_file,
+        replace_PQR_with_m101,
+    )
+
 
 # %% Parsing functions
 
@@ -178,7 +191,7 @@ def hit2df(
     columns = columns_2004
 
     # Use cache file if possible
-    fcache = cache_file_name(fname, engine=engine)
+    fcache = DataFileManager(engine).cache_file(fname)
     if cache and exists(fcache):
         relevant_if_metadata_above = (
             {"wavenum_max": load_wavenum_min} if load_wavenum_min else {}
@@ -1171,7 +1184,7 @@ def download_all_hitran_isotopes(molecule, directory):
             molecule, tempdir
         )
 
-        writer = self.get_hdf5_manager()
+        writer = self.get_datafile_manager()
 
         # Create HDF5 cache file for all isotopes
         Nlines = 0
@@ -1417,9 +1430,10 @@ def fetch_hitran(
     df = ldb.load(
         local_file,
         columns=columns,
-        isotope=isotope,
-        load_wavenum_min=load_wavenum_min,  # for relevant files, get only the right range
-        load_wavenum_max=load_wavenum_max,
+        within=[("iso", isotope)] if isotope is not None else [],
+        # for relevant files, get only the right range :
+        lower_bound=[("wav", load_wavenum_min)] if load_wavenum_min is not None else [],
+        upper_bound=[("wav", load_wavenum_max)] if load_wavenum_max is not None else [],
         output=output,
     )
 
diff --git a/radis/lbl/loader.py b/radis/lbl/loader.py
index e99a64093..69338e37c 100644
--- a/radis/lbl/loader.py
+++ b/radis/lbl/loader.py
@@ -51,7 +51,7 @@
 
 import warnings
 from copy import deepcopy
-from os.path import exists, expanduser, join
+from os.path import exists, expanduser, join, splitext
 from time import time
 from uuid import uuid1
 
@@ -63,7 +63,6 @@
 from radis.db.molecules import getMolecule
 from radis.db.molparam import MolParams
 from radis.db.references import doi
-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.geisa import fetch_geisa
@@ -1289,7 +1288,7 @@ def fetch_databank(
                 )
 
             # Download, setup local databases, and fetch (use existing if possible)
-            if memory_mapping_engine not in ["vaex", "feather"]:
+            if memory_mapping_engine not in ["vaex", "feather", "pytables"]:
                 raise NotImplementedError(
                     f"{memory_mapping_engine} with ExoMol files. Define radis.config['MEMORY_MAPPING_ENGINE'] = 'vaex' or 'feather'"
                 )
@@ -1779,8 +1778,14 @@ def _check_database_params(
 
             filtered_path = [fname for fname in path]
             for fname in path:
-                if cache_file_name(fname) in path and cache_file_name(fname) != fname:
-                    filtered_path.remove(cache_file_name(fname))
+                for likely_fname_cache in [
+                    splitext(fname)[0] + ".h5",
+                    fname + ".h5",
+                    splitext(fname)[0] + ".hdf5",
+                    fname + ".hdf5",
+                ]:
+                    if likely_fname_cache in path and likely_fname_cache != fname:
+                        filtered_path.remove(likely_fname_cache)
             new_paths += filtered_path
 
             # Raise errors if no file / print which files were selected.
diff --git a/radis/misc/basics.py b/radis/misc/basics.py
index 91640fe27..e81054d7a 100644
--- a/radis/misc/basics.py
+++ b/radis/misc/basics.py
@@ -500,9 +500,10 @@ def is_number(s):
 
     try:
         float(s)
-        return True
-    except ValueError:
+    except (ValueError, TypeError):
         return False
+    else:
+        return True
 
 
 def to_str(a):
diff --git a/radis/misc/utils.py b/radis/misc/utils.py
index fd76262a4..7cbf2eeef 100644
--- a/radis/misc/utils.py
+++ b/radis/misc/utils.py
@@ -181,6 +181,8 @@ def get_files_from_regex(path):
     """Returns a list of absolute paths of all the files whose names match the
     input regular expression."""
     directory_name = dirname(path)
+    if directory_name == "":
+        directory_name = "."
     regex = basename(path)
 
     file_names = []
diff --git a/radis/misc/warning.py b/radis/misc/warning.py
index a3818acc4..bce18936c 100644
--- a/radis/misc/warning.py
+++ b/radis/misc/warning.py
@@ -75,6 +75,10 @@ class EmptyDatabaseError(ValueError):
     pass
 
 
+class InconsistentDatabaseError(ValueError):
+    pass
+
+
 # Warnings
 
 


Test patch:
diff --git a/radis/test/files/geisa_O2_fragment.par b/radis/test/files/geisa_O2_fragment.par
index 7a2e45214..e86948a18 100644
--- a/radis/test/files/geisa_O2_fragment.par
+++ b/radis/test/files/geisa_O2_fragment.par
@@ -4,4 +4,4 @@
     0.000003 1.2240D-400.0356  664.2595                       X0                       X0                Q 21Q 21     g0.72 66  7SA3 71 7.246E-30 0.0390 0.000000-.9999 -0.999999-9.9999D-01-.9999-.99 0.000000-.9999-9.9999-.99-.990.0000000.000000-.99-.99
     0.000003 1.2510D-400.0420  664.2605                       X0                       X0                Q 21Q 21     d0.71 66  7G03 71-9.999E-01 0.0412 0.000000-.9999 -0.999999-9.9999D-01-.9999-.99 0.000000-.9999-9.9999-.99-.990.0000000.000000-.99-.99
     0.000003 2.0900D-410.0352  933.5309                       X0                       X0                Q 25Q 25     g0.72 66  7SA3 71 3.862E-30 0.0370 0.000000-.9999 -0.999999-9.9999D-01-.9999-.99 0.000000-.9999-9.9999-.99-.990.0000000.000000-.99-.99
-    0.000003 2.1730D-410.0390  933.5327                       X0                       X0                Q 25Q 25     d0.71 66  7G03 71-9.999E-01 0.0383 0.000000-.9999 -0.999999-9.9999D-01-.9999-.99 0.000000-.9999-9.9999-.99-.990.0000000.000000-.99-.99
\ No newline at end of file
+    0.000003 2.1730D-410.0390  933.5327                       X0                       X0                Q 25Q 25     d0.71 66  7G03 71-9.999E-01 0.0383 0.000000-.9999 -0.999999-9.9999D-01-.9999-.99 0.000000-.9999-9.9999-.99-.990.0000000.000000-.99-.99
diff --git a/radis/test/io/test_exomol.py b/radis/test/io/test_exomol.py
index 9dd6d56c4..1828f72e0 100644
--- a/radis/test/io/test_exomol.py
+++ b/radis/test/io/test_exomol.py
@@ -95,7 +95,7 @@ def test_calc_exomol_vs_hitemp(verbose=True, plot=True, *args, **kwargs):
     }
 
     s_exomol = calc_spectrum(
-        **conditions, databank="exomol", name="EXOMOL/HITEMP (H2 broadened)"
+        **conditions, databank="exomol", name="EXOMOL/HITEMP (default broadening)"
     )
     s_hitemp = calc_spectrum(
         **conditions,
diff --git a/radis/test/io/test_hdf5.py b/radis/test/io/test_hdf5.py
index 031a8b328..87a0713be 100644
--- a/radis/test/io/test_hdf5.py
+++ b/radis/test/io/test_hdf5.py
@@ -7,7 +7,7 @@
 
 import pytest
 
-from radis.io.hdf5 import HDF5Manager, hdf2df
+from radis.io.hdf5 import DataFileManager, hdf2df
 from radis.io.hitemp import fetch_hitemp
 from radis.misc.config import getDatabankEntries
 
@@ -39,25 +39,25 @@ def test_hdf5_io_engines(*args, **kwargs):
             f.create_dataset(c, data=df0[c])
 
     # Test Pytables engine : add_metadata ; load ; read_metadata
-    manager = HDF5Manager(engine="pytables")
+    manager = DataFileManager(engine="pytables")
     manager.add_metadata("test_pytables.h5", metadata0)
 
-    df = manager.load("test_pytables.h5")
+    df = manager.read("test_pytables.h5")
     assert (df == df0).all().all()
     assert manager.read_metadata("test_pytables.h5") == metadata0
 
     # Test h5py engine : add_metadata ; load ; read_metadata
-    manager = HDF5Manager(engine="h5py")
+    manager = DataFileManager(engine="h5py")
     manager.add_metadata("test_h5py.h5", metadata0, key=None)
-    df = manager.load("test_h5py.h5", key=None)
+    df = manager.read("test_h5py.h5", key=None)
     assert (df == df0).all().all()
     assert manager.read_metadata("test_h5py.h5") == metadata0
 
     # Test vaex engine : add_metadata ; load ; read_
     # .. also able to read h5py files
-    manager = HDF5Manager(engine="vaex")
+    manager = DataFileManager(engine="vaex")
     manager.add_metadata("test_h5py.h5", metadata0, key=None)
-    df = manager.load("test_h5py.h5", key=None)
+    df = manager.read("test_h5py.h5", key=None)
     # this time; df is a vaex DataFrame
     # ... it keeps the file open so we should close the file handle
     assert (df.to_pandas_df() == df0).all().all()
diff --git a/radis/test/io/test_hitemp.py b/radis/test/io/test_hitemp.py
index 02c9eff2c..06a1883bf 100644
--- a/radis/test/io/test_hitemp.py
+++ b/radis/test/io/test_hitemp.py
@@ -12,6 +12,29 @@
 from radis.misc.config import getDatabankList
 
 
+@pytest.mark.needs_connection
+@pytest.mark.fast
+def test_no_redownload(*args, **kwargs):
+    """Test that modification time for 2 successive hitemp fetch is the same, i.e.,
+    there is no redownload"""
+
+    # 1st call
+    df, local_files = fetch_hitemp(
+        "OH",
+        return_local_path=True,
+    )
+    from os.path import getmtime
+
+    first_mod_time = getmtime(local_files[0])
+
+    # 2nd call
+    df, local_files = fetch_hitemp(
+        "OH",
+        return_local_path=True,
+    )
+    assert first_mod_time == getmtime(local_files[0])
+
+
 @pytest.mark.fast
 def test_relevant_files_filter():
     from radis.io.hitemp import keep_only_relevant
@@ -425,6 +448,7 @@ def test_parse_hitemp_missing_labels_issue280(*args, **kwargs):
 
 
 if __name__ == "__main__":
+    test_no_redownload()
     test_relevant_files_filter()
     test_fetch_hitemp_OH_pytables()
     test_fetch_hitemp_OH_vaex()