Repository: radis/radis (language: python)

Issue:
Partition function calculation stuck in Spyder with "tabulation" mode (vaex)
### 🐛 Describe the bug
<!--- A clear and concise description of the bug (or unexpected behavior) -->

As already seen in #316 , the just-in-time tabulation of partition functions with `parsum_mode='tabulation'` is stuck if ran within a Spyder iPython console.


### 🖥 Steps To Reproduce
<!--- Steps or code snippet to reproduce the behavior.-->
<!--- Best : use test files so that everyone can run your example : https://radis.readthedocs.io/en/latest/dev/developer.html#test-files -->

Run the following snippet in a Spyder console : 

```python
import radis
radis.DEBUG_MODE = True
from radis import SpectrumFactory
sf = SpectrumFactory(wmin=2283.7,
                      wmax=2285,
                      isotope='1,2,3',
                      verbose=3,    # more for more details
                      parsum_mode="tabulation"     # 👈 comment to run in normal mode
                      )
sf.load_databank('HITEMP-CO2-TEST')        # predefined in ~/radis.json
s = sf.non_eq_spectrum(Tvib=3000, Trot=2000, path_length=1)
s.print_perf_profile()

```

The stuck happens when computing a mean with vaex.

Even more minimal error : 
```python
import numpy as np
import vaex
df = vaex.from_dict({"E":np.random.rand(1000)})
print(df)
print("before EQ mean")
E_bins = df.mean("E", binby="E", shape=10)
print("after EQ mean")       # 👈 I never get to here
```

Script works well if 
- ran from an Ipython ( 7.22.0) console in the terminal 
- ran from Python in the terminal
- ran from VS Code
- ran from Jupyter Lab (3.0.14)
- ran in Jupyter QtConsole



As pointed out by @minouHub , opening a new Spyder IPython console unlocks the process. If console has been "unlocked", I can even re-run the script another time before it gets stuck again. (?)

Last thing : although no graphics are plot anywhere in this example, this problem appears when my iPython console has **Qt5** output graphics. It does not appear with **Inline** graphics (tried both consistently 3 times)



> spyder                    5.0.5            py38haa95532_2
> spyder-kernels            2.0.5            py38haa95532_0
> ipython                   7.22.0           py38hd4e2768_0

on Windows.

Partition function calculation stuck in Spyder with "tabulation" mode (vaex)
### 🐛 Describe the bug
<!--- A clear and concise description of the bug (or unexpected behavior) -->

As already seen in #316 , the just-in-time tabulation of partition functions with `parsum_mode='tabulation'` is stuck if ran within a Spyder iPython console.


### 🖥 Steps To Reproduce
<!--- Steps or code snippet to reproduce the behavior.-->
<!--- Best : use test files so that everyone can run your example : https://radis.readthedocs.io/en/latest/dev/developer.html#test-files -->

Run the following snippet in a Spyder console : 

```python
import radis
radis.DEBUG_MODE = True
from radis import SpectrumFactory
sf = SpectrumFactory(wmin=2283.7,
                      wmax=2285,
                      isotope='1,2,3',
                      verbose=3,    # more for more details
                      parsum_mode="tabulation"     # 👈 comment to run in normal mode
                      )
sf.load_databank('HITEMP-CO2-TEST')        # predefined in ~/radis.json
s = sf.non_eq_spectrum(Tvib=3000, Trot=2000, path_length=1)
s.print_perf_profile()

```

The stuck happens when computing a mean with vaex.

Even more minimal error : 
```python
import numpy as np
import vaex
df = vaex.from_dict({"E":np.random.rand(1000)})
print(df)
print("before EQ mean")
E_bins = df.mean("E", binby="E", shape=10)
print("after EQ mean")       # 👈 I never get to here
```

Script works well if 
- ran from an Ipython ( 7.22.0) console in the terminal 
- ran from Python in the terminal
- ran from VS Code
- ran from Jupyter Lab (3.0.14)
- ran in Jupyter QtConsole



As pointed out by @minouHub , opening a new Spyder IPython console unlocks the process. If console has been "unlocked", I can even re-run the script another time before it gets stuck again. (?)

Last thing : although no graphics are plot anywhere in this example, this problem appears when my iPython console has **Qt5** output graphics. It does not appear with **Inline** graphics (tried both consistently 3 times)



> spyder                    5.0.5            py38haa95532_2
> spyder-kernels            2.0.5            py38haa95532_0
> ipython                   7.22.0           py38hd4e2768_0

on Windows.

Broadening of ExoMol molecules
Hello guys,

### 🎯 Todo
I noticed the broadening of the ExoMol molecules are not implemented. Yet, the  default broadening is used and the broadening functions in `radis/io/exomolapi.py` are commented. Is someone working on it at the moment?

### 💡 Possible implementations
We work on it! 

Minou






---
<sub>*Tip : you can vote for [New Features](https://github.com/radis/radis/issues/79) on https://feathub.com/radis/radis*</sub>



Discussion / hints:
Reported to the Spyder team https://github.com/spyder-ide/spyder/issues/16183 
Hello,
with enabled `parsum_mode="tabulation"` in plot_multi_temperature_fit.py example my spyder 5.0.3 crashes without explanation and running from terminal gives an error with the following traceback:

```
ERROR:MainThread:vaex.execution:error in task, flush task queue
Traceback (most recent call last):
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\execution.py", line 219, in execute_async task._result = parts[0].get_result()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\cpu.py", line 532, in get_result result = result.copy()

numpy.core._exceptions._ArrayMemoryError: Unable to allocate 206. MiB for an array with shape (300, 300, 300) and data type int64

Traceback (most recent call last):
  File "examples/plot_multi_temperature_fit.py", line 70, in <module> s_best, best = sf.fit_spectrum(
  File "c:\users\blekhshtein\radis\radis\lbl\factory.py", line 1860, in fit_spectrum return fit_spectrum(
  File "c:\users\blekhshtein\radis\radis\tools\fitting.py", line 217, in fit_spectrum compute_los_model(
  File "c:\users\blekhshtein\radis\radis\tools\fitting.py", line 211, in <lambda> compute_los_model = lambda fit_parameters: model(
  File "c:\users\blekhshtein\radis\radis\tools\fitting.py", line 131, in Tvib12Tvib3Trot_NonLTEModel s = factory.non_eq_spectrum(**kwargs)
  File "c:\users\blekhshtein\radis\radis\lbl\factory.py", line 1328, in non_eq_spectrum self._calc_populations_noneq_multiTvib(
  File "c:\users\blekhshtein\radis\radis\lbl\base.py", line 2742, in _calc_populations_noneq_multiTvib / Q(Tvib, Trot)
  File "c:\users\blekhshtein\radis\radis\lbl\base.py", line 2655, in Q Q_dict[iso] = parsum.at_noneq_3Tvib(
  File "c:\users\blekhshtein\radis\radis\levels\partfunc.py", line 809, in at_noneq_3Tvib return self._noneq_3Tvib_tabulation_eval(
  File "c:\users\blekhshtein\radis\radis\levels\partfunc.py", line 962, in _noneq_3Tvib_tabulation_eval self._noneq_3Tvib_tabulation_setup(
  File "c:\users\blekhshtein\radis\radis\levels\partfunc.py", line 897, in _noneq_3Tvib_tabulation_setup g_bins_neq = df.sum(
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 943, in sum return self._compute_agg('sum', expression, binby, limits, shape, selection, delay, edges, progress, array_type=array_type)
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 773, in _compute_agg return self._delay(delay, var)
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 1537, in _delay self.execute()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 375, in execute just_run(self.execute_async())
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\asyncio.py", line 35, in just_run return loop.run_until_complete(coro)
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\nest_asyncio.py", line 70, in run_until_complete return f.result()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\asyncio\futures.py", line 178, in result raise self._exception
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\asyncio\tasks.py", line 280, in __step result = coro.send(None)
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 380, in execute_async await self.executor.execute_async()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\execution.py", line 219, in execute_async task._result = parts[0].get_result()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\cpu.py", line 532, in get_result result = result.copy()
numpy.core._exceptions.MemoryError: Unable to allocate 206. MiB for an array with shape (300, 300, 300) and data type int64
```
This error is not related to the vaex-stuck problem above. 
Here are some insights : 

-----

The table used to compute the  partition functions is too large for your system. 

You can adjust it by reducing the [N_bins](https://github.com/radis/radis/blob/74c9883d0063060a3bca799a0753eb03b2312b3a/radis/levels/partfunc.py#L228) parameter. 
It can be done directly in the source code, or  at runtime by changing `sf.parsum_calc[molecule][isotope][electronic_state].N_bins` . 

It will cost you a bit of accuracy. Below, it turns out that reducing from the default 300 to 100 has no visible impact on a typical spectrum (adjust with your typical conditions). But this makes the tabulation grid 3^3 = 27 times smaller. Should fix your memory problems. 

: 
```
from radis import SpectrumFactory
sf = SpectrumFactory(
    2284.2,
    2284.6,
    wstep=0.001,  # cm-1
    pressure=20 * 1e-3,  # bar
    db_use_cached=True,
    lvl_use_cached=True,
    cutoff=1e-25,
    isotope="1,2",
    path_length=10,  # cm-1
    mole_fraction=0.1 * 28.97 / 44.07,
    broadening_max_width=1,  # cm-1
    medium="vacuum",
    export_populations=None,  # 'vib',
    parsum_mode="tabulation"
)
sf.load_databank("HITEMP-CO2-TEST")

Tvib = (1600, 1600, 2500)
Trot = 1200

sf.load_databank('HITEMP-CO2-TEST')        # predefined in ~/radis.
s0 = sf.non_eq_spectrum(Tvib, Trot)

# now, adjust number of bins in partitoin functino calculator : 
for iso in sf._get_isotope_list():
    sf.parsum_calc['CO2'][iso]["X"].N_bins = 100    #👈 reduce this further if you want
s = sf.non_eq_spectrum(Tvib, Trot)

from radis import plot_diff
plot_diff(s0, s)
```
![image](https://user-images.githubusercontent.com/16088743/128718729-4d5c5f1c-287e-456a-959f-98792018a0d0.png)

Thank you for the information.
I thought it is related to the topic, sorry for misposting

All good, this lead to an improvement. We can coarsen the grid for higher dimensions with no impact on accuracy. Proposed in #339 
See https://github.com/spyder-ide/spyder/issues/16183#issuecomment-947664125 for latest developments
Reported to the Spyder team https://github.com/spyder-ide/spyder/issues/16183 
Hello,
with enabled `parsum_mode="tabulation"` in plot_multi_temperature_fit.py example my spyder 5.0.3 crashes without explanation and running from terminal gives an error with the following traceback:

```
ERROR:MainThread:vaex.execution:error in task, flush task queue
Traceback (most recent call last):
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\execution.py", line 219, in execute_async task._result = parts[0].get_result()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\cpu.py", line 532, in get_result result = result.copy()

numpy.core._exceptions._ArrayMemoryError: Unable to allocate 206. MiB for an array with shape (300, 300, 300) and data type int64

Traceback (most recent call last):
  File "examples/plot_multi_temperature_fit.py", line 70, in <module> s_best, best = sf.fit_spectrum(
  File "c:\users\blekhshtein\radis\radis\lbl\factory.py", line 1860, in fit_spectrum return fit_spectrum(
  File "c:\users\blekhshtein\radis\radis\tools\fitting.py", line 217, in fit_spectrum compute_los_model(
  File "c:\users\blekhshtein\radis\radis\tools\fitting.py", line 211, in <lambda> compute_los_model = lambda fit_parameters: model(
  File "c:\users\blekhshtein\radis\radis\tools\fitting.py", line 131, in Tvib12Tvib3Trot_NonLTEModel s = factory.non_eq_spectrum(**kwargs)
  File "c:\users\blekhshtein\radis\radis\lbl\factory.py", line 1328, in non_eq_spectrum self._calc_populations_noneq_multiTvib(
  File "c:\users\blekhshtein\radis\radis\lbl\base.py", line 2742, in _calc_populations_noneq_multiTvib / Q(Tvib, Trot)
  File "c:\users\blekhshtein\radis\radis\lbl\base.py", line 2655, in Q Q_dict[iso] = parsum.at_noneq_3Tvib(
  File "c:\users\blekhshtein\radis\radis\levels\partfunc.py", line 809, in at_noneq_3Tvib return self._noneq_3Tvib_tabulation_eval(
  File "c:\users\blekhshtein\radis\radis\levels\partfunc.py", line 962, in _noneq_3Tvib_tabulation_eval self._noneq_3Tvib_tabulation_setup(
  File "c:\users\blekhshtein\radis\radis\levels\partfunc.py", line 897, in _noneq_3Tvib_tabulation_setup g_bins_neq = df.sum(
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 943, in sum return self._compute_agg('sum', expression, binby, limits, shape, selection, delay, edges, progress, array_type=array_type)
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 773, in _compute_agg return self._delay(delay, var)
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 1537, in _delay self.execute()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 375, in execute just_run(self.execute_async())
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\asyncio.py", line 35, in just_run return loop.run_until_complete(coro)
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\nest_asyncio.py", line 70, in run_until_complete return f.result()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\asyncio\futures.py", line 178, in result raise self._exception
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\asyncio\tasks.py", line 280, in __step result = coro.send(None)
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\dataframe.py", line 380, in execute_async await self.executor.execute_async()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\execution.py", line 219, in execute_async task._result = parts[0].get_result()
  File "C:\Users\Blekhshtein\anaconda3\envs\radis-env\lib\site-packages\vaex\cpu.py", line 532, in get_result result = result.copy()
numpy.core._exceptions.MemoryError: Unable to allocate 206. MiB for an array with shape (300, 300, 300) and data type int64
```
This error is not related to the vaex-stuck problem above. 
Here are some insights : 

-----

The table used to compute the  partition functions is too large for your system. 

You can adjust it by reducing the [N_bins](https://github.com/radis/radis/blob/74c9883d0063060a3bca799a0753eb03b2312b3a/radis/levels/partfunc.py#L228) parameter. 
It can be done directly in the source code, or  at runtime by changing `sf.parsum_calc[molecule][isotope][electronic_state].N_bins` . 

It will cost you a bit of accuracy. Below, it turns out that reducing from the default 300 to 100 has no visible impact on a typical spectrum (adjust with your typical conditions). But this makes the tabulation grid 3^3 = 27 times smaller. Should fix your memory problems. 

: 
```
from radis import SpectrumFactory
sf = SpectrumFactory(
    2284.2,
    2284.6,
    wstep=0.001,  # cm-1
    pressure=20 * 1e-3,  # bar
    db_use_cached=True,
    lvl_use_cached=True,
    cutoff=1e-25,
    isotope="1,2",
    path_length=10,  # cm-1
    mole_fraction=0.1 * 28.97 / 44.07,
    broadening_max_width=1,  # cm-1
    medium="vacuum",
    export_populations=None,  # 'vib',
    parsum_mode="tabulation"
)
sf.load_databank("HITEMP-CO2-TEST")

Tvib = (1600, 1600, 2500)
Trot = 1200

sf.load_databank('HITEMP-CO2-TEST')        # predefined in ~/radis.
s0 = sf.non_eq_spectrum(Tvib, Trot)

# now, adjust number of bins in partitoin functino calculator : 
for iso in sf._get_isotope_list():
    sf.parsum_calc['CO2'][iso]["X"].N_bins = 100    #👈 reduce this further if you want
s = sf.non_eq_spectrum(Tvib, Trot)

from radis import plot_diff
plot_diff(s0, s)
```
![image](https://user-images.githubusercontent.com/16088743/128718729-4d5c5f1c-287e-456a-959f-98792018a0d0.png)

Thank you for the information.
I thought it is related to the topic, sorry for misposting

All good, this lead to an improvement. We can coarsen the grid for higher dimensions with no impact on accuracy. Proposed in #339 
See https://github.com/spyder-ide/spyder/issues/16183#issuecomment-947664125 for latest developments
See SPEARS https://www.sciencedirect.com/science/article/abs/pii/S0022407321004507#, btw 

Fix patch:
diff --git a/docs/lbl/lbl.rst b/docs/lbl/lbl.rst
index daddfcb70..28d9fffd2 100644
--- a/docs/lbl/lbl.rst
+++ b/docs/lbl/lbl.rst
@@ -364,11 +364,15 @@ and :py:mod:`~astropy.units` ::
                          cutoff=1e-25,              # cm/molecule
                          broadening_max_width=10,   # cm-1
                          )
-    sf.load_databank('HITRAN-CO2-TEST')        # this database must be defined in ~/radis.json
+    sf.load_databank('HITRAN-CO2-TEST', load_columns='noneq')        # this database must be defined in ~/radis.json
     s1 = sf.eq_spectrum(Tgas=300 * u.K)
     s2 = sf.eq_spectrum(Tgas=2000 * u.K)
     s3 = sf.non_eq_spectrum(Tvib=2000 * u.K, Trot=300 * u.K)
 
+Note that for non-LTE calculations, specific columns must be loaded. This is done by using the
+``load_columns='noneq'`` parameter. See :py:meth:`~radis.lbl.loader.DatabankLoader.load_databank`
+for more information.
+
 
 .. _label_lbl_config_file:
 
@@ -737,6 +741,21 @@ points for each linewidth.
     but is still experimental in 0.9.30. Feedback welcome!
 
 
+Sparse wavenumber grid
+----------------------
+
+To compute large band spectra with a small number of lines, RADIS includes
+a sparse wavenumber implementation of the DIT algorithm, which is
+activated based on a scarcity criterion (``Nlines/Ngrid_points > 1``).
+
+The sparse version can be forced to be activated or deactivated. This behavior
+is done by setting the `SPARSE_WAVENUMBER` key of the :py:attr:`radis.config`
+dictionary, or of the ~/radis.json user file.
+
+See the :ref:`HITRAN full-range example <example_hitran_full_range>` for an
+example.
+
+
 Database loading
 ----------------
 
diff --git a/docs/spectrum/spectrum.rst b/docs/spectrum/spectrum.rst
index b15398ad4..8ba0b4873 100644
--- a/docs/spectrum/spectrum.rst
+++ b/docs/spectrum/spectrum.rst
@@ -68,6 +68,7 @@ be used to batch-generate multiple spectra using the same line database::
 
     from radis import SpectrumFactory
     sf = SpectrumFactory(...)
+    sf.fetch_databank("hitemp", load_columns='noneq')  # or 'hitran', 'exomol', etc.
     s = sf.eq_spectrum(...)
     s2 = sf.non_eq_spectrum(...)
 
diff --git a/environment.yml b/environment.yml
index de3bd871e..6684da8b7 100644
--- a/environment.yml
+++ b/environment.yml
@@ -19,17 +19,16 @@ dependencies:
 - pandas>=1.0.5
 - plotly>=2.5.1  # for line survey HTML output
 - psutil # for getting user RAM
-- pytables  # for pandas to HDF5 export
-- pytest  # to run test suite
 - scipy>=1.4.0
 - seaborn   # other matplotlib themes
 - termcolor  # terminal colors
 - specutils
+- vaex>=4.4.0   # load HDF5 files  (version for custom HDF5 groups)
 - pip
 - pip:
   - beautifulsoup4 # parse ExoMol website
   - lxml           # parser used for ExoMol website
-  - pyarrow       # for the feather format (temporarily needed for ExoMol)
+  - hjson          # Json with comments (for default_radis.json)
   - publib>=0.3.2  # Plotting styles for Matplotlib
   - hitran-api     # HAPI, used to access TIPS partition functions
   - json-tricks>=3.15.0  # to deal with non jsonable formats
@@ -37,5 +36,4 @@ dependencies:
   - progressbar2    # used in vaex
   - tuna           # to generate visual/interactive performance profiles
   - tables  # for pandas to HDF5 export
-  - vaex    # load HDF5 files
   - habanero  # CrossRef API to retrieve data from doi
diff --git a/examples/plot_exomol_spectrum.py b/examples/plot_exomol_spectrum.py
index 79a5ff06a..67094b716 100644
--- a/examples/plot_exomol_spectrum.py
+++ b/examples/plot_exomol_spectrum.py
@@ -32,7 +32,7 @@
 #%% See line data:
 from radis.io.exomol import fetch_exomol
 
-df = fetch_exomol("SiO", database="EBJT", isotope="1")
+df = fetch_exomol("SiO", database="EBJT", isotope="1", load_wavenum_max=5000)
 print(df)
 
 
diff --git a/examples/plot_hitemp_OH_database.py b/examples/plot_hitemp_OH_database.py
index b4f6e2254..920973167 100644
--- a/examples/plot_hitemp_OH_database.py
+++ b/examples/plot_hitemp_OH_database.py
@@ -31,6 +31,5 @@
 # function. The already downloaded database will be used:
 #
 
-fetch_hitemp("OH", load_wavenum_min=31500, load_wavenum_max=33000, isotope="1").plot(
-    "wav", "int"
-)
+df = fetch_hitemp("OH", load_wavenum_min=31500, load_wavenum_max=33000, isotope="1")
+df.plot("wav", "int")
diff --git a/examples/plot_hitran_full_range.py b/examples/plot_hitran_full_range.py
new file mode 100644
index 000000000..9a5959056
--- /dev/null
+++ b/examples/plot_hitran_full_range.py
@@ -0,0 +1,78 @@
+# -*- coding: utf-8 -*-
+"""
+.. _example_hitran_full_range:
+
+===============================
+Calculate a full range spectrum
+===============================
+
+We compute the transmittance of an homogeneous, 1-km atmosphere layer
+at 300 K and 1 atm, with 420 ppm CO2 and 2% H2O.
+
+This example makes use of a lot of the RADIS optimization and convenience functions :
+
+- The database is retrieved from the latest HITRAN data using [HAPI]_, parsed
+  with RADIS and stored to memory-mapping formats. After the first calls,
+  retrieving the database will be very fast
+
+- The RADIS algorithm includes a sparse wavenumber implementation, which is
+  activated by default in the following example where the waverange is very
+  large. This behavior can be changed by setting the `SPARSE_WAVENUMBER` key
+  of the :py:attr:`radis.config` dictionary, or of the ~/radis.json user file.
+
+- The wavenumber grid is automatically computed using the `wstep='auto'` parameter.
+
+- A linestrength cutoff reduces the initial number of lines.
+
+Eventually, the full-range spectrum (about 4.8M points) is computed within a minute,
+with about 140k lines resolved, i.e. a performance of about 1e10 gridpoints*lines/s.
+
+Extra optimization could be achieved by cutting the spectrum in small intervals,
+but requires an a-priori knowledge of the absorption & emission ranges of the
+spectrum. See the :ref:`Calculate a large spectrum by part <example_large_range_by_part>`
+
+
+"""
+
+
+#%%
+
+import astropy.units as u
+
+from radis import calc_spectrum
+
+s = calc_spectrum(
+    wmin=0.5 * u.um,
+    wmax=15 * u.um,  # cm-1
+    mole_fraction={"CO2": 420e-6, "H2O": 0.02},
+    isotope="1,2,3",
+    pressure=1.01325,  # bar
+    Tgas=300,  # K
+    path_length=1e5,  # 1 km in cm
+    verbose=2,
+    databank="hitran",
+    wstep="auto",
+)
+
+
+#%% Plot low and high resolution spectrum
+s.apply_slit(10, "nm")
+import matplotlib.pyplot as plt
+
+plt.figure(figsize=(16, 6))
+s.plot("transmittance_noslit", wunit="nm", color="k", alpha=0.1, nfig="same")
+s.plot("transmittance", wunit="nm", color="k", nfig="same")
+
+
+#%% Print some details about the computed spectrum :
+print("Number of grid points: ", len(s))
+print(
+    "Number of lines: ", s.c["total_lines"]
+)  # some are discarded by linestrength cutoff
+print("Number of lines: ", s.c["lines_calculated"])
+print("Whether the sparse wavenumber algorithm was activated: ", s.c["sparse_dlm"])
+print(f"Lineshape truncation used: {s.c['truncation']:.1f}cm-1")
+print(f"Total calculation time: {s.c['calculation_time']:.1f}s")
+print(
+    f"Gridpoints * lines/s : {len(s)*s.c['lines_calculated']/s.c['calculation_time']:.1e}"
+)
diff --git a/examples/plot_merged_large_spectrum.py b/examples/plot_merged_large_spectrum.py
index 068abcdb8..3691b79db 100644
--- a/examples/plot_merged_large_spectrum.py
+++ b/examples/plot_merged_large_spectrum.py
@@ -1,12 +1,20 @@
 # -*- coding: utf-8 -*-
 """
-==========================
-Calculate a large spectrum
-==========================
+.. _example_large_range_by_part:
 
-Calculate a full-range spectrum on several partial spectral ranges, and combine them.
+==================================
+Calculate a large spectrum by part
+==================================
+
+In may be faster to calculate a full-range spectrum on several partial spectral ranges,
+and combine them.
+
+Uses the :py:func:`~radis.los.slabs.MergeSlabs` function for that.
+
+Starting from 0.11, RADIS introduces a sparse waverange implementation
+that should make it possible to directly compute a full range spectrum.
+See the :ref:`HITRAN full-range example <example_hitran_full_range>`
 
-Uses the :py:func:`~radis.los.slabs.MergeSlabs` function
 """
 
 from radis import MergeSlabs, calc_spectrum
@@ -32,3 +40,7 @@
 s = MergeSlabs(*spectra, resample="full", out="transparent")
 print(s)
 s.plot("transmittance_noslit", wunit="nm")
+
+import matplotlib.pyplot as plt
+
+plt.ylim(0, 1)
diff --git a/examples/plot_specutils_processing.py b/examples/plot_specutils_processing.py
index 358b4a6ab..54c4c364a 100644
--- a/examples/plot_specutils_processing.py
+++ b/examples/plot_specutils_processing.py
@@ -27,7 +27,7 @@
     .take("radiance")
 )
 s.trim()  # removes nans created by the slit convolution boundary effects
-noise = np.random.normal(0.0, s.max() * 0.03, len(s))
+noise = np.random.normal(0.0, s.max().value * 0.03, len(s))
 s_exp = s + noise
 
 s_exp.plot()
diff --git a/radis/default_radis.json b/radis/default_radis.json
index d503d2194..cf22119cd 100644
--- a/radis/default_radis.json
+++ b/radis/default_radis.json
@@ -1,23 +1,31 @@
+# default values for RADIS configuration parameters
+# the values of your ~/radis.json file will overwrite these.
 {
-    "database": {},
-    "spectroscopic_constants": {
+    "database": {},                          # list of registered databases on this computer. See expected format in https://radis.readthedocs.io/en/latest/lbl/lbl.html#label-lbl-config-file
+    "spectroscopic_constants": {             # path to spectroscopic constants for non-LTE calculations. By default look for them in radis/db/[MOLECULE]/*
         "CO2": "molecules_data.json",
         "CO": "molecules_data.json"
     },
     "plot": {
-        "plotlib": "publib",
-        "context": "paper",
-        "style": "origin"
+        "plotlib": "publib",                 # library for plot style. "publib", "seaborn".
+        "context": "paper",                  # "paper", "talk", "poster"
+        "style": "origin"                    # default style. See Publib (https://github.com/erwanp/publib) or Seaborn (https://seaborn.pydata.org/)
     },
+    # global parameters
+    # -----------------
     "DEBUG_MODE": false,
-    "AUTO_UPDATE_DATABASE": false,
-    "AUTO_UPDATE_SPEC": false,
-    "OLDEST_COMPATIBLE_VERSION": ",0.9.1",
-    "USE_CYTHON": true,
-    "GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD": 3,
-    "GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD": 1,
-    "SPARSE_WAVERANGE": "auto",
-    "molparams": {
+    "AUTO_UPDATE_DATABASE": false,           # auto-convert 'pytables' files to 'vaex'
+    "AUTO_UPDATE_SPEC": false,               # auto-update saved Spectrum objects to latest format.
+    "OLDEST_COMPATIBLE_VERSION": "0.9.1",    # automatically regenerate cache files generated with versions anterior to this one
+    "USE_CYTHON": true,                      # use Cython module if available (else, default to Python)
+    "GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD": 3,    # raise a warning if less than THIS number of grid points per lineshape
+    "GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD": 1,   # raise an error if less than THIS number of grid points per lineshape
+    "MEMORY_MAPPING_ENGINE": "vaex",         #  warning: on Spyder this may result in freezes. see https://github.com/spyder-ide/spyder/issues/16183.
+    "SPARSE_WAVERANGE": "auto",              # sparse DLM algorithm. May be smaller on dense spectra. If "auto", a scarcity criterion is used (Nlines/Ngrids > 1)
+
+    # molecular parameters
+    # --------------------
+    "molparams": {                           # molecular parameters for non-HITRAN species (the latter being defined in https://github.com/radis/radis/blob/master/radis/db/molparam.txt)
         "isotope_full_name": {
             "AlCl": {
                 "1": "27Al-35Cl",
@@ -829,4 +837,4 @@
             }
         }
     }
-}
\ No newline at end of file
+}
diff --git a/radis/io/cache_files.py b/radis/io/cache_files.py
index 56b7dbc6f..f721918fd 100644
--- a/radis/io/cache_files.py
+++ b/radis/io/cache_files.py
@@ -55,6 +55,7 @@
 def load_h5_cache_file(
     cachefile,
     use_cached,
+    columns=None,
     valid_if_metadata_is={},
     relevant_if_metadata_above={},
     relevant_if_metadata_below={},
@@ -79,6 +80,8 @@ def load_h5_cache_file(
         if using the cache file, check if the file is deprecated. If it
         is deprecated, regenerate the file unless ``'force'`` was used
         (in that case, raise an error)
+    columns: list, or ``None``
+        columns to load
     valid_if_metadata_is: dict
         values are compared to cache file attributes. If they dont match,
         the file is considered deprecated. See ``use_cached`` to know
@@ -177,7 +180,7 @@ def load_h5_cache_file(
     try:
         # Load file :
         manager = HDF5Manager(engine)
-        df = manager.load(cachefile, key="df")
+        df = manager.load(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
@@ -236,7 +239,7 @@ def check_cache_file(
     expected_metadata={},
     compare_as_close=[],
     verbose=True,
-    engine="auto",
+    engine="guess",
 ):
     """Quick function that check status of cache file generated by RADIS:
 
@@ -265,6 +268,8 @@ def check_cache_file(
         compare with ``np.isclose(a,b)`` rather than ``a==b``
     verbose: boolean
         print stuff
+     engine: ``'h5py'``, ``'pytables'``, ``'vaex'``, ``'guess'``
+        which HDF5 library to use. If ``'guess'``, try to guess.
 
     Returns
     -------
@@ -330,7 +335,7 @@ def check_not_deprecated(
     compare_as_close=[],
     current_version=None,
     last_compatible_version=radis.config["OLDEST_COMPATIBLE_VERSION"],
-    engine="auto",
+    engine="guess",
 ):
     """Make sure cache file is not deprecated: checks that ``metadata`` is the
     same, and that the version under which the file was generated is valid.
@@ -359,11 +364,10 @@ def check_not_deprecated(
         If the file was generated in a non-compatible version, an error is raised.
         (useful parameter to force regeneration of certain cache files after a
          breaking change in a new version)
-    engine: ``'h5py'``, ``'pytables'``, ``'vaex'``, ``'auto'``
-        which HDF5 library to use. If ``'auto'``, try to guess. Note: ``'vaex'``
-        uses ``'h5py'`` compatible HDF5.
+    engine: ``'h5py'``, ``'pytables'``, ``'vaex'``, ``'guess'``
+        which HDF5 library to use. If ``'guess'``, try to guess.
     """
-    if engine == "auto":
+    if engine == "guess":
         engine = HDF5Manager.guess_engine(file)
 
     # Get metadata :
@@ -482,47 +486,46 @@ def check_relevancy(
     relevant_if_metadata_above,
     relevant_if_metadata_below,
     verbose=True,
-    key="df",
-    engine="auto",
+    key="default",
+    engine="guess",
 ):
     """Make sure cache file is relevant.
 
-     Use case: checks that  wavenumber min and
-     wavenumber max in ``metadata`` are relevant for the specified spectral
-     range.
+    Use case: checks that  wavenumber min and
+    wavenumber max in ``metadata`` are relevant for the specified spectral
+    range.
 
-     Parameters
-     ----------
-     file: str
-         a `` .h5``  line database cache file
-     load_only_wavenum_above, relevant_if_metadata_below: dict
-         only load the cached file if the metadata values are above/below
-         the specific values for each key.
-     relevant_if_metadata_above, relevant_if_metadata_below: dict
-         file is relevant if the file metadata value for each key of the dictionary
-         is above/below the value in the dictionary
+    Parameters
+    ----------
+    file: str
+        a `` .h5``  line database cache file
+    load_only_wavenum_above, relevant_if_metadata_below: dict
+        only load the cached file if the metadata values are above/below
+        the specific values for each key.
+    relevant_if_metadata_above, relevant_if_metadata_below: dict
+        file is relevant if the file metadata value for each key of the dictionary
+        is above/below the value in the dictionary
 
-     Other Parameters
-     ----------------
-     key: str
-         dataset key in storer.
-    engine: ``'h5py'``, ``'pytables'``, ``'vaex'``, ``'auto'``
-        which HDF5 library to use. If ``'auto'``, try to guess. Note: ``'vaex'``
-        uses ``'h5py'`` compatible HDF5.
+    Other Parameters
+    ----------------
+    key: str
+        dataset key in storer.
+    engine: ``'h5py'``, ``'pytables'``, ``'vaex'``, ``'guess'``
+       which HDF5 library to use. If ``'guess'``, try to guess.
 
-     Examples
-     --------
-     You want to compute a spectrum in between 2300 and 2500 cm-1. A line database
-     file is relevant only if its metadata says that ``'wavenum_max' > 2300``
-     and ``'wavenum_min'`` < 2500 cm-1.
+    Examples
+    --------
+    You want to compute a spectrum in between 2300 and 2500 cm-1. A line database
+    file is relevant only if its metadata says that ``'wavenum_max' > 2300``
+    and ``'wavenum_min'`` < 2500 cm-1.
 
-         check_relevancy('path/to/file', relevant_if_metadata_above={'wavenum_max':2300},
-                         relevant_if_metadata_below={'wavenum_min':2500})
+        check_relevancy('path/to/file', relevant_if_metadata_above={'wavenum_max':2300},
+                        relevant_if_metadata_below={'wavenum_min':2500})
 
-         the specified value.
+        the specified value.
 
     """
-    if engine == "auto":
+    if engine == "guess":
         engine = HDF5Manager.guess_engine(file)
 
     # Get metadata :
@@ -564,7 +567,7 @@ def save_to_hdf(
     fname,
     metadata,
     version=None,
-    key="df",
+    key="default",
     overwrite=True,
     verbose=True,
     engine="pytables",
@@ -573,8 +576,8 @@ def save_to_hdf(
 
      Parameters
      ----------
-     df: a pandas DataFrame
-         data will be stored in the key ``'df'``
+     df: a pandas/vaex DataFrame
+         data will be stored in this key.
      fname: str
          ``.h5`` file where to store.
      metadata: dict
@@ -690,7 +693,7 @@ def cache_file_name(fname, engine="pytables"):
     engine: ``'h5py'``, ``'pytables'``, ``'vaex'``
        which HDF5 library to use. Default ``pytables``
     """
-    if engine in ["pytables"]:
+    if engine in ["pytables", "pytables-fixed"]:
         return splitext(fname)[0] + ".h5"
     elif engine in ["h5py", "vaex"]:
         return splitext(fname)[0] + ".hdf5"
diff --git a/radis/io/cdsd.py b/radis/io/cdsd.py
index 99a5a3c89..f40c57774 100644
--- a/radis/io/cdsd.py
+++ b/radis/io/cdsd.py
@@ -116,10 +116,12 @@ def cdsd2df(
     fname,
     version="hitemp",
     cache=True,
+    load_columns=None,
     verbose=True,
     drop_non_numeric=True,
     load_wavenum_min=None,
     load_wavenum_max=None,
+    engine="pytables",
 ):
     """Convert a CDSD-HITEMP [1]_ or CDSD-4000 [2]_ file to a Pandas dataframe.
 
@@ -135,6 +137,12 @@ def cdsd2df(
         improves performances a lot (but changes in the database are not
         taken into account). If ``False``, no database is used. If 'regen', temp
         file are reconstructed. Default ``True``.
+    load_columns: list
+        columns to load. If ``None``, loads everything
+
+        .. note::
+            this is only relevant if loading from a cache file. To generate
+            the cache file, all columns are loaded anyway.
 
     Other Parameters
     ----------------
@@ -147,6 +155,8 @@ def cdsd2df(
         if not ``'None'``, only load the cached file if it contains data for
         wavenumbers above/below the specified value. See :py:func`~radis.io.cache_files.load_h5_cache_file`.
         Default ``'None'``.
+    engine: 'pytables', 'vaex'
+        format for Hdf5 cache file. Default `pytables`
 
     Returns
     -------
@@ -231,7 +241,7 @@ def cdsd2df(
         raise ValueError("Unknown CDSD version: {0}".format(version))
 
     # Use cache file if possible
-    fcache = cache_file_name(fname)
+    fcache = cache_file_name(fname, engine=engine)
     if cache and exists(fcache):
         relevant_if_metadata_above = (
             {"wavenum_max": load_wavenum_min} if load_wavenum_min else {}
@@ -242,12 +252,14 @@ def cdsd2df(
         df = load_h5_cache_file(
             fcache,
             cache,
+            columns=load_columns,
             valid_if_metadata_is=metadata,
             relevant_if_metadata_above=relevant_if_metadata_above,
             relevant_if_metadata_below=relevant_if_metadata_below,
             current_version=radis.__version__,
             last_compatible_version=radis.config["OLDEST_COMPATIBLE_VERSION"],
             verbose=verbose,
+            engine=engine,
         )
         if df is not None:
             return df
@@ -284,6 +296,7 @@ def cdsd2df(
                 key="df",
                 overwrite=True,
                 verbose=verbose,
+                engine=engine,
             )
         except PermissionError:
             if verbose:
diff --git a/radis/io/dbmanager.py b/radis/io/dbmanager.py
index 820405e38..4528bf28b 100644
--- a/radis/io/dbmanager.py
+++ b/radis/io/dbmanager.py
@@ -22,7 +22,7 @@
     from radis.io.cache_files import check_not_deprecated
 
 from datetime import date
-from os.path import join
+from os.path import join, split
 
 import numpy as np
 import pandas as pd
@@ -58,6 +58,9 @@ class DatabaseManager(object):
     molecule: str
     local_databases: str
         path to local database
+    engine: 'vaex', 'pytables', 'h5py', or 'default'
+        memory-mapping library to use with this database. If 'default' use
+        the value from ~/radis.json
 
     Other Parameters
     ----------------
@@ -88,19 +91,38 @@ def __init__(
         name,
         molecule,
         local_databases,
+        engine,
         verbose=False,
         parallel=True,
         nJobs=-2,
         batch_size="auto",
     ):
+        if engine == "default":
+            from radis import config
+
+            engine = config["MEMORY_MAPPING_ENGINE"]  # 'pytables', 'vaex', 'feather'
+
+        # vaex processes are stuck if ran from Spyder. See https://github.com/spyder-ide/spyder/issues/16183
+        if engine == "vaex" and any("SPYDER" in name for name in os.environ):
+            from radis.misc.log import printwarn
+
+            printwarn(
+                "Spyder IDE detected while using memory_mapping_engine='vaex'.\nVaex is the fastest way to read database files in RADIS, but Vaex processes may be stuck if ran from Spyder. See https://github.com/spyder-ide/spyder/issues/16183. You may consider using another IDE, or using a different `memory_mapping_engine` such as 'pytables' or 'feather'. You can change the engine in Spectrum.fetch_databank() calls, or globally by setting the 'MEMORY_MAPPING_ENGINE' key in your ~/radis.json  (note: starting another iPython console somehow releases the freeze in Spyder)  \n"
+            )
 
         self.name = name
         self.molecule = molecule
         self.local_databases = local_databases
+        # create folder if needed
+        from radis.misc.basics import make_folders
+
+        make_folders(*split(abspath(local_databases)))
+
         self.downloadable = False  # by default
         self.format = ""
+        self.engine = engine
 
-        self.tempdir = join(self.local_databases, "downloads(can_be_deleted)")
+        self.tempdir = join(self.local_databases, "downloads__can_be_deleted")
         self.ds = DataSource(self.tempdir)
 
         self.verbose = verbose
@@ -112,7 +134,7 @@ def __init__(
             4  #: type: int. If there are less files, don't use parallel mode.
         )
 
-    def get_filenames(self, engine):
+    def get_filenames(self):
         """Get names of all files in the database (even if not downloaded yet)
 
         See Also
@@ -120,6 +142,7 @@ def get_filenames(self, engine):
         :py:meth:`~radis.io.linedb.get_files_to_download`"""
         verbose = self.verbose
         local_databases = self.local_databases
+        engine = self.engine
 
         # First; checking if the database is registered in radis.json
         if self.is_registered():
@@ -132,7 +155,7 @@ def get_filenames(self, engine):
                 )
             except AssertionError as err:
                 raise DeprecatedFileWarning(
-                    "Database entry {self.name} not valid anymore. See above."
+                    f"Database entry {self.name} not valid anymore. See above which keys are missing/wrong. Update or delete the entry in your ~/radis.json"
                 ) from err
             else:
                 local_files = entries["path"]
@@ -190,11 +213,12 @@ def fetch_urlnames(self) -> list:
             "This function should be overwritten by the DatabaseManager subclass"
         )
 
-    def check_deprecated_files(self, local_files, engine, remove=True):
+    def check_deprecated_files(self, local_files, auto_remove=True):
         """Check metadata of files and remove the deprecated ones
 
-        Unless remove=False: Then raise an error"""
+        Unless auto_remove=False: Then raise an error"""
         verbose = self.verbose
+        engine = self.engine
         for local_file in local_files:
             try:
                 check_not_deprecated(
@@ -204,7 +228,7 @@ def check_deprecated_files(self, local_files, engine, remove=True):
                     engine=engine,
                 )
             except DeprecatedFileWarning as err:
-                if not remove:
+                if not auto_remove:
                     raise err
                 else:  # delete file to regenerate it in the end of the script
                     if verbose:
@@ -268,11 +292,11 @@ def register(self, dict_entries):
             self.verbose,
         )
 
-    def get_hdf5_manager(self, engine):
-        return HDF5Manager(engine=engine)
+    def get_hdf5_manager(self):
+        return HDF5Manager(engine=self.engine)
 
-    def download_and_parse(self, urlnames, local_files, engine="pytables"):
-        all_local_files, _ = self.get_filenames(engine)
+    def download_and_parse(self, urlnames, local_files):
+        all_local_files, _ = self.get_filenames()
 
         verbose = self.verbose
         molecule = self.molecule
@@ -317,7 +341,6 @@ def download_and_parse_one_file(urlname, local_file, Ndownload):
                 pbar_Ntot_estimate_factor=pbar_Ntot_estimate_factor,
                 pbar_Nlines_already=Nlines_total,
                 pbar_last=(Ndownload == Ntotal_downloads),
-                engine=engine,
             )
             # except Exception as err:
             #     raise IOError("Problem parsing `{0}`. Check the error above. It may arise if the file wasn't properly downloaded. Try to delete it".format(self.ds._findfile(urlname))) from err
@@ -374,7 +397,6 @@ def load(
         columns,
         load_wavenum_min,
         load_wavenum_max,
-        engine="pytables",
     ):
         """
         Other Parameters
@@ -382,6 +404,7 @@ def load(
         columns: list of str
             list of columns to load. If ``None``, returns all columns in the file.
         """
+        engine = self.engine
         if engine == "pytables":
             df_all = []
             for local_file in local_files:
@@ -393,7 +416,7 @@ def load(
                         load_wavenum_min=load_wavenum_min,
                         load_wavenum_max=load_wavenum_max,
                         verbose=self.verbose,
-                        engine="pytables",
+                        engine=engine,
                     )
                 )
             return pd.concat(df_all)
@@ -407,7 +430,7 @@ def load(
                 load_wavenum_min=load_wavenum_min,
                 load_wavenum_max=load_wavenum_max,
                 verbose=self.verbose,
-                engine="vaex",
+                engine=engine,
             )
         else:
             raise NotImplementedError(engine)
@@ -420,12 +443,22 @@ def plot(self, local_files, isotope, wavenum_min, wavenum_max):
         )
         df.plot("wav", "int")
 
-    def get_nrows(self, local_file, engine="pytables"):
+    def get_nrows(self, local_file):
+        """ Get number of rows (without loading all DataFrame)"""
+        engine = self.engine
         if engine == "pytables":
             with pd.HDFStore(local_file, "r") as store:
                 nrows = store.get_storer("df").nrows
 
-        elif engine in ["vaex", "h5py"]:
+        elif engine == "vaex":
+            import vaex
+
+            # by default vaex does not load everything
+            df = vaex.open(local_file)
+            nrows = len(df)
+            df.close()
+
+        elif engine in ["h5py"]:
             raise NotImplementedError
         else:
             raise ValueError(engine)
diff --git a/radis/io/exomol.py b/radis/io/exomol.py
index 763979a15..ea9e58933 100644
--- a/radis/io/exomol.py
+++ b/radis/io/exomol.py
@@ -218,11 +218,13 @@ def fetch_exomol(
     isotope="1",
     load_wavenum_min=None,
     load_wavenum_max=None,
+    columns=None,
     cache=True,
     verbose=True,
     clean_cache_files=True,
     return_local_path=False,
     return_partition_function=False,
+    engine="default",
 ):
     """Stream ExoMol file from EXOMOL website. Unzip and build a HDF5 file directly.
 
@@ -261,6 +263,8 @@ def fetch_exomol(
 
     load_wavenum_min, load_wavenum_max: float (cm-1)
         load only specific wavenumbers.
+    columns: list of str
+        list of columns to load. If ``None``, returns all columns in the file.
 
     Other Parameters
     ----------------
@@ -274,6 +278,8 @@ def fetch_exomol(
         if ``True``, also returns the path of the local database file.
     return_partition_function: bool
         if ``True``, also returns a :py:class:`~radis.levels.partfunc.PartFuncExoMol` object.
+    engine: 'vaex', 'feather'
+        which memory-mapping library to use. If 'default' use the value from ~/radis.json
 
     Returns
     -------
@@ -303,6 +309,9 @@ def fetch_exomol(
     :py:func:`~radis.io.hdf5.hdf2df`
 
     """
+    # TODO: implement columns= ... to load only specific columns.
+    # refactor with "self._quantumNumbers" (which serves the same purpose)
+
     # Ensure isotope format:
     try:
         isotope = int(isotope)
@@ -318,7 +327,7 @@ def fetch_exomol(
     )
 
     _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"
-    if database is None:
+    if database is None or database == "default":
         if len(known_exomol_databases) == 1:
             database = known_exomol_databases[0]
             if verbose:
@@ -350,7 +359,7 @@ def fetch_exomol(
         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])
+    mdb = MdbExomol(local_path, [load_wavenum_min, load_wavenum_max], engine=engine)
 
     # Attributes of the DataFrame
     from radis.db.classes import HITRAN_MOLECULES
@@ -412,6 +421,12 @@ class MdbExomol(object):
 
     """
 
+    # TODO : inherit from DatabaseManager or similar
+
+    # @dev: In exojax this class is defined in exojax/spec/moldb.py
+    # see https://github.com/HajimeKawahara/exojax/blob/develop/src/exojax/spec/moldb.py
+    # It uses Jax arrays (jnp). Here RADIS uses Numpy arrays.
+
     def __init__(
         self,
         path,
@@ -420,6 +435,7 @@ def __init__(
         crit=-np.inf,
         bkgdatm="H2",
         broadf=True,
+        engine="vaex",
     ):
         """Molecular database for Exomol form
 
@@ -435,7 +451,21 @@ def __init__(
            The trans/states files can be very large. For the first time to read it, we convert it to the feather-format. After the second-time, we use the feather format instead.
 
         """
-        explanation = "Note: Couldn't find the feather format. We convert data to the feather format. After the second time, it will become much faster."
+        if engine == "default":
+            import radis
+
+            engine = radis.config["MEMORY_MAPPING_ENGINE"]
+
+        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)
         t0 = self.path.parents[0].stem
@@ -498,12 +528,26 @@ def __init__(
             self.alpha_ref_def = 0.07
 
         # load states
-        if self.states_file.with_suffix(".feather").exists():
-            states = pd.read_feather(self.states_file.with_suffix(".feather"))
-        else:
-            print(explanation)
-            states = exomolapi.read_states(self.states_file, dic_def)
-            states.to_feather(self.states_file.with_suffix(".feather"))
+        if engine == "feather":
+            if self.states_file.with_suffix(".feather").exists():
+                ndstates = 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."
+                )
+                ndstates = exomolapi.read_states(self.states_file, dic_def)
+                ndstates.to_feather(self.states_file.with_suffix(".feather"))
+        elif engine == "vaex":
+            if self.states_file.with_suffix(".bz2.hdf5").exists():
+                states = vaex.open(self.states_file.with_suffix(".bz2.hdf5"))
+                ndstates = vaex.array_types.to_numpy(states)
+            else:
+                print(
+                    "Note: Caching states data to the hdf5 format with vaex. After the second time, it will become much faster."
+                )
+                states = exomolapi.read_states(self.states_file, dic_def)
+                ndstates = vaex.array_types.to_numpy(states)
+
         # load pf
         pf = exomolapi.read_pf(self.pf_file)
         self.gQT = pf["QT"].to_numpy()  # grid QT
@@ -516,15 +560,41 @@ def __init__(
             if not self.trans_file.exists():
                 self.download(molec, [".trans.bz2"])
 
-            if self.trans_file.with_suffix(".feather").exists():
-                trans = pd.read_feather(self.trans_file.with_suffix(".feather"))
-            else:
-                print(explanation)
-                trans = exomolapi.read_trans(self.trans_file)
-                trans.to_feather(self.trans_file.with_suffix(".feather"))
+            if engine == "vaex":
+                if self.trans_file.with_suffix(".hdf5").exists():
+                    trans = vaex.open(self.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]
+                    ndtrans = vaex.array_types.to_numpy(trans)
+
+                    # 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)
+                    ndtrans = vaex.array_types.to_numpy(trans)
+
+                    # mask needs to be applied
+                    mask_needed = True
+            elif engine == "feather":
+                if self.trans_file.with_suffix(".feather").exists():
+                    ndtrans = 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."
+                    )
+                    ndtrans = exomolapi.read_trans(self.trans_file)
+                    ndtrans.to_feather(self.trans_file.with_suffix(".feather"))
 
-            #!!TODO:restrict NOW the trans size to avoid useless overload of memory and CPU
-            # trans = trans[(trans['nu'] > self.nurange[0] - self.margin) & (trans['nu'] < self.nurange[1] + self.margin)]
             # compute gup and elower
             (
                 self._A,
@@ -533,9 +603,46 @@ def __init__(
                 self._gpp,
                 self._jlower,
                 self._jupper,
+                mask_zeronu,
                 self._quantumNumbers,
-            ) = exomolapi.pickup_gE(states, trans, dic_def)
-        else:
+            ) = exomolapi.pickup_gE(ndstates, ndtrans, self.trans_file, dic_def)
+
+            # Compute linestrengths or retrieve them from cache
+            self.Tref = 296.0
+            self.QTref = np.array(self.QT_interp(self.Tref))
+
+            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]
+            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,
+                )
+
+                # 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["nu_positive"] = mask_zeronu
+                trans = trans[trans.nu_positive].extract()
+                trans.drop("nu_positive", inplace=True)
+
+                trans["nu_lines"] = self.nu_lines
+                trans["Sij0"] = self.Sij0
+
+                if engine == "vaex":
+                    trans.export(self.trans_file.with_suffix(".hdf5"))
+
+        else:  # dic_def["numinf"] is not None
             imin = (
                 np.searchsorted(dic_def["numinf"], nurange[0], side="right") - 1
             )  # left side
@@ -551,15 +658,42 @@ def __init__(
                     self.download(
                         molec, extension=[".trans.bz2"], numtag=dic_def["numtag"][i]
                     )
-                if trans_file.with_suffix(".feather").exists():
-                    trans = pd.read_feather(trans_file.with_suffix(".feather"))
-                else:
-                    print(explanation)
-                    trans = exomolapi.read_trans(trans_file)
-                    trans.to_feather(trans_file.with_suffix(".feather"))
-                    #!!TODO:restrict NOW the trans size to avoid useless overload of memory and CPU
-                    # trans = trans[(trans['nu'] > self.nurange[0] - self.margin) & (trans['nu'] < self.nurange[1] + self.margin)]
-
+                if engine == "feather":
+                    if trans_file.with_suffix(".feather").exists():
+                        trans = pd.read_feather(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."
+                        )
+                        ndtrans = exomolapi.read_trans(trans_file)
+                        ndtrans.to_feather(trans_file.with_suffix(".feather"))
+                        #!!TODO:restrict NOW the trans size to avoid useless overload of memory and CPU
+                        # trans = trans[(trans['nu'] > self.nurange[0] - self.margin) & (trans['nu'] < self.nurange[1] + self.margin)]
+                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]
+                        ndtrans = vaex.array_types.to_numpy(trans)
+
+                        # 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)
+                        ndtrans = vaex.array_types.to_numpy(trans)
+
+                        # mask needs to be applied
+                        mask_needed = True
                 self.trans_file.append(trans_file)
 
                 # compute gup and elower
@@ -571,9 +705,42 @@ def __init__(
                         self._gpp,
                         self._jlower,
                         self._jupper,
+                        mask_zeronu,
                         self._quantumNumbers,
-                    ) = exomolapi.pickup_gE(states, trans, dic_def)
-                else:
+                    ) = exomolapi.pickup_gE(ndstates, ndtrans, 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]
+                    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,
+                        )
+
+                        # 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["nu_positive"] = mask_zeronu
+                        trans = trans[trans.nu_positive].extract()
+                        trans.drop("nu_positive", inplace=True)
+
+                        trans["nu_lines"] = self.nu_lines
+                        trans["Sij0"] = self.Sij0
+                else:  # k!=0
+
                     (
                         Ax,
                         nulx,
@@ -581,8 +748,35 @@ def __init__(
                         gppx,
                         jlowerx,
                         jupperx,
-                        _quantumNumbersx,
-                    ) = exomolapi.pickup_gE(states, trans, dic_def)
+                        mask_zeronu,
+                        quantumNumbersx,
+                    ) = exomolapi.pickup_gE(ndstates, ndtrans, trans_file)
+                    if engine == "vaex" and trans_file.with_suffix(".hdf5").exists():
+                        Sij0x = ndtrans[:, 4]
+                    else:
+                        ##Line strength:
+                        from radis.lbl.base import (  # TODO: move elsewhere
+                            linestrength_from_Einstein,
+                        )
+
+                        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,
+                        )
+
+                        # 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["nu_positive"] = mask_zeronu
+                        trans = trans[trans.nu_positive].extract()
+                        trans.drop("nu_positive", inplace=True)
+
+                        trans["nu_lines"] = nulx
+                        trans["Sij0"] = Sij0x
+
                     self._A = np.hstack([self._A, Ax])
                     self.nu_lines = np.hstack([self.nu_lines, nulx])
                     self._elower = np.hstack([self._elower, elowerx])
@@ -590,23 +784,14 @@ def __init__(
                     self._jlower = np.hstack([self._jlower, jlowerx])
                     self._jupper = np.hstack([self._jupper, jupperx])
                     self._quantumNumbers = np.hstack(
-                        [self._quantumNumbers, _quantumNumbersx]
+                        [self._quantumNumbers, quantumNumbersx]
                     )
 
-        self.Tref = 296.0
-        self.QTref = np.array(self.QT_interp(self.Tref))
-
-        from radis.lbl.base import linestrength_from_Einstein  # TODO: move elsewhere
-
-        self.Sij0 = linestrength_from_Einstein(
-            A=self._A,
-            gu=self._gpp,
-            El=self._elower,
-            Ia=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 (
+                        engine == "vaex"
+                        and not trans_file.with_suffix(".hdf5").exists()
+                    ):
+                        trans.export(trans_file.with_suffix(".hdf5"))
 
         ### MASKING ###
         mask = (
@@ -615,27 +800,31 @@ def __init__(
             * (self.Sij0 > self.crit)
         )
 
-        self.masking(mask)
+        self.masking(mask, mask_needed)
 
-    def masking(self, mask):
+    def masking(self, mask, mask_needed=True):
         """applying mask and (re)generate jnp.arrays
 
         Args:
            mask: mask to be applied. self.mask is updated.
+           mask_needed: whether mask needs to be applied or not
 
         """
-        # TODO : replace with HDF5 masking ?
+        # TODO : with Vaex-HDF5, selection should happen on disk.
 
         # numpy float 64 Do not convert them jnp array
-        self.nu_lines = self.nu_lines[mask]
-        self.Sij0 = self.Sij0[mask]
-        self._A = self._A[mask]
-        self._elower = self._elower[mask]
-        self._gpp = self._gpp[mask]
-        self._jlower = self._jlower[mask]
-        self._jupper = self._jupper[mask]
-
-        # jnp arraysgamma_
+        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)
@@ -647,9 +836,6 @@ def masking(self, mask):
         ##Broadening parameters
         self.set_broadening()
 
-        for key, array in self._quantumNumbers.items():
-            self._quantumNumbers[key] = array[mask]
-
     def set_broadening(self, alpha_ref_def=None, n_Texp_def=None):
         """setting broadening parameters
 
diff --git a/radis/io/exomolapi.py b/radis/io/exomolapi.py
index 405faf283..326ec31c5 100644
--- a/radis/io/exomolapi.py
+++ b/radis/io/exomolapi.py
@@ -111,7 +111,7 @@ def read_pf(pff):
     return dat
 
 
-def read_trans(transf):
+def read_trans(transf, engine="vaex"):
     """Exomol IO for a transition file
     Note:
         i_upper=Upper state counting number
@@ -122,26 +122,46 @@ def read_trans(transf):
 
     Args:
         transf: transition file
+        engine: parsing engine to use ('vaex', 'csv')
     Returns:
         transition data in pandas DataFrame
 
     """
-    try:
-        dat = pd.read_csv(
-            transf,
-            compression="bz2",
-            sep="\s+",
-            names=("i_upper", "i_lower", "A", "nu_lines"),
-        )
-    except:
-        dat = pd.read_csv(
-            transf, sep="\s+", names=("i_upper", "i_lower", "A", "nu_lines")
-        )
+    if engine == "vaex":
+        import vaex
+
+        try:
+            dat = vaex.from_csv(
+                transf,
+                compression="bz2",
+                sep="\s+",
+                names=("i_upper", "i_lower", "A", "nu_lines"),
+                convert=True,
+            )
+        except:
+            dat = vaex.read_csv(
+                transf,
+                sep="\s+",
+                names=("i_upper", "i_lower", "A", "nu_lines"),
+                convert=True,
+            )
+    elif engine == "csv":
+        try:
+            dat = pd.read_csv(
+                transf,
+                compression="bz2",
+                sep="\s+",
+                names=("i_upper", "i_lower", "A", "nu_lines"),
+            )
+        except:
+            dat = pd.read_csv(
+                transf, sep="\s+", names=("i_upper", "i_lower", "A", "nu_lines")
+            )
 
     return dat
 
 
-def read_states(statesf, dic_def):
+def read_states(statesf, dic_def, engine="vaex"):
     """Exomol IO for a state file
     Note:
         i=state counting number
@@ -153,6 +173,7 @@ def read_states(statesf, dic_def):
     Args:
         statesf: state file
         dic_def: Info from def file to read extra quantum numbers
+        engine: parsing engine to use ('vaex', 'csv')
     Returns:
         states data in pandas DataFrame
 
@@ -163,41 +184,63 @@ def read_states(statesf, dic_def):
     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((np.arange(4), 5 + 2 + np.arange(N)))
-        names = ("i", "E", "g", "J", "Lande", "lifetime") + tuple(quantum_labels)
+        usecol = np.concatenate((usecol, 5 + 2 + np.arange(N)))
+        names = names + ("Lande", "lifetime") + tuple(quantum_labels)
     elif dic_def["Landé"]:
-        usecol = np.concatenate((np.arange(4), 5 + 1 + np.arange(N)))
-        names = ("i", "E", "g", "J", "Lande") + tuple(quantum_labels)
+        usecol = np.concatenate((usecol, 5 + 1 + np.arange(N)))
+        names = names + ("Lande") + tuple(quantum_labels)
     elif dic_def["lifetime"]:
-        usecol = np.concatenate((np.arange(4), 5 + 1 + np.arange(N)))
-        names = ("i", "E", "g", "J", "lifetime") + tuple(quantum_labels)
+        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((np.arange(4), 5 + np.arange(N)))
-        names = ("i", "E", "g", "J") + tuple(quantum_labels)
-
-    try:
-        dat = pd.read_csv(
-            statesf, compression="bz2", sep="\s+", usecols=usecol, names=names
-        )
-    except:  #!!!TODO What was the expected error?
-        dat = pd.read_csv(statesf, sep="\s+", usecols=usecol, names=names)
+        usecol = np.concatenate((usecol, 5 + np.arange(N)))
+        names = names + tuple(quantum_labels)
+
+    if engine == "vaex":
+        import vaex
+
+        try:
+            dat = vaex.from_csv(
+                statesf,
+                compression="bz2",
+                sep="\s+",
+                usecols=usecol,
+                names=names,
+                convert=True,
+            )
+        except:
+            dat = vaex.read_csv(
+                statesf, sep="\s+", usecols=usecol, names=names, convert=True
+            )
+    elif engine == "csv":
+        try:
+            dat = pd.read_csv(
+                statesf, compression="bz2", sep="\s+", usecols=usecol, names=names
+            )
+        except:  #!!!TODO What was the expected error?
+            dat = pd.read_csv(statesf, sep="\s+", usecols=usecol, names=names)
+    else:
+        raise NotImplementedError(engine)
 
     return dat
 
 
-def pickup_gE(states, trans, dic_def, trans_lines=False):
+def pickup_gE(ndstates, ndtrans, trans_file, dic_def, trans_lines=False):
     """extract g_upper (gup), E_lower (elower), and J_lower and J_upper from states DataFrame and insert them to transition DataFrame.
 
     Args:
-       states: states pandas DataFrame
-       trans: transition pandas DataFrame
+       ndstates: states numpy array    - the i, E, g, J are in the 4 first columns
+       ndtrans: 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
 
 
     Returns:
-       A, nu_lines, elower, gup, jlower, jupper, **quantum_labels
+       A, nu_lines, elower, gup, jlower, jupper, mask, **quantum_labels
 
     Note:
        We first convert pandas DataFrame to ndarray. The state counting numbers in states DataFrame is used as indices of the new array for states (newstates). We remove the state count numbers as the column of newstate, i.e. newstates[:,k] k=0: E, 1: g, 2: J. Then, we can directly use the state counting numbers as mask.
@@ -205,8 +248,6 @@ def pickup_gE(states, trans, dic_def, trans_lines=False):
 
     """
     ### Step 1. Essential quantum number for spectra
-    ndstates = states.to_numpy()[:, :4]  # the i, E, g, J are in the 4 first columns
-    ndtrans = trans.to_numpy()
 
     iorig = np.array(ndstates[:, 0], dtype=int)
     maxii = int(np.max(iorig) + 1)
@@ -216,21 +257,13 @@ def pickup_gE(states, trans, dic_def, trans_lines=False):
     i_upper = np.array(ndtrans[:, 0], dtype=int)
     i_lower = np.array(ndtrans[:, 1], dtype=int)
 
-    # use the state counting numbers (upper and lower) as masking.
-    elower = newstates[i_lower, 0]
-    eupper = newstates[i_upper, 0]
-    gup = newstates[i_upper, 1]
-    jlower = np.array(newstates[i_lower, 2], dtype=int)
-    jupper = np.array(newstates[i_upper, 2], dtype=int)
-
     ### Step 2. Extra quantum numbers (e/f parity, vib and rot numbers)
     #!!! Todo: Awfull management of memory, i know! @minou
     n0 = (
         4 + dic_def["Landé"] + dic_def["lifetime"]
     )  # we need to shift if these quantities exist
-    ndstates_extra = states.to_numpy()[
-        :, n0:
-    ]  # the i, E, g, J are in the 4 first columns
+    # 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:
@@ -241,14 +274,51 @@ def pickup_gE(states, trans, dic_def, trans_lines=False):
                 ]
                 # 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
-
-    # See Issue #16
+    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."
+            )
+
+    # See Issue exojax#16
     # import matplotlib.pyplot as plt
     # nu_lines_t=ndtrans[:,3]
     # plt.plot(nu_lines_t-nu_lines,".",alpha=0.03)
@@ -257,7 +327,7 @@ def pickup_gE(states, trans, dic_def, trans_lines=False):
     # plt.savefig("nudiff.png", bbox_inches="tight", pad_inches=0.0)
     # plt.show()
 
-    return A, nu_lines, elower, gup, jlower, jupper, dico_quantumNumbers
+    return A, nu_lines, elower, gup, jlower, jupper, mask, dico_quantumNumbers
 
 
 # def pickup_gEslow(states, trans):
@@ -298,150 +368,138 @@ def pickup_gE(states, trans, dic_def, trans_lines=False):
 #     return A, nu_lines, elower, gup
 
 
-# def read_broad(broadf):
-#     """Reading braodening file (.broad)
-
-#     Args:
-#        broadf: .broad file
-
-#     Return:
-#        broadening info in bdat form (pandas), defined by this instance.
-
-#     Note:
-#        See Table 16 in https://arxiv.org/pdf/1603.05890.pdf
-
-#     """
-#     bdat = pd.read_csv(
-#         broadf,
-#         sep="\s+",
-#         names=(
-#             "code",
-#             "alpha_ref",
-#             "n_Texp",
-#             "jlower",
-#             "jupper",
-#             "kalower",
-#             "kclower",
-#             "kaupper",
-#             "kcupper",
-#             "v1lower",
-#             "v2lower",
-#             "v3lower",
-#             "v1upper",
-#             "v2upper",
-#             "v3upper",
-#         ),
-#     )
-#     return bdat
-
-
-# def check_bdat(bdat):
-#     """cheking codes in .broad
-
-#     Args:
-#        bdat: exomol .broad data given by exomolapi.read_broad
-
-#     Returns:
-#        code level: None, a0, a1, other codes unavailable currently,
-
-#     """
-
-#     def checkcode(code):
-#         cmask = bdat["code"] == code
-#         if len(bdat["code"][cmask]) > 0:
-#             return True
-#         else:
-#             return False
-
-#     codelv = None
-#     for code in ["a0", "a1"]:
-#         if checkcode(code):
-#             codelv = code
-
-#     return codelv
+def read_broad(broadf):
+    """Reading braodening file (.broad)
+    Args:
+       broadf: .broad file
 
+    Return:
+       broadening info in bdat form (pandas), defined by this instance.
+    Note:
+       See Table 16 in https://arxiv.org/pdf/1603.05890.pdf
+    """
+    bdat = pd.read_csv(
+        broadf,
+        sep="\s+",
+        names=(
+            "code",
+            "alpha_ref",
+            "n_Texp",
+            "jlower",
+            "jupper",
+            "kalower",
+            "kclower",
+            "kaupper",
+            "kcupper",
+            "v1lower",
+            "v2lower",
+            "v3lower",
+            "v1upper",
+            "v2upper",
+            "v3upper",
+        ),
+    )
+
+    return bdat
+
+
+def check_bdat(bdat):
+    """cheking codes in .broad
+    Args:
+       bdat: exomol .broad data given by exomolapi.read_broad
+    Returns:
+       code level: None, a0, a1, other codes unavailable currently,
+    """
 
-# def make_j2b(bdat, alpha_ref_default=0.07, n_Texp_default=0.5, jlower_max=None):
-#     """compute j2b (code a0, map from jlower to alpha_ref)
+    def checkcode(code):
+        cmask = bdat["code"] == code
+        if len(bdat["code"][cmask]) > 0:
+            return True
+        else:
+            return False
 
-#     Args:
-#        bdat: exomol .broad data given by exomolapi.read_broad
-#        alpha_ref_default: default value
-#        n_Texp_default: default value
-#        jlower_max: maximum number of jlower
+    codelv = None
+    for code in ["a0", "a1"]:
+        if checkcode(code):
+            codelv = code
 
-#     Returns:
-#        j2alpha_ref[jlower] provides alpha_ref for jlower
-#        j2n_Texp[jlower]  provides nT_exp for jlower
+    return codelv
 
-#     """
-#     # a0
-#     cmask = bdat["code"] == "a0"
-#     jlower_arr = np.array(bdat["jlower"][cmask], dtype=int)
-#     alpha_ref_arr = np.array(bdat["alpha_ref"][cmask])
-#     n_Texp_arr = np.array(bdat["n_Texp"][cmask])
-
-#     if jlower_max is None:
-#         Nblower = np.max(jlower_arr) + 1
-#     else:
-#         Nblower = np.max([jlower_max, np.max(jlower_arr)]) + 1
-#     j2alpha_ref = np.ones(Nblower) * alpha_ref_default
-#     j2n_Texp = np.ones(Nblower) * n_Texp_default
-
-#     j2alpha_ref[jlower_arr] = alpha_ref_arr
-#     j2n_Texp[jlower_arr] = n_Texp_arr
-
-#     Ndef = Nblower - (np.max(jlower_arr) + 1)
-#     if Ndef > 0:
-#         print(
-#             "default broadening parameters are used for ",
-#             Ndef,
-#             " J lower states in ",
-#             Nblower,
-#             " states",
-#         )
-
-#     return j2alpha_ref, j2n_Texp
-
-
-# def make_jj2b(bdat, j2alpha_ref_def, j2n_Texp_def, jupper_max=None):
-#     """compute jj2b (code a1, map from (jlower, jupper) to alpha_ref and n_Texp)
 
-#     Args:
-#        bdat: exomol .broad data given by exomolapi.read_broad
-#        j2alpha_ref_def: default value from a0
-#        j2n_Texp_def: default value from a0
-#        jupper_max: maximum number of jupper
+def make_j2b(bdat, alpha_ref_default=0.07, n_Texp_default=0.5, jlower_max=None):
+    """compute j2b (code a0, map from jlower to alpha_ref)
 
-#     Returns:
-#        jj2alpha_ref[jlower,jupper] provides alpha_ref for (jlower, jupper)
-#        jj2n_Texp[jlower,jupper]  provides nT_exp for (jlower, jupper)
+    Args:
+       bdat: exomol .broad data given by exomolapi.read_broad
+       alpha_ref_default: default value
+       n_Texp_default: default value
+       jlower_max: maximum number of jlower
+    Returns:
+       j2alpha_ref[jlower] provides alpha_ref for jlower
+       j2n_Texp[jlower]  provides nT_exp for jlower
+    """
+    # a0
+    cmask = bdat["code"] == "a0"
+    jlower_arr = np.array(bdat["jlower"][cmask], dtype=int)
+    alpha_ref_arr = np.array(bdat["alpha_ref"][cmask])
+    n_Texp_arr = np.array(bdat["n_Texp"][cmask])
+
+    if jlower_max is None:
+        Nblower = np.max(jlower_arr) + 1
+    else:
+        Nblower = np.max([jlower_max, np.max(jlower_arr)]) + 1
+    j2alpha_ref = np.ones(Nblower) * alpha_ref_default
+    j2n_Texp = np.ones(Nblower) * n_Texp_default
+
+    j2alpha_ref[jlower_arr] = alpha_ref_arr
+    j2n_Texp[jlower_arr] = n_Texp_arr
+
+    Ndef = Nblower - (np.max(jlower_arr) + 1)
+    if Ndef > 0:
+        print(
+            "default broadening parameters are used for ",
+            Ndef,
+            " J lower states in ",
+            Nblower,
+            " states",
+        )
 
-#     Note:
-#        The pair of (jlower, jupper) for which broadening parameters are not given, jj2XXX contains None.
+    return j2alpha_ref, j2n_Texp
 
-#     """
-#     # a1
-#     cmask = bdat["code"] == "a1"
-#     jlower_arr = np.array(bdat["jlower"][cmask], dtype=int)
-#     jupper_arr = np.array(bdat["jupper"][cmask], dtype=int)
-#     alpha_ref_arr = np.array(bdat["alpha_ref"][cmask])
-#     n_Texp_arr = np.array(bdat["n_Texp"][cmask])
 
-#     # Nblower = len(j2alpha_ref_def)
+def make_jj2b(bdat, j2alpha_ref_def, j2n_Texp_def, jupper_max=None):
+    """compute jj2b (code a1, map from (jlower, jupper) to alpha_ref and n_Texp)
 
-#     if jupper_max is None:
-#         Nbupper = np.max(jupper_arr) + 1
-#     else:
-#         Nbupper = np.max([jupper_max, np.max(jupper_arr)]) + 1
+    Args:
+       bdat: exomol .broad data given by exomolapi.read_broad
+       j2alpha_ref_def: default value from a0
+       j2n_Texp_def: default value from a0
+       jupper_max: maximum number of jupper
+    Returns:
+       jj2alpha_ref[jlower,jupper] provides alpha_ref for (jlower, jupper)
+       jj2n_Texp[jlower,jupper]  provides nT_exp for (jlower, jupper)
+    Note:
+       The pair of (jlower, jupper) for which broadening parameters are not given, jj2XXX contains None.
+    """
+    # a1
+    cmask = bdat["code"] == "a1"
+    jlower_arr = np.array(bdat["jlower"][cmask], dtype=int)
+    jupper_arr = np.array(bdat["jupper"][cmask], dtype=int)
+    alpha_ref_arr = np.array(bdat["alpha_ref"][cmask])
+    n_Texp_arr = np.array(bdat["n_Texp"][cmask])
+
+    if jupper_max is None:
+        Nbupper = np.max(jupper_arr) + 1
+    else:
+        Nbupper = np.max([jupper_max, np.max(jupper_arr)]) + 1
 
-#     jj2alpha_ref = j2alpha_ref_def[:, np.newaxis] * np.ones(Nbupper)
-#     jj2n_Texp = j2n_Texp_def[:, np.newaxis] * np.ones(Nbupper)
+    jj2alpha_ref = j2alpha_ref_def[:, np.newaxis] * np.ones(Nbupper)
+    jj2n_Texp = j2n_Texp_def[:, np.newaxis] * np.ones(Nbupper)
 
-#     jj2alpha_ref[jlower_arr, jupper_arr] = alpha_ref_arr
-#     jj2n_Texp[jlower_arr, jupper_arr] = n_Texp_arr
+    jj2alpha_ref[jlower_arr, jupper_arr] = alpha_ref_arr
+    jj2n_Texp[jlower_arr, jupper_arr] = n_Texp_arr
 
-#     return jj2alpha_ref, jj2n_Texp
+    return jj2alpha_ref, jj2n_Texp
 
 
 if __name__ == "__main__":
diff --git a/radis/io/hdf5.py b/radis/io/hdf5.py
index 9d456f82f..d66714f6e 100644
--- a/radis/io/hdf5.py
+++ b/radis/io/hdf5.py
@@ -7,7 +7,7 @@
 
 import os
 import sys
-from os.path import exists
+from os.path import exists, splitext
 from time import time
 
 import h5py
@@ -56,8 +56,42 @@ def update_pytables_to_vaex(fname, remove_initial=False, verbose=True, key="df")
 
 
 class HDF5Manager(object):
-    def __init__(self, engine="pytables"):
+    def __init__(self, engine=None):
+        """Class to handle all memory-mapping-librairies with one common API
+
+        All functions may not be fully implemetned, will raise a NotImplementedError
+        if that's not the case.
+
+        Librairies ::
+
+            'vaex'     > HDF5,  column-based
+            'pytables' > Pandas's HDF5,  row-based
+            'h5py'     > HDF5
+
+        Functions ::
+
+            add_metadata
+            read_metadata
+            write
+            load
+            guess_engine
+
+        Examples
+        --------
+        ::
+
+            file = 'CO.hdf5'
+            from radis.io.hdf5 import HDF5Manager
+            engine = HDF5Manager.guess_engine(file)
+            mgr = HDF5Manager(engine)
+            mgr.read_metadata(file)
+
+
+        """
         self.engine = engine
+        self._temp_batch_files = (
+            []
+        )  # list of batch files when writing by part in vaex mode
 
     def open(self, file, mode="w"):
         if self.engine == "pytables":
@@ -111,9 +145,23 @@ def write(
             import vaex
 
             if append == True:
-                raise ValueError(
-                    "Cannot append with 'vaex' engine. Load all files separately using vaex.open('many') then export to a single file"
-                )
+                # In vaex we cannot append. Here we write lots of small files then combine them.
+                # self.combine_temp_batch_files() should be combined at the end.
+                base, ext = splitext(file)
+                i = 0
+                temp_batch_file = base + "_temp" + str(i).zfill(6) + ext
+                while temp_batch_file in self._temp_batch_files:
+                    i += 1
+                    temp_batch_file = base + "_temp" + str(i).zfill(6) + ext
+                file = temp_batch_file
+                # Check no remaining one from a non-cleaned previous run:
+                if exists(file):
+                    from radis.misc.printer import printr
+
+                    printr(f"Temp file {temp_batch_file} already exists: deleting it")
+                    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
@@ -122,6 +170,39 @@ def write(
             raise NotImplementedError(self.engine)
             # h5py is not designed to write Pandas DataFrames
 
+    def combine_temp_batch_files(self, file, key="default"):
+        """Combine all batch files in ``self._temp_batch_files`` into one.
+        Removes all batch files.
+        """
+        if self.engine == "vaex":
+            if len(self._temp_batch_files) == 0:
+                raise ValueError(f"No batch temp files were written for {file}")
+            if key == "default":
+                key = r"/table"
+            import vaex
+
+            df = vaex.open(self._temp_batch_files, group=key)
+            df.export_hdf5(file, group=key, mode="w")
+            df.close()
+        self._close_temp_batch_files()
+
+    def _close_temp_batch_files(self):
+        for i in range(len(self._temp_batch_files) - 1, -1, -1):
+            os.remove(self._temp_batch_files[i])
+            del self._temp_batch_files[i]
+
+    def __del__(self):
+        """ clean before deleting"""
+        if len(self._temp_batch_files) > 0:
+            from radis.misc.printer import printr
+
+            printr(
+                "Warning : {0} files were written but not combined. Deleting them.\n\n{1}".format(
+                    len(self._temp_batch_files), self._temp_batch_files
+                )
+            )
+            self._close_temp_batch_files()
+
     def load(
         self,
         fname,
@@ -323,7 +404,11 @@ def read_metadata(self, fname: str, key="default") -> dict:
             if key == "default":
                 key = "df"
             with pd.HDFStore(fname, mode="r", complib="blosc", complevel=9) as f:
-                metadata = f.get_storer(key).attrs.metadata
+                try:
+                    metadata = f.get_storer(key).attrs.metadata
+                except KeyError as err:
+                    print(f"Error reading metadata from {fname}")
+                    raise err
 
         elif self.engine == "h5py":
             if key == "default":
@@ -332,7 +417,11 @@ def read_metadata(self, fname: str, key="default") -> dict:
                 if key is None:  # read metadta at root level
                     metadata = dict(hf.attrs)
                 else:
-                    metadata = dict(hf[key].attrs)
+                    try:
+                        metadata = dict(hf[key].attrs)
+                    except KeyError as err:
+                        print(f"Error reading metadata from {fname}")
+                        raise err
 
         elif self.engine == "vaex":
             if key == "default":
@@ -350,7 +439,11 @@ def read_metadata(self, fname: str, key="default") -> dict:
                     if key is None:  # add metadta at root level
                         metadata = dict(hf.attrs)
                     else:
-                        metadata = dict(hf[key].attrs)
+                        try:
+                            metadata = dict(hf[key].attrs)
+                        except KeyError as err:
+                            print(f"Error reading metadata from {fname}")
+                            raise err
 
         else:
             raise NotImplementedError(self.engine)
@@ -359,16 +452,40 @@ def read_metadata(self, fname: str, key="default") -> dict:
 
     @classmethod
     def guess_engine(self, file, verbose=True):
-        """Guess which HDF5 library ``file`` is compatible with"""
+        """Guess which HDF5 library ``file`` is compatible with
+
+        .. note::
+            it still take about 1 ms for this functino to execute. For extreme
+            performance you want to directly give the correct engine
+
+        Examples
+        --------
+        ::
+
+            file = 'CO.hdf5'
+            from radis.io.hdf5 import HDF5Manager
+            engine = HDF5Manager.guess_engine(file)
+            mgr = HDF5Manager(engine)
+            mgr.read_metadata(file)
+
+        """
         # See if it looks like PyTables
-        #  TODO : move in Manager
-        with pd.HDFStore(file, mode="r") as store:
-            if store.get_storer("df"):
-                engine = "pytables"
-            else:
-                engine = "h5py"
+        import tables
+
+        if tables.is_pytables_file(file):
+            engine = "pytables"
+        else:
+            # Try Vaex
+            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")
+        # raise
         return engine
 
 
@@ -380,7 +497,7 @@ def hdf2df(
     load_wavenum_max=None,
     verbose=True,
     store_kwargs={},
-    engine="auto",
+    engine="guess",
 ):
     """Load a HDF5 line databank into a Pandas DataFrame.
 
@@ -403,7 +520,7 @@ def hdf2df(
     store_kwargs: dict
         arguments forwarded to :py:meth:`~pandas.io.pytables.read_hdf`
     engine: ``'h5py'``, ``'pytables'``, ``'vaex'``, ``'auto'``
-        which HDF5 library to use. If ``'auto'``, try to guess. Note: ``'vaex'``
+        which HDF5 library to use. If ``'guess'``, try to guess. Note: ``'vaex'``
         uses ``'h5py'`` compatible HDF5.
 
     Returns
@@ -438,7 +555,7 @@ def hdf2df(
 
     """
 
-    if engine == "auto":
+    if engine == "guess":
         engine = HDF5Manager.guess_engine(fname)
 
     t0 = time()
@@ -457,13 +574,13 @@ def hdf2df(
             where.append(f'iso in {isotope.split(",")}')
 
     elif engine == "vaex":
-        #  Selection is done after opening the file time in vaex
-        where = []
+        # 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)
 
@@ -472,45 +589,28 @@ def hdf2df(
 
         # Selection
         selection = True
-        # limitation : so far the df.select(df.iso in isotope.split(",")) syntax
-        # fails in Vaex (and worse: still returns something.)
-        # For the moment, only implement with one isotope.
-        # TODO : add all cases manually...
+        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:
-            if not isinstance(isotope, int):
-                try:
-                    isotope = int(isotope)
-                except:
-                    raise NotImplementedError(
-                        f"When reading HDF5 in vaex mode, selection works for a single isotope only (got {isotope})"
-                    )
-        if (
-            load_wavenum_min is not None
-            and load_wavenum_max is not None
-            and isotope is not None
-        ):
-            df.select(
-                (df.wav > load_wavenum_min)
-                & (df.wav < load_wavenum_max)
-                & (df.iso == isotope)
-            )
-        elif load_wavenum_min is not None and load_wavenum_max is not None:
-            df.select((df.wav > load_wavenum_min) & (df.wav < load_wavenum_max))
-        elif load_wavenum_min is not None and isotope is not None:
-            df.select((df.wav > load_wavenum_min) & (df.iso == isotope))
-        elif load_wavenum_max is not None and isotope is not None:
-            df.select((df.wav < load_wavenum_max) & (df.iso == isotope))
-        elif load_wavenum_min is not None:
-            df.select(df.wav > load_wavenum_min)
-        elif load_wavenum_max is not None:
-            df.select(df.wav < load_wavenum_max)
-        elif isotope is not None:
-            df.select(df.iso == isotope)
-        else:
-            selection = False
+            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
 
         # Load
-        df = df.to_pandas_df(column_names=columns, selection=selection)
+        if columns:  # load only these columns (if they exist)
+            columns = [c for c in columns if c in df.columns]
+        df = df.to_pandas_df(column_names=columns)
 
     # Read and add metadata in the DataFrame
     metadata = manager.read_metadata(fname)
diff --git a/radis/io/hitemp.py b/radis/io/hitemp.py
index b1c2598a5..f12f64b0c 100644
--- a/radis/io/hitemp.py
+++ b/radis/io/hitemp.py
@@ -17,8 +17,6 @@
 
 import numpy as np
 
-import radis
-
 try:
     from .dbmanager import DatabaseManager
     from .hdf5 import update_pytables_to_vaex
@@ -107,12 +105,18 @@ def __init__(
         name,
         molecule,
         local_databases,
+        engine="default",
         verbose=True,
         chunksize=100000,
         parallel=True,
     ):
         super().__init__(
-            name, molecule, local_databases, verbose=verbose, parallel=parallel
+            name,
+            molecule,
+            local_databases,
+            engine,
+            verbose=verbose,
+            parallel=parallel,
         )
         self.chunksize = chunksize
         self.downloadable = True
@@ -269,7 +273,6 @@ def parse_to_local_file(
         pbar_Ntot_estimate_factor=None,
         pbar_Nlines_already=0,
         pbar_last=True,
-        engine="pytables",
     ):
         """Uncompress ``urlname`` into ``local_file``.
         Also add metadata
@@ -306,7 +309,7 @@ def parse_to_local_file(
         wmin = np.inf
         wmax = 0
 
-        writer = self.get_hdf5_manager(engine)
+        writer = self.get_hdf5_manager()
 
         with opener.open(urlname) as gfile:  # locally downloaded file
 
@@ -359,6 +362,7 @@ def parse_to_local_file(
                 b = np.zeros(
                     chunksize, dtype=dt
                 )  # receives the HITRAN 160-character data.
+        writer.combine_temp_batch_files(local_file)  # used for vaex mode only
         if pbar_last:
             pb.update(
                 Nlines_tot,
@@ -372,6 +376,8 @@ def parse_to_local_file(
         assert Nlines == Nlines_raw
 
         # Add metadata
+        from radis import __version__
+
         writer.add_metadata(
             local_file,
             {
@@ -380,18 +386,16 @@ def parse_to_local_file(
                 "download_date": self.get_today(),
                 "download_url": urlname,
                 "total_lines": Nlines_raw,
-                "version": radis.__version__,
+                "version": __version__,
             },
         )
 
         return Nlines
 
-    def register(self, engine):
-        """register in ~/radis.json
-
-        engine required to differentiate '.h5', '.hdf5'"""
+    def register(self):
+        """register in ~/radis.json"""
 
-        local_files, urlnames = self.get_filenames(engine=engine)
+        local_files, urlnames = self.get_filenames()
         info = f"HITEMP {self.molecule} lines ({self.wmin:.1f}-{self.wmax:.1f} cm-1) with TIPS-2017 (through HAPI) for partition functions"
 
         if self.molecule in ["CO2", "H2O"]:
@@ -431,7 +435,7 @@ def fetch_hitemp(
     chunksize=100000,
     clean_cache_files=True,
     return_local_path=False,
-    engine="pytables",
+    engine="default",
     parallel=True,
 ):
     """Stream HITEMP file from HITRAN website. Unzip and build a HDF5 file directly.
@@ -457,8 +461,9 @@ def fetch_hitemp(
 
     Other Parameters
     ----------------
-    cache: bool, or ``'regen'``
-        if ``True``, use existing HDF5 file. If ``False`` or ``'regen'``, rebuild it..
+    cache: ``True``, ``False``, ``'regen'`` or ``'force'``
+        if ``True``, use existing HDF5 file. If ``False`` or ``'regen'``, rebuild it.
+        If ``'force'``, raise an error if cache file cannot be used (useful for debugging).
         Default ``True``.
     verbose: bool
     chunksize: int
@@ -468,8 +473,8 @@ def fetch_hitemp(
         if ``True`` clean downloaded cache files after HDF5 are created.
     return_local_path: bool
         if ``True``, also returns the path of the local database file.
-    engine: 'pytables', 'vaex'
-        which HDF5 library to use.
+    engine: 'pytables', 'vaex', 'default'
+        which HDF5 library to use. If 'default' use the value from ~/radis.json
     parallel: bool
         if ``True``, uses joblib.parallel to load database with multiple processes
 
@@ -511,10 +516,10 @@ def fetch_hitemp(
 
     """
 
-    if databank_name == "HITEMP-{molecule}":
+    if r"{molecule}" in databank_name:
         databank_name = databank_name.format(**{"molecule": molecule})
 
-    local_databases = abspath(local_databases.replace("~", expanduser("~")))
+    local_databases = abspath(expanduser(local_databases))
 
     ldb = HITEMPDatabaseManager(
         databank_name,
@@ -523,10 +528,11 @@ def fetch_hitemp(
         verbose=verbose,
         chunksize=chunksize,
         parallel=parallel,
+        engine=engine,
     )
 
     # Get list of all expected local files for this database:
-    local_files, urlnames = ldb.get_filenames(engine=engine)
+    local_files, urlnames = ldb.get_filenames()
 
     # Delete files if needed:
     relevant_files = ldb.keep_only_relevant(
@@ -536,8 +542,7 @@ def fetch_hitemp(
         ldb.remove_local_files(relevant_files)
     ldb.check_deprecated_files(
         ldb.get_existing_files(relevant_files),
-        engine=engine,
-        remove=True if cache != "force" else False,
+        auto_remove=True if cache != "force" else False,
     )
 
     # Get missing files
@@ -546,11 +551,11 @@ def fetch_hitemp(
         download_files, load_wavenum_min, load_wavenum_max
     )
     # do not re-download files if they exist in another format :
-    if engine in ["vaex", "auto"]:
+    if engine in ["vaex", "auto", "default"]:
         # ... convert files if asked:
-        import radis
+        from radis import config
 
-        if radis.config["AUTO_UPDATE_DATABASE"]:
+        if config["AUTO_UPDATE_DATABASE"]:
             converted = []
             for f in download_files:
                 if exists(f.replace(".hdf5", ".h5")):
@@ -569,7 +574,7 @@ def fetch_hitemp(
             urlnames = ldb.fetch_urlnames()
         filesmap = dict(zip(local_files, urlnames))
         download_urls = [filesmap[k] for k in download_files]
-        ldb.download_and_parse(download_urls, download_files, engine=engine)
+        ldb.download_and_parse(download_urls, download_files)
 
         # Done: add final checks
         if molecule not in ["CO2", "H2O"]:
@@ -587,9 +592,9 @@ def fetch_hitemp(
 
     # Register
     if not ldb.is_registered():
-        ldb.register(engine=engine)  # engine required to differentiate '.h5', '.hdf5'
+        ldb.register()
 
-    if clean_cache_files:
+    if len(download_files) > 0 and clean_cache_files:
         ldb.clean_download_files()
 
     # Load and return
@@ -606,7 +611,6 @@ def fetch_hitemp(
         isotope=isotope,
         load_wavenum_min=load_wavenum_min,  # for relevant files, get only the right range
         load_wavenum_max=load_wavenum_max,
-        engine=engine,
     )
 
     return (df, files_loaded) if return_local_path else df
diff --git a/radis/io/hitran.py b/radis/io/hitran.py
index cb012c9be..4ef161861 100644
--- a/radis/io/hitran.py
+++ b/radis/io/hitran.py
@@ -20,17 +20,17 @@
 """
 
 
+import os
 import sys
 
 # from radis.test.utils import getTestFile
 import time
 from collections import OrderedDict
-from os.path import exists, getmtime
+from os.path import exists, getmtime, join, split
 
 import pandas as pd
 from numpy import int64
 
-import radis
 from radis.db.classes import (  # get_molecule_identifier,
     HITRAN_CLASS1,
     HITRAN_CLASS2,
@@ -51,6 +51,7 @@
     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,
@@ -99,6 +100,7 @@ def cast_to_int64_with_missing_values(dg, keys):
 def hit2df(
     fname,
     cache=True,
+    load_columns=None,
     verbose=True,
     drop_non_numeric=True,
     load_wavenum_min=None,
@@ -117,6 +119,12 @@ def hit2df(
         improves performances a lot (but changes in the database are not
         taken into account). If False, no database is used. If ``'regen'``, temp
         file are reconstructed. Default ``True``.
+    load_columns: list
+        columns to load. If ``None``, loads everything
+
+        .. note::
+            this is only relevant if loading from a cache file. To generate
+            the cache file, all columns are loaded anyway.
 
     Other Parameters
     ----------------
@@ -144,16 +152,8 @@ def hit2df(
     .. [1] `HITRAN 1996, Rothman et al., 1998 <https://www.sciencedirect.com/science/article/pii/S0022407398000788>`__
 
 
-
-    Notes
-    -----
-
-    Performances: see CDSD-HITEMP parser
-
-
     See Also
     --------
-
     :func:`~radis.io.cdsd.cdsd2df`
     """
     metadata = {}
@@ -176,14 +176,17 @@ def hit2df(
         relevant_if_metadata_below = (
             {"wavenum_min": load_wavenum_max} if load_wavenum_max else {}
         )  # not relevant if wavenum_min of file is > wavenum max required
+        from radis import __version__, config
+
         df = load_h5_cache_file(
             fcache,
             cache,
+            columns=load_columns,
             valid_if_metadata_is=metadata,
             relevant_if_metadata_above=relevant_if_metadata_above,
             relevant_if_metadata_below=relevant_if_metadata_below,
-            current_version=radis.__version__,
-            last_compatible_version=radis.config["OLDEST_COMPATIBLE_VERSION"],
+            current_version=__version__,
+            last_compatible_version=config["OLDEST_COMPATIBLE_VERSION"],
             verbose=verbose,
             engine=engine,
         )
@@ -274,13 +277,14 @@ def hit2df(
                     fcache, new_metadata
                 )
             )
+        from radis import __version__
+
         try:
             save_to_hdf(
                 df,
                 fcache,
                 metadata=new_metadata,
-                version=radis.__version__,
-                key="df",
+                version=__version__,
                 overwrite=True,
                 verbose=verbose,
                 engine=engine,
@@ -1036,6 +1040,340 @@ def parse_global_quanta(df, mol, verbose=True):
     return df
 
 
+#%%
+
+
+class HITRANDatabaseManager(DatabaseManager):
+    def __init__(
+        self,
+        name,
+        molecule,
+        local_databases,
+        engine="default",
+        verbose=True,
+        parallel=True,
+    ):
+        super().__init__(
+            name,
+            molecule,
+            local_databases,
+            engine=engine,
+            verbose=verbose,
+            parallel=parallel,
+        )
+        self.downloadable = True
+        self.base_url = None
+        self.Nlines = None
+
+    def get_filenames(self):
+        if self.engine == "vaex":
+            return [join(self.local_databases, f"{self.molecule}.hdf5")]
+        else:
+            raise NotImplementedError()
+
+    def download_and_parse(
+        self,
+        local_file,
+        cache=True,
+    ):
+        """Download from HITRAN and parse into ``local_file``.
+        Also add metadata
+
+        Overwrites :py:meth:`radis.io.dbmanager.DatabaseManager.download_and_parse`
+        which downloads from a list of URL, because here we use [HAPI]_ to
+        download the files.
+
+        Parameters
+        ----------
+        opener: an opener with an .open() command
+        gfile : file handler. Filename: for info"""
+
+        from hapi import db_begin, fetch
+
+        from radis import hit2df
+        from radis.db.classes import get_molecule_identifier
+
+        if isinstance(local_file, list):
+            assert (
+                len(local_file) == 1
+            )  # fetch_hitran stores all lines of a given molecule in one file
+            local_file = local_file[0]
+
+        wmin = 1
+        wmax = 40000
+
+        def download_all_hitran_isotopes(molecule, directory):
+            """Blindly try to download all isotpes 1 - 9 for the given molecule
+
+            .. warning::
+                this won't be able to download higher isotopes (ex : isotope 10-11-12 for CO2)
+                Neglected for the moment, they're irrelevant for most calculations anyway
+
+            """
+            # create temp folder :
+            from radis.misc.basics import make_folders
+
+            make_folders(*split(abspath(directory)))
+
+            db_begin(directory)
+            isotope_list = []
+            data_file_list = []
+            header_file_list = []
+            for iso in range(1, 10):
+                file = f"{molecule}_{iso}"
+                if exists(join(directory, file + ".data")):
+                    if cache == "regen":
+                        # remove without printing message
+                        os.remove(join(directory, file + ".data"))
+                    else:
+                        from radis.misc.printer import printr
+
+                        printr(
+                            "File already exist: {0}. Deleting it.`".format(
+                                join(directory, file + ".data")
+                            )
+                        )
+                        os.remove(join(directory, file + ".data"))
+                try:
+                    fetch(file, get_molecule_identifier(molecule), iso, wmin, wmax)
+                except KeyError:
+                    # Isotope not defined:
+                    continue
+                else:
+                    isotope_list.append(iso)
+                    data_file_list.append(file + ".data")
+                    header_file_list.append(file + ".header")
+            return isotope_list, data_file_list, header_file_list
+
+        molecule = self.molecule
+        wmin_final = 100000
+        wmax_final = -1
+
+        # make_folders(dirname(self.local_databases), basename(self.local_databases))
+        # make_folders(self.local_databases, ["downloads"])
+
+        # Use HAPI only to download the files, then we'll parse them with RADIS's
+        # parsers, and convert to RADIS's fast HDF5 file formats.
+        isotope_list, data_file_list, header_file_list = download_all_hitran_isotopes(
+            molecule, self.tempdir
+        )
+
+        writer = self.get_hdf5_manager()
+
+        # Create HDF5 cache file for all isotopes
+        Nlines = 0
+        for iso, data_file in zip(isotope_list, data_file_list):
+            df = hit2df(
+                join(self.tempdir, data_file),
+                cache=False,  # do not generate cache yet
+            )
+            # if engine == 'vaex':
+            #     df.executor.async_method = "awaitio"   # Temp fix for https://github.com/spyder-ide/spyder/issues/16183
+
+            wmin_final = min(wmin_final, df.wav.min())
+            wmax_final = max(wmax_final, df.wav.max())
+            Nlines += len(df)
+
+            writer.write(
+                local_file, df, append=True
+            )  # create temporary files if required
+
+        # Open all HDF5 cache files and export in a single file with Vaex
+        writer.combine_temp_batch_files(local_file)  # used for vaex mode only
+
+        self.wmin = wmin_final
+        self.wmax = wmax_final
+
+        # Add metadata
+        from radis import __version__
+
+        writer.add_metadata(
+            local_file,
+            {
+                "wavenumber_min": self.wmin,
+                "wavenumber_max": self.wmax,
+                "download_date": self.get_today(),
+                "download_url": "downloaded by HAPI, parsed & store with RADIS",
+                "total_lines": Nlines,
+                "version": __version__,
+            },
+        )
+
+        # # clean downloaded files  TODO
+        # for file in data_file_list + header_file_list:
+        #     os.remove(join(self.local_databases, "downloads", file))
+
+    def register(self):
+        """register in ~/radis.json"""
+
+        from radis.db import MOLECULES_LIST_NONEQUILIBRIUM
+
+        local_files = self.get_filenames()
+        info = f"HITRAN {self.molecule} lines ({self.wmin:.1f}-{self.wmax:.1f} cm-1) with TIPS-2021 (through HAPI) for partition functions"
+
+        dict_entries = {
+            "info": info,
+            "path": local_files,
+            "format": "hdf5-radisdb",
+            "parfuncfmt": "hapi",
+            "wavenumber_min": self.wmin,
+            "wavenumber_max": self.wmax,
+            "download_date": self.get_today(),
+        }
+
+        # Add energy level calculation
+        if self.molecule in MOLECULES_LIST_NONEQUILIBRIUM:
+            dict_entries[
+                "info"
+            ] += " and RADIS spectroscopic constants for rovibrational energies (nonequilibrium)"
+            dict_entries["levelsfmt"] = "radis"
+
+        super().register(dict_entries)
+
+
+from os.path import abspath, expanduser
+
+
+# TODO: implement parallel=True for all isotopes ?
+def fetch_hitran(
+    molecule,
+    local_databases="~/.radisdb/hitran/",
+    databank_name="HITRAN-{molecule}",
+    isotope=None,
+    load_wavenum_min=None,
+    load_wavenum_max=None,
+    columns=None,
+    cache=True,
+    verbose=True,
+    clean_cache_files=True,
+    return_local_path=False,
+    engine="default",
+    parallel=True,
+):
+    """Download all HITRAN lines from HITRAN website. Unzip and build a HDF5 file directly.
+
+    Returns a Pandas DataFrame containing all lines.
+
+    Parameters
+    ----------
+    molecule: `"H2O", "CO2", "N2O", "CO", "CH4", "NO", "NO2", "OH"`
+        HITEMP molecule. See https://hitran.org/hitemp/
+    local_databases: str
+        where to create the RADIS HDF5 files. Default ``"~/.radisdb/"``
+    databank_name: str
+        name of the databank in RADIS :ref:`Configuration file <label_lbl_config_file>`
+        Default ``"HITRAN-{molecule}"``
+    isotope: str
+        load only certain isotopes : ``'2'``, ``'1,2'``, etc. If ``None``, loads
+        everything. Default ``None``.
+    load_wavenum_min, load_wavenum_max: float (cm-1)
+        load only specific wavenumbers.
+    columns: list of str
+        list of columns to load. If ``None``, returns all columns in the file.
+
+    Other Parameters
+    ----------------
+    cache: ``True``, ``False``, ``'regen'`` or ``'force'``
+        if ``True``, use existing HDF5 file. If ``False`` or ``'regen'``, rebuild it.
+        If ``'force'``, raise an error if cache file cannot be used (useful for debugging).
+        Default ``True``.
+    verbose: bool
+    clean_cache_files: bool
+        if ``True`` clean downloaded cache files after HDF5 are created.
+    return_local_path: bool
+        if ``True``, also returns the path of the local database file.
+    engine: 'pytables', 'vaex', 'default'
+        which HDF5 library to use. If 'default' use the value from ~/radis.json
+    parallel: bool
+        if ``True``, uses joblib.parallel to load database with multiple processes
+
+    Returns
+    -------
+    df: pd.DataFrame
+        Line list
+        A HDF5 file is also created in ``local_databases`` and referenced
+        in the :ref:`RADIS config file <label_lbl_config_file>` with name
+        ``databank_name``
+    local_path: str
+        path of local database file if ``return_local_path``
+
+    Examples
+    --------
+    ::
+
+        from radis import fetch_hitemp
+        df = fetch_hitemp("CO")
+        print(df.columns)
+        >>> Index(['id', 'iso', 'wav', 'int', 'A', 'airbrd', 'selbrd', 'El', 'Tdpair',
+            'Pshft', 'ierr', 'iref', 'lmix', 'gp', 'gpp', 'Fu', 'branch', 'jl',
+            'syml', 'Fl', 'vu', 'vl'],
+            dtype='object')
+
+    .. minigallery:: radis.io.hitemp.fetch_hitemp
+
+    Notes
+    -----
+    if using ``load_only_wavenum_above/below`` or ``isotope``, the whole
+    database is anyway downloaded and uncompressed to ``local_databases``
+    fast access .HDF5 files (which will take a long time on first call). Only
+    the expected wavenumber range & isotopes are returned. The .HFD5 parsing uses
+    :py:func:`~radis.io.hdf5.hdf2df`
+
+    See Also
+    --------
+    :py:func:`~radis.io.hdf5.hdf2df`, :py:meth:`~radis.lbl.loader.DatabankLoader.fetch_databank`
+
+    """
+
+    if r"{molecule}" in databank_name:
+        databank_name = databank_name.format(**{"molecule": molecule})
+
+    local_databases = abspath(local_databases.replace("~", expanduser("~")))
+
+    ldb = HITRANDatabaseManager(
+        databank_name,
+        molecule=molecule,
+        local_databases=local_databases,
+        engine=engine,
+        verbose=verbose,
+        parallel=parallel,
+    )
+
+    # Get expected local files for this database:
+    local_file = ldb.get_filenames()
+
+    # Delete files if needed:
+    if cache == "regen":
+        ldb.remove_local_files(local_file)
+    ldb.check_deprecated_files(
+        ldb.get_existing_files(local_file),
+        auto_remove=True if cache != "force" else False,
+    )
+
+    # Download files
+    download_files = ldb.get_missing_files(local_file)
+    if download_files:
+        ldb.download_and_parse(download_files, cache=cache)
+
+    # Register
+    if not ldb.is_registered():
+        ldb.register()
+
+    if len(download_files) > 0 and clean_cache_files:
+        ldb.clean_download_files()
+
+    # Load and return
+    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,
+    )
+
+    return (df, local_file) if return_local_path else df
+
+
 # ======================================================
 # %% Test
 
diff --git a/radis/io/query.py b/radis/io/query.py
index 04c3e212a..c9e90d804 100644
--- a/radis/io/query.py
+++ b/radis/io/query.py
@@ -24,9 +24,12 @@
 import numpy as np
 import pandas as pd
 
-import radis
 from radis.db.classes import get_molecule, get_molecule_identifier
-from radis.io.cache_files import check_cache_file, get_cache_file, save_to_hdf
+
+try:
+    from .cache_files import check_cache_file, get_cache_file, save_to_hdf
+except ImportError:
+    from radis.io.cache_files import check_cache_file, get_cache_file, save_to_hdf
 from radis.misc import is_float
 from radis.misc.printer import printr
 
@@ -34,7 +37,14 @@
 
 
 def fetch_astroquery(
-    molecule, isotope, wmin, wmax, verbose=True, cache=True, expected_metadata={}
+    molecule,
+    isotope,
+    wmin,
+    wmax,
+    verbose=True,
+    cache=True,
+    expected_metadata={},
+    engine="pytables-fixed",
 ):
     """Download a HITRAN line database to a Pandas DataFrame.
 
@@ -116,13 +126,11 @@ def fetch_astroquery(
                 expected_metadata=expected_metadata,
                 compare_as_close=["wmin", "wmax"],
                 verbose=verbose,
-                engine="pytables-fixed",
+                engine=engine,
             )
             if exists(fcache):
                 try:
-                    return get_cache_file(
-                        fcache, verbose=verbose, engine="pytables-fixed"
-                    )
+                    return get_cache_file(fcache, verbose=verbose, engine=engine)
                 except Exception as err:
                     if verbose:
                         printr(
@@ -247,16 +255,18 @@ def fetch_astroquery(
                     fcache, new_metadata
                 )
             )
+        from radis import __version__
+
         try:
             save_to_hdf(
                 df,
                 fcache,
                 metadata=new_metadata,
-                version=radis.__version__,
+                version=__version__,
                 key="df",
                 overwrite=True,
                 verbose=verbose,
-                engine="pytables-fixed",
+                engine=engine,
             )
         except PermissionError:
             if verbose:
diff --git a/radis/lbl/base.py b/radis/lbl/base.py
index d444c1f8e..b7ac44071 100644
--- a/radis/lbl/base.py
+++ b/radis/lbl/base.py
@@ -1378,6 +1378,12 @@ def _add_ju(self, df):
         P, Q, R are replaced by -1, 0, 1 in the DataFrame to ensure that all
         terms are numeric (improves performances)
         """
+        if "branch" not in df:
+            raise KeyError(
+                f"`branch` not defined in database columns: ({list(df.columns)}). "
+                + "You can add it with `load_columns=['branch', ...]` or simply `load_columns='noneq'` / `load_columns='all'` in fetch_databank / load_databank()"
+            )
+
         if df.dtypes["branch"] != np.int64:
             raise DeprecationWarning(
                 "For performance purpose, numeric (-1, 0, 1) "
@@ -1443,8 +1449,9 @@ def _calc_noneq_parameters(self, vib_distribution, singleTvibmode):
                 if "globu" in df:
                     error_details = ". However, it looks like `globu` is defined. Maybe HITRAN-like database wasn't fully parsed? See radis.io.hitran.hit2df"
                 raise KeyError(
-                    f"`{k}` not defined in database ({list(df.columns)})"
+                    f"`{k}` not defined in database ({list(df.columns)}). "
                     + error_details
+                    + f"Make sure you properly load parameters required for non-LTE calculations by adding `load_columns=['{k}', ...]` or simply `load_columns='noneq'` in fetch_databank / load_databank()"
                 )
 
         # Make sure database has pre-computed non equilibrium quantities
@@ -1687,9 +1694,10 @@ def get_molar_mass(self, df):
                 ]
             except KeyError:
                 raise NotImplementedError(
-                    "Molar mass of {0} (isotope {1}) is unknown. You can manually add it in your radis.json file {'molparams':{'molar_mass':{'molecule':{'ISOTOPE':...}}}}; or in the SpectrumFactory._EXTRA_MOLAR_MASS[molecule][isotope] = M (g/mol) dictionary. Please also report on GitHub so we can update !".format(
+                    "Molar mass of {0} (isotope {1}) is unknown.".format(
                         df.attrs["molecule"], df.attrs["iso"]
                     )
+                    + " You can manually add it in your radis.json file {'molparams':{'molar_mass':{'molecule':{'ISOTOPE':...}}}}; or in the SpectrumFactory._EXTRA_MOLAR_MASS[molecule][isotope] = M (g/mol) dictionary. Please also report on GitHub so we can update !"
                 )
 
         if "iso" in df.columns:
diff --git a/radis/lbl/calc.py b/radis/lbl/calc.py
index ebad7593e..eff277ce1 100644
--- a/radis/lbl/calc.py
+++ b/radis/lbl/calc.py
@@ -134,16 +134,33 @@ def calc_spectrum(
 
     databank: str or dict
         can be either:
-        - ``'hitran'``, to fetch automatically the latest HITRAN version
-          through :py:func:`~radis.io.query.fetch_astroquery` (based on [HAPI]_ ).
-        - ``'hitemp'``, to fetch automatically the latest HITEMP version
-          through :py:func:`~radis.io.hitemp.fetch_hitemp`.
+        - ``'hitran'``, to fetch the latest HITRAN version
+          through :py:func:`~radis.io.hitran.fetch_hitran` (download full database
+          with  [HAPI]_) or :py:func:`~radis.io.query.fetch_astroquery` (download
+          only the required range). To use one mode or the other, use ::
+
+            databank=('hitran', 'full')     # download and cache full database, all isotopes
+            databank=('hitran', 'range')    # download and cache required range, required isoope
+
+        - ``'hitemp'``, to fetch the latest HITEMP version
+          through :py:func:`~radis.io.hitemp.fetch_hitemp`. Downloads all lines
+          and all isotopes.
+        - ``'exomol'``, to fetch the latest ExoMol database
+          through :py:func:`~radis.io.hitemp.fetch_exomol`. To download a specific
+          database use (more info in fetch_exomol) ::
+
+            databank=('exomol', 'EBJT')   # 'EBJT' is a specific ExoMol database name
+
         - the name of a a valid database file, in which case the format is inferred.
           For instance, ``'.par'`` is recognized as ``hitran/hitemp`` format.
-          Accepts wildcards ``'*'`` to select multiple files.
+          Accepts wildcards ``'*'`` to select multiple files ::
+
+            databank='PATH/TO/co_*.par'
+
         - the name of a spectral database registered in your ``~/radis.json``
-          :ref:`configuration file <label_lbl_config_file>`.
-          This allows to use multiple database files.
+          :ref:`configuration file <label_lbl_config_file>` ::
+
+            databank='MY_SPECTRAL_DATABASE'
 
         Default ``'hitran'``. See :class:`~radis.lbl.loader.DatabankLoader` for more
         information on line databases, and :data:`~radis.misc.config.DBFORMAT` for
@@ -633,8 +650,16 @@ def _is_at_equilibrium():
         export_lines=export_lines,
         **kwargs
     )
-    if databank in ["fetch", "hitran", "hitemp", "exomol",] or (
-        isinstance(databank, tuple) and databank[0] == "exomol"
+    if (
+        databank
+        in [
+            "fetch",
+            "hitran",
+            "hitemp",
+            "exomol",
+        ]
+        or (isinstance(databank, tuple) and databank[0] == "exomol")
+        or (isinstance(databank, tuple) and databank[0] == "hitran")
     ):  # mode to get databank without relying on  Line databases
         # Line database :
         if databank in ["fetch", "hitran"]:
@@ -642,6 +667,14 @@ def _is_at_equilibrium():
                 "source": "hitran",
                 "parfuncfmt": "hapi",  # use HAPI (TIPS) partition functions for equilibrium
             }
+        elif isinstance(databank, tuple) and databank[0] == "hitran":
+            conditions = {
+                "source": "hitran",
+                "database": databank[
+                    1
+                ],  # 'full' or 'partial', cf LoaderFactory.fetch_databank()
+                "parfuncfmt": "hapi",  # use HAPI (TIPS) partition functions for equilibrium
+            }
         elif databank in ["hitemp"]:
             conditions = {
                 "source": "hitemp",
@@ -655,7 +688,7 @@ def _is_at_equilibrium():
         elif isinstance(databank, tuple) and databank[0] == "exomol":
             conditions = {
                 "source": "exomol",
-                "exomol_database": databank[1],
+                "database": databank[1],
                 "parfuncfmt": "exomol",  # download & use Exo partition functions for equilibrium}
             }
         # Partition functions :
@@ -671,9 +704,19 @@ def _is_at_equilibrium():
             # constants (not all molecules are supported!)
             conditions["levelsfmt"] = "radis"
             conditions["lvl_use_cached"] = use_cached
+
+        # Columns to load
+        if export_lines:
+            conditions["load_columns"] = "all"
+        elif not _equilibrium:
+            conditions["load_columns"] = "noneq"
+        else:
+            conditions["load_columns"] = "equilibrium"
+
         conditions["load_energies"] = not _equilibrium
         # Details to identify lines
         conditions["parse_local_global_quanta"] = (not _equilibrium) or export_lines
+        # Finally, LOAD :
         sf.fetch_databank(**conditions)
     elif exists(databank):
         conditions = {
@@ -713,14 +756,34 @@ def _is_at_equilibrium():
                 + "and define the format there. More information on "
                 + "https://radis.readthedocs.io/en/latest/lbl/lbl.html#configuration-file"
             )
+
+        # Columns to load
+        if export_lines:
+            conditions["load_columns"] = "all"
+        elif not _equilibrium:
+            conditions["load_columns"] = "noneq"
+        else:
+            conditions["load_columns"] = "equilibrium"
+
         conditions["load_energies"] = not _equilibrium
+        # Finally, LOAD :
         sf.load_databank(**conditions)
 
     else:  # manual mode: get from user-defined line databases defined in ~/radis.json
+
+        # Columns to load
+        if export_lines:
+            load_columns = "all"
+        elif not _equilibrium:
+            load_columns = "noneq"
+        else:
+            load_columns = "equilibrium"
+
         sf.load_databank(
             databank,
             load_energies=not _equilibrium,  # no need to load/calculate energies at eq.
             drop_columns=drop_columns,
+            load_columns=load_columns,
         )
 
     #    # Get optimisation strategies
diff --git a/radis/lbl/loader.py b/radis/lbl/loader.py
index 041072cdd..eac33bca3 100644
--- a/radis/lbl/loader.py
+++ b/radis/lbl/loader.py
@@ -58,6 +58,7 @@
 import numpy as np
 import pandas as pd
 
+from radis import config
 from radis.db.classes import get_molecule
 from radis.db.molecules import getMolecule
 from radis.db.molparam import MolParams
@@ -67,7 +68,12 @@
 from radis.io.exomol import fetch_exomol
 from radis.io.hdf5 import hdf2df
 from radis.io.hitemp import fetch_hitemp
-from radis.io.hitran import hit2df, parse_global_quanta, parse_local_quanta
+from radis.io.hitran import (
+    fetch_hitran,
+    hit2df,
+    parse_global_quanta,
+    parse_local_quanta,
+)
 from radis.io.query import fetch_astroquery
 from radis.io.tools import drop_object_format_columns, replace_PQR_with_m101
 from radis.levels.partfunc import (
@@ -199,7 +205,7 @@
     "Pshft",
     "El",
 ]
-""" dict: drop all columns but these if using ``drop_columns='all'`` in load_databank
+""" list: drop all columns but these if using ``drop_columns='all'`` in load_databank
 Note: nonequilibrium calculations wont be possible anymore and it wont be possible
 to identify lines with :py:meth:`~radis.spectrum.spectrum.Spectrum.line_survey`
 
@@ -209,6 +215,31 @@
 - 'cdsd-hitemp' (CDSD HITEMP): :data:`~radis.io.cdsd.columns_hitemp`,
 - 'cdsd-4000': (CDSD 4000) :data:`~radis.io.cdsd.columns_4000`,
 """
+required_non_eq = [
+    "branch",
+    "jl",
+    "vl",
+    "vu",
+    "v1l",
+    "v2l",
+    "l2l",
+    "v3l",
+    "v1u",
+    "v2u",
+    "l2u",
+    "v3u",
+    "polyl",
+    "wangl",
+    "rankl",
+    "polyu",
+    "wangu",
+    "ranku",
+]
+"""list: column names required for non-equilibrium calculations.
+See load_column= key of fetch_databank() and load_databank() """
+# TODO refactor : directly go & parse the identifications names (globu, locu, etc.)
+# in radis.io.hitran ?
+# For the moment we just try to be exhaustive
 
 # @dev: Sanity checks
 # (make sure all variables are defined everywhere)
@@ -504,7 +535,7 @@ class MiscParams(ConditionDict):
         "warning_linestrength_cutoff",
         "total_lines",
         "zero_padding",
-        "hdf5_engine",
+        "memory_mapping_engine",
         "add_at_used",  # function used in DIT ; a Cython and a pure-Python version exist
     ]
 
@@ -525,7 +556,9 @@ def __init__(self):
         )
         self.warning_linestrength_cutoff = None  #: float [0-1]: raise a warning if the sum of linestrength cut is above that
         self.total_lines = 0  #: int : number of lines in database.
-        self.hdf5_engine = "pytables"  # 'pytables', 'vaex' (/!\ experimental in 0.9.30)
+        self.memory_mapping_engine = config[
+            "MEMORY_MAPPING_ENGINE"
+        ]  # 'pytables', 'vaex', 'feather'
 
         self.add_at_used = (
             ""  # function used in DIT ; a Cython and a pure-Python version exist
@@ -686,7 +719,6 @@ def __init__(self):
         # HARDCODED molar mass; for WIP ExoMol implementation, until MolParams
         # is an attribute and can be updated with definitions from ExoMol.
         # https://github.com/radis/radis/issues/321
-        from radis import config
 
         self._EXTRA_MOLAR_MASS = config["molparams"]["molar_mass"]
         """Extra molar mass when not found in HITRAN molecular parameter database
@@ -793,7 +825,6 @@ def init_databank(self, *args, **kwargs):
             ``True``, includes off-range, neighbouring lines that contribute
             because of lineshape broadening. The ``neighbour_lines``
             parameter is used to determine the limit. Default ``True``.
-        *Other arguments are related to how to open the files*
         drop_columns: list
             columns names to drop from Line DataFrame after loading the file.
             Not recommended to use, unless you explicitely want to drop information
@@ -802,6 +833,21 @@ def init_databank(self, *args, **kwargs):
             are dropped. See :data:`~radis.lbl.loader.drop_auto_columns_for_dbformat`
             and :data:`~radis.lbl.loader.drop_auto_columns_for_levelsfmt`.
             Default ``'auto'``.
+        load_columns: list, ``'all'``, ``'equilibrium'``, ``'noneq'``
+            columns names to load.
+            If ``'equilibrium'``, only load the columns required for equilibrium
+            calculations. If ``'noneq'``, also load the columns required for
+            non-LTE calculations. See :data:`~radis.lbl.loader.drop_all_but_these`.
+            If ``'all'``, load everything. Note that for performances, it is
+            better to load only certain columsn rather than loading them all
+            and dropping them with ``drop_columns``.
+            Default ``'equilibrium'``.
+
+            .. warning::
+                if using ``'equilibrium'``, not all parameters will be available
+                for a Spectrum :py:func:`~radis.spectrum.spectrum.Spectrum.line_survey`.
+
+        *Other arguments are related to how to open the files*
 
         Notes
         -----
@@ -816,6 +862,7 @@ def init_databank(self, *args, **kwargs):
         - Load from local files: :meth:`~radis.lbl.loader.DatabankLoader.load_databank`
         - Reload databank: :meth:`~radis.lbl.loader.DatabankLoader._check_line_databank`
         """
+        # TODO : refactor drop_columns/load_columns
 
         # Check inputs
         (
@@ -829,6 +876,7 @@ def init_databank(self, *args, **kwargs):
             db_use_cached,
             lvl_use_cached,
             drop_columns,
+            load_columns,
             load_energies,
             include_neighbouring_lines,
         ) = self._check_database_params(*args, **kwargs)
@@ -861,7 +909,7 @@ def init_databank(self, *args, **kwargs):
     def fetch_databank(
         self,
         source="hitran",
-        exomol_database=None,
+        database="default",
         parfunc=None,
         parfuncfmt="hapi",
         levels=None,
@@ -872,7 +920,8 @@ def fetch_databank(
         drop_non_numeric=True,
         db_use_cached=True,
         lvl_use_cached=True,
-        hdf5_engine="default",
+        memory_mapping_engine="default",
+        load_columns="equilibrium",
         parallel=True,
     ):
         """Fetch the latest databank files from HITRAN or HITEMP with the
@@ -883,9 +932,21 @@ def fetch_databank(
         source: ``'hitran'``, ``'hitemp'``, ``'exomol'``
             [Download database lines from the latest HITRAN (see [HITRAN-2016]_),
             HITEMP (see [HITEMP-2010]_  )] or EXOMOL see [ExoMol-2020]_  ) databases.
-        exomol_database: None
-            if fetching from ''`exomol`'', choose which database to use. Keep
-            ``None`` to use the recommended one. See all available databases
+        database: ``'full'``, ``'range'``, name of an ExoMol database, or ``'default'``
+            if fetching from HITRAN, ``'full'`` download the full database and register
+            it, ``'range'`` download only the lines in the range of the molecule.
+
+            .. note::
+                ``'range'`` will be faster, but will require a new download each time
+                you'll change the range. ``'full'`` is slower and takes more memory, but
+                will be downloaded only once.
+
+            Default is ``'full'``.
+
+            If fetching from HITEMP, only ``'full'`` is available.
+
+            if fetching from ''`exomol`'', use this parameter to choose which database
+            to use. Keep ``'default'`` to use the recommended one. See all available databases
             with :py:func:`radis.io.exomol.get_exomol_database_list`
 
 
@@ -929,12 +990,26 @@ def fetch_databank(
             will be left untouched.
         db_use_cached: bool, or ``'regen'``
             use cached
-        hdf5_engine: ``'pytables'``, ``'vaex'``
+        memory_mapping_engine: ``'pytables'``, ``'vaex'``, ``'feather'``
             which library to use to read HDF5 files (they are incompatible: ``'pytables'`` is
-            row-major while ``'vaex'`` is column-major)
+            row-major while ``'vaex'`` is column-major) or other memory-mapping formats
+            If ``'default'``, use the value from ~/radis.json `["MEMORY_MAPPING_ENGINE"]`
         parallel: bool
             if ``True``, uses joblib.parallel to load database with multiple processes
             (works only for HITEMP files)
+        load_columns: list, ``'all'``, ``'equilibrium'``, ``'noneq'``
+            columns names to load.
+            If ``'equilibrium'``, only load the columns required for equilibrium
+            calculations. If ``'noneq'``, also load the columns required for
+            non-LTE calculations. See :data:`~radis.lbl.loader.drop_all_but_these`.
+            If ``'all'``, load everything. Note that for performances, it is
+            better to load only certain columsn rather than loading them all
+            and dropping them with ``drop_columns``.
+            Default ``'equilibrium'``.
+
+            .. warning::
+                if using ``'equilibrium'``, not all parameters will be available
+                for a Spectrum :py:func:`~radis.spectrum.spectrum.Spectrum.line_survey`.
 
         Notes
         -----
@@ -968,22 +1043,25 @@ def fetch_databank(
             raise NotImplementedError("source: {0}".format(source))
         if source == "hitran":
             dbformat = "hitran"
+            if database == "default":
+                database = "full"
         elif source == "hitemp":
             dbformat = (
                 "hitemp-radisdb"  # downloaded in RADIS local databases ~/.radisdb
             )
+            if database == "default":
+                database = "full"
         elif source == "exomol":
             dbformat = "exomol-radisdb"  # downloaded in RADIS local databases ~/.radisdb  # Note @EP : still WIP.
-        if exomol_database != None:
-            assert source == "exomol"
+
         if [parfuncfmt, source].count("exomol") == 1:
             self.warn(
                 f"Using lines from {source} but partition functions from {parfuncfmt}"
                 + "for consistency we recommend using lines and partition functions from the same database",
                 "AccuracyWarning",
             )
-        if hdf5_engine == "default":
-            hdf5_engine = self.misc.hdf5_engine
+        if memory_mapping_engine == "default":
+            memory_mapping_engine = self.misc.memory_mapping_engine
 
         # Get inputs
         molecule = self.input.molecule
@@ -1013,6 +1091,20 @@ def fetch_databank(
         self.params.db_use_cached = db_use_cached
         self.params.lvl_use_cached = lvl_use_cached
 
+        # %% Which columns to load
+        if load_columns == "equilibrium":
+            columns = list(drop_all_but_these)
+        elif load_columns == "noneq":
+            columns = list(set(drop_all_but_these) | set(required_non_eq))
+        elif load_columns == "all":
+            columns = None  # see fetch_hitemp, fetch_hitran, etc.
+        elif isinstance(load_columns, list):
+            columns = list(set(drop_all_but_these) | set(load_columns))
+        else:
+            raise ValueError(
+                f"Expected a list or 'all' for `load_columns`, got `load_columns={load_columns}"
+            )
+
         # %% Init Line database
         # ---------------------
         self._reset_references()  # bibliographic references
@@ -1021,48 +1113,77 @@ def fetch_databank(
             self.reftracker.add(doi["HITRAN-2016"], "line database")  # [HITRAN-2016]_
             self.reftracker.add(doi["Astroquery"], "data retrieval")  # [Astroquery]_
 
-            if hdf5_engine != "pytables":
-                raise NotImplementedError(f"{hdf5_engine} with ExoMol files")
+            if database == "full":
 
-            # Query one isotope at a time
-            if isotope == "all":
-                raise ValueError(
-                    "Please define isotope explicitely (cannot use 'all' with fetch_databank('hitran'))"
-                )
-            isotope_list = self._get_isotope_list()
+                if isotope == "all":
+                    isotope_list = None
+                else:
+                    isotope_list = ",".join([str(k) for k in self._get_isotope_list()])
 
-            frames = []  # lines for all isotopes
-            for iso in isotope_list:
-                df = fetch_astroquery(
-                    molecule, iso, wavenum_min, wavenum_max, verbose=self.verbose
+                df, local_paths = fetch_hitran(
+                    molecule,
+                    isotope=isotope_list,
+                    load_wavenum_min=wavenum_min,
+                    load_wavenum_max=wavenum_max,
+                    columns=columns,
+                    cache=db_use_cached,
+                    verbose=self.verbose,
+                    return_local_path=True,
+                    engine=memory_mapping_engine,
+                    parallel=parallel,
                 )
-                if len(df) > 0:
-                    frames.append(df)
-                else:
-                    self.warn(
-                        "No line for isotope n°{}".format(iso),
-                        "EmptyDatabaseWarning",
-                        level=2,
+                self.params.dbpath = ",".join(local_paths)
+
+            elif database == "range":
+
+                # Query one isotope at a time
+                if isotope == "all":
+                    raise ValueError(
+                        "Please define isotope explicitely (cannot use 'all' with fetch_databank('hitran'))"
                     )
+                isotope_list = self._get_isotope_list()
 
-            # Merge
-            if frames == []:
-                raise EmptyDatabaseError(
-                    f"{molecule} has no lines on range "
-                    + "{0:.2f}-{1:.2f} cm-1".format(wavenum_min, wavenum_max)
-                )
-            if len(frames) > 1:
-                # Note @dev : may be faster/less memory hungry to keep lines separated for each isotope. TODO : test both versions
-                for df in frames:
-                    assert "iso" in df.columns
-                df = pd.concat(frames, ignore_index=True)  # reindex
+                frames = []  # lines for all isotopes
+                for iso in isotope_list:
+                    df = fetch_astroquery(
+                        molecule, iso, wavenum_min, wavenum_max, verbose=self.verbose
+                    )
+                    if len(df) > 0:
+                        frames.append(df)
+                    else:
+                        self.warn(
+                            "No line for isotope n°{}".format(iso),
+                            "EmptyDatabaseWarning",
+                            level=2,
+                        )
+
+                # Merge
+                if frames == []:
+                    raise EmptyDatabaseError(
+                        f"{molecule} has no lines on range "
+                        + "{0:.2f}-{1:.2f} cm-1".format(wavenum_min, wavenum_max)
+                    )
+                if len(frames) > 1:
+                    # Note @dev : may be faster/less memory hungry to keep lines separated for each isotope. TODO : test both versions
+                    for df in frames:
+                        assert "iso" in df.columns
+                    df = pd.concat(frames, ignore_index=True)  # reindex
+                else:
+                    df = frames[0]
+                self.params.dbpath = "fetched from hitran"
             else:
-                df = frames[0]
-            self.params.dbpath = "fetched from hitran"
+                raise ValueError(
+                    f"Got `database={database}`. When fetching HITRAN, choose `database='full'` to download all database (once for all) or `database='range'` to download only the lines in the current range."
+                )
 
         elif source == "hitemp":
             self.reftracker.add(doi["HITEMP-2010"], "line database")  # [HITEMP-2010]_
 
+            if database != "full":
+                raise ValueError(
+                    f"Got `database={database}`. When fetching HITEMP, only the `database='full'` option is available."
+                )
+
             # Download, setup local databases, and fetch (use existing if possible)
 
             if isotope == "all":
@@ -1075,10 +1196,11 @@ def fetch_databank(
                 isotope=isotope_list,
                 load_wavenum_min=wavenum_min,
                 load_wavenum_max=wavenum_max,
+                columns=columns,
                 cache=db_use_cached,
                 verbose=self.verbose,
                 return_local_path=True,
-                engine=hdf5_engine,
+                engine=memory_mapping_engine,
                 parallel=parallel,
             )
             self.params.dbpath = ",".join(local_paths)
@@ -1092,9 +1214,14 @@ def fetch_databank(
         elif source == "exomol":
             self.reftracker.add(doi["ExoMol-2020"], "line database")  # [ExoMol-2020]
 
+            if database in ["full", "range"]:
+                raise ValueError(
+                    f"Got `database={database}`. When fetching ExoMol, only the `database=` key to retrieve a speciifc database. Use `database='default'` to get the recommended database. See more informatino in radis.io.fetch_exomol()"
+                )
+
             # Download, setup local databases, and fetch (use existing if possible)
-            if hdf5_engine != "pytables":
-                raise NotImplementedError(f"{hdf5_engine} with ExoMol files")
+            if memory_mapping_engine not in ["vaex", "feather"]:
+                raise NotImplementedError(f"{memory_mapping_engine} with ExoMol files")
 
             if isotope == "all":
                 raise ValueError(
@@ -1110,14 +1237,16 @@ def fetch_databank(
             for iso in isotope_list:
                 df, local_path, Z_exomol = fetch_exomol(
                     molecule,
-                    database=exomol_database,
+                    database=database,
                     isotope=iso,
                     load_wavenum_min=wavenum_min,
                     load_wavenum_max=wavenum_max,
+                    columns=columns,
                     cache=db_use_cached,
                     verbose=self.verbose,
                     return_local_path=True,
                     return_partition_function=True,
+                    engine=memory_mapping_engine,
                 )
                 # @dev refactor : have a DatabaseClass from which we load lines and partition functions
                 if len(df) > 0:
@@ -1228,6 +1357,7 @@ def load_databank(
         load_energies=False,
         include_neighbouring_lines=True,
         drop_columns="auto",
+        load_columns="equilibrium",
     ):
         """Loads databank from shortname in the :ref:`Configuration file.
         <label_lbl_config_file>` (`~/radis.json`), or by manually setting all
@@ -1303,6 +1433,20 @@ def load_databank(
             available in :py:meth:`~radis.spectrum.spectrum.Spectrum` method.
             Warning: nonequilibrium calculations are not possible in this mode.
             Default ``'auto'``.
+        load_columns: list, ``'all'``, ``'equilibrium'``, ``'noneq'``
+            columns names to load.
+            If ``'equilibrium'``, only load the columns required for equilibrium
+            calculations. If ``'noneq'``, also load the columns required for
+            non-LTE calculations. See :data:`~radis.lbl.loader.drop_all_but_these`.
+            If ``'all'``, load everything. Note that for performances, it is
+            better to load only certain columsn rather than loading them all
+            and dropping them with ``drop_columns``.
+            Default ``'equilibrium'``.
+
+            .. warning::
+                if using ``'equilibrium'``, not all parameters will be available
+                for a Spectrum :py:func:`~radis.spectrum.spectrum.Spectrum.line_survey`.
+
 
         See Also
         --------
@@ -1335,6 +1479,7 @@ def load_databank(
             db_use_cached,
             lvl_use_cached,
             drop_columns,
+            load_columns,
             load_energies,
             include_neighbouring_lines,
         ) = self._check_database_params(
@@ -1348,8 +1493,9 @@ def load_databank(
             db_use_cached=db_use_cached,
             lvl_use_cached=lvl_use_cached,
             load_energies=load_energies,
-            include_neighbouring_lines=include_neighbouring_lines,
             drop_columns=drop_columns,
+            load_columns=load_columns,
+            include_neighbouring_lines=include_neighbouring_lines,
         )
         # Let's store all params so they can be parsed by "get_conditions()"
         # and saved in output spectra information
@@ -1378,6 +1524,7 @@ def load_databank(
             levelsfmt=levelsfmt,
             db_use_cached=db_use_cached,
             drop_columns=drop_columns,
+            load_columns=load_columns,
             include_neighbouring_lines=include_neighbouring_lines,
         )
         self.misc.total_lines = len(self.df0)  # will be stored in Spectrum metadata
@@ -1415,8 +1562,9 @@ def _check_database_params(
         db_use_cached=None,
         lvl_use_cached=None,
         load_energies=False,
-        include_neighbouring_lines=True,
         drop_columns="auto",
+        load_columns="required",
+        include_neighbouring_lines=True,
     ):
         """Check that database parameters are valid, in particular that paths
         exist. Loads all parameters if a Database from radis.json config file was
@@ -1553,6 +1701,7 @@ def _check_database_params(
             db_use_cached,
             lvl_use_cached,
             drop_columns,
+            load_columns,
             load_energies,
             include_neighbouring_lines,
         )
@@ -1847,7 +1996,8 @@ def _load_databank(
         dbformat,
         levelsfmt,
         db_use_cached,
-        drop_columns="auto",
+        drop_columns,
+        load_columns,
         include_neighbouring_lines=True,
     ) -> pd.DataFrame:
         """Loads all available database files and keep the relevant one.
@@ -1874,7 +2024,20 @@ def _load_databank(
             are dropped, including all information about lines that could be otherwise
             available in :py:meth:`~radis.spectrum.spectrum.Spectrum` method.
             Warning: nonequilibrium calculations are not possible in this mode.
-            Default ``'auto'``.
+        load_columns: list, ``'all'``, ``'equilibrium'``, ``'noneq'``
+            columns names to load.
+            If ``'equilibrium'``, only load the columns required for equilibrium
+            calculations. If ``'noneq'``, also load the columns required for
+            non-LTE calculations. See :data:`~radis.lbl.loader.drop_all_but_these`.
+            If ``'all'``, load everything. Note that for performances, it is
+            better to load only certain columsn rather than loading them all
+            and dropping them with ``drop_columns``.
+            Default ``'equilibrium'``.
+
+            .. warning::
+                if using ``'equilibrium'``, not all parameters will be available
+                for a Spectrum :py:func:`~radis.spectrum.spectrum.Spectrum.line_survey`.
+
 
         Other Parameters
         ----------------
@@ -1907,6 +2070,20 @@ def _load_databank(
                 + drop_auto_columns_for_levelsfmt[levelsfmt]
             )
 
+        # Which columns to load
+        if load_columns == "equilibrium":
+            columns = list(drop_all_but_these)
+        elif load_columns == "noneq":
+            columns = list(set(drop_all_but_these) | set(required_non_eq))
+        elif load_columns == "all":
+            columns = None  # see fetch_hitemp, fetch_hitran, etc.
+        elif isinstance(load_columns, list):
+            columns = list(set(drop_all_but_these) | set(load_columns))
+        else:
+            raise ValueError(
+                f"Expected a list or 'all' for `load_columns`, got `load_columns={load_columns}"
+            )
+
         # subroutine load_and_concat
         # --------------------------------------
         def load_and_concat(files):
@@ -1944,11 +2121,14 @@ def load_and_concat(files):
                             filename,
                             version="hitemp" if dbformat == "cdsd-hitemp" else "4000",
                             cache=db_use_cached,
+                            # load_columns=columns,  # not possible with "pytables-fixed"
                             verbose=verbose,
                             drop_non_numeric=True,
                             load_wavenum_min=wavenum_min,
                             load_wavenum_max=wavenum_max,
+                            engine="pytables",
                         )
+                        # TODO: implement load_columns
                     elif dbformat in ["hitran", "hitemp"]:
                         if dbformat == "hitran":
                             self.reftracker.add(
@@ -1961,10 +2141,12 @@ def load_and_concat(files):
                         df = hit2df(
                             filename,
                             cache=db_use_cached,
+                            # load_columns=columns,  # not possible with "pytables-fixed"
                             verbose=verbose,
                             drop_non_numeric=True,
                             load_wavenum_min=wavenum_min,
                             load_wavenum_max=wavenum_max,
+                            engine="pytables",
                         )
                     elif dbformat in ["hdf5-radisdb", "hitemp-radisdb"]:
                         if dbformat == "hitemp-radisdb":
@@ -1978,6 +2160,7 @@ def load_and_concat(files):
                             )
                         df = hdf2df(
                             filename,
+                            columns=columns,
                             # cache=db_use_cached,
                             verbose=verbose,
                             # drop_non_numeric=True,
@@ -1986,7 +2169,7 @@ def load_and_concat(files):
                             else None,
                             load_wavenum_min=wavenum_min,
                             load_wavenum_max=wavenum_max,
-                            engine="pytables",
+                            engine="guess",
                         )
                     elif dbformat in ["exomol"]:
                         # self.reftracker.add("10.1016/j.jqsrt.2020.107228", "line database")  # [ExoMol-2020]
diff --git a/radis/levels/partfunc.py b/radis/levels/partfunc.py
index a6875a77e..939a1ec33 100644
--- a/radis/levels/partfunc.py
+++ b/radis/levels/partfunc.py
@@ -54,6 +54,7 @@
 # TODO: store molecule_data.json in the H5 file metadata. If not done already.
 
 
+import os
 import sys
 from os.path import exists
 from warnings import warn
@@ -220,6 +221,15 @@ def __init__(self, electronic_state, mode="full summation", verbose=False):
 
         if not mode in ["full summation", "tabulation"]:
             raise ValueError("Choose mode = one of 'full summation', 'tabulation'")
+
+        # vaex processes are stuck if ran from Spyder. See https://github.com/spyder-ide/spyder/issues/16183
+        if mode == "tabulation" and any("SPYDER" in name for name in os.environ):
+            from radis.misc.log import printwarn
+
+            printwarn(
+                "Spyder IDE detected while using parsum_mode='tabulation', which depends on `vaex`. This is the fastest way to compute partition functions, but Vaex processes may be stuck if ran from Spyder. See https://github.com/radis/radis/issues/338. You may consider using another IDE, or use `parsum_mode='full summation'` for the moment (note: starting another iPython console somehow releases the freeze in Spyder) \n"
+            )
+
         self.mode = mode
         self._tab_at = None  # tabulated function
         self._tab_at_noneq = None  # tabulated function
@@ -301,6 +311,8 @@ def _eq_tabulation_setup(self, N_bins):
         import vaex  # import delayed until now (takes ~2s to import)
 
         df = vaex.from_pandas(self.df)
+        # Vaex processs may get stuck in Spyder IDE. https://github.com/radis/radis/issues/338
+        # Temp fix : >>> df.executor.async_method = "awaitio"     (doesn't always work here)
 
         epsilon = 1e-4  # prevent log(0)
         df["logE"] = np.log(df["E"] + epsilon)  # to bin on a log grid
@@ -495,6 +507,8 @@ def _noneq_tabulation_setup(self, N_bins, vib_distribution, rot_distribution):
         import vaex  # import delayed until now (takes ~2s to import)
 
         df = vaex.from_pandas(self.df)
+        # Vaex processs may get stuck in Spyder IDE. https://github.com/radis/radis/issues/338
+        # Temp fix : >>> df.executor.async_method = "awaitio"     (doesn't always work here)
 
         epsilon = 1e-4  # prevent log(0)
         df["logEvib"] = np.log(df["Evib"] + epsilon)  # to bin on a log grid
@@ -850,6 +864,8 @@ def _noneq_3Tvib_tabulation_setup(self, N_bins, vib_distribution, rot_distributi
         import vaex  # import delayed until now (takes ~2s to import)
 
         df = vaex.from_pandas(self.df)
+        # Vaex processs may get stuck in Spyder IDE. https://github.com/radis/radis/issues/338
+        # Temp fix : >>> df.executor.async_method = "awaitio"     (doesn't always work here)
 
         epsilon = 1e-4  # prevent log(0)
 
diff --git a/radis/misc/config.py b/radis/misc/config.py
index d2f6e0259..0ee71508f 100644
--- a/radis/misc/config.py
+++ b/radis/misc/config.py
@@ -39,11 +39,12 @@
 import configparser
 import json
 import os
-import shutil
 import warnings
 from json import JSONDecodeError
 from os.path import dirname, exists, expanduser, join
 
+import hjson
+
 from radis.misc.basics import compare_dict, compare_lists, stdpath
 from radis.misc.utils import DatabankNotFound, getProjectRoot
 from radis.misc.warning import DatabaseAlreadyExists
@@ -66,7 +67,11 @@ def get_config(configpath=CONFIG_PATH_JSON):
     jsonfile = CONFIG_PATH_DEFAULT
     with open(jsonfile) as f:
         try:
-            config = json.load(f)
+            config = hjson.load(f)
+            # read with Hjson to allow comments in CONFIG_PATH_DEFAULT
+            # we do not allow comments in (user) CONFIG_PATH_JSON
+            # because it is read/written automatically. So far there is no
+            # two-way JSON editor with comments. See https://github.com/radis/radis/issues/328
         except JSONDecodeError as err:
             raise JSONDecodeError(
                 "Error reading '{0}' (line {2} col {3}): \n{1}".format(
@@ -391,14 +396,40 @@ def init_radis_json(config_path_json, config_path_old=CONFIG_PATH_OLD):
         convertRadisToJSON(config_path_json, config_path_old)
     else:
         # Create it from the default config path
-        try:
-            shutil.copy2(CONFIG_PATH_DEFAULT, config_path_json)
-        except Exception as err:
-            print(
-                f"Couldn't create RADIS configuration in {config_path_json} ({str(err)})"
-            )
-        else:
-            print(f"RADIS configuration file created in {config_path_json}")
+
+        with open(CONFIG_PATH_DEFAULT) as f:
+            try:
+                default_config = hjson.load(f)
+                # read with Hjson to allow comments in CONFIG_PATH_DEFAULT
+                # we do not allow comments in (user) CONFIG_PATH_JSON
+                # because it is read/written automatically. So far there is no
+                # two-way JSON editor with comments. See https://github.com/radis/radis/issues/328
+            except JSONDecodeError as err:
+                raise JSONDecodeError(
+                    "Error reading '{0}' (line {2} col {3}): \n{1}".format(
+                        CONFIG_PATH_DEFAULT, err.msg, err.lineno, err.colno
+                    ),
+                    err.doc,
+                    err.pos,
+                ) from err
+
+        # Creating json file
+        config = {
+            "_HELLO!": "This is your RADIS configuration file. Any key here will overwrite the values of the radis/default_radis.json file on your local computer (see online version of the latest release: https://github.com/radis/radis/blob/master/radis/default_radis.json). More information on the online docs : https://radis.readthedocs.io/en/latest/lbl/lbl.html#label-lbl-config-file"
+        }
+        config.update(default_config)
+        with open(config_path_json, "w") as outfile:
+            json.dump(config, outfile, indent=3)
+        # outfile.close()
+
+        # try:
+        #     shutil.copy2(CONFIG_PATH_DEFAULT, config_path_json)
+        # except Exception as err:
+        #     print(
+        #         f"Couldn't create RADIS configuration in {config_path_json} ({str(err)})"
+        #     )
+        # else:
+        #     print(f"RADIS configuration file created in {config_path_json}")
 
     # Check it user file exists
     if not exists(config_path_json):
diff --git a/radis/spectrum/spectrum.py b/radis/spectrum/spectrum.py
index a96b774e3..d183e3d64 100644
--- a/radis/spectrum/spectrum.py
+++ b/radis/spectrum/spectrum.py
@@ -999,8 +999,7 @@ def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
         # Get wavespace (in correct unit, and correct medium)
         if wunit == "default":
             wunit = self.get_waveunit()
-        else:
-            wunit = cast_waveunit(wunit)
+        wunit = cast_waveunit(wunit)
         if wunit == "cm-1":
             w = self.get_wavenumber(copy=copy)
         elif wunit == "nm":
diff --git a/radis/spectrum/utils.py b/radis/spectrum/utils.py
index fe99f4fea..4d1a1851d 100644
--- a/radis/spectrum/utils.py
+++ b/radis/spectrum/utils.py
@@ -125,7 +125,7 @@
 
 
 def cast_waveunit(unit, force_match=True):
-    """Standardize unit formats."""
+    """Standardize unit formats, return either "nm", "nm_vac" or "cm-1"."""
     if unit in WAVELEN_UNITS:
         return "nm"
     if unit in WAVELENVAC_UNITS:
diff --git a/setup.py b/setup.py
index 0d3b460cd..25524ad28 100644
--- a/setup.py
+++ b/setup.py
@@ -245,12 +245,12 @@ def run_setup(with_binary):
             "hitran-api",
             "beautifulsoup4",  # parse ExoMol website
             "lxml",  # parser used for ExoMol website
-            "pyarrow",  # for the feather format (temporarily needed for ExoMol)
             "numpy",
             "scipy>=1.4.0",
             "matplotlib",  # ">=3.4.0" to suppress the Ruler warning, but only available for Python >= 3.7
             "seaborn",  # other matplotlib themes
             "cython",
+            "hjson",
             "pandas>=1.0.5",
             "plotly>=2.5.1",
             "progressbar2",  # used in vaex
@@ -259,7 +259,6 @@ def run_setup(with_binary):
             "astropy",  # Unit aware calculations
             "publib>=0.3.2",  # Plotting styles for Matplotlib
             "plotly>=2.5.1",  # for line survey HTML output
-            "termcolor",  # terminal colors
             "configparser",
             "astroquery>=0.3.9",  # to fetch HITRAN databases
             "json-tricks>=3.15.0",  # to deal with non jsonable formats
@@ -269,7 +268,7 @@ def run_setup(with_binary):
             "numba",  # just-in-time compiler
             "psutil",  # for getting user RAM
             "tuna",  # to generate visual/interactive performance profiles
-            "vaex",  # HDF5,
+            "vaex>=4.4.0",  # load HDF5 files  (version for custom HDF5 groups)
             "h5py",  # HDF5
             "habanero",  # CrossRef API to retrieve data from doi
         ],


Test patch:
diff --git a/radis/test/io/test_hitemp.py b/radis/test/io/test_hitemp.py
index 7350a6f61..678114927 100644
--- a/radis/test/io/test_hitemp.py
+++ b/radis/test/io/test_hitemp.py
@@ -98,15 +98,18 @@ def test_fetch_hitemp_partial_download_CO2(verbose=True, *args, **kwargs):
     # ... This is done at the HITEMPDatabaseManager level
     # ... unitest for this part:
     ldb = HITEMPDatabaseManager(
-        name="HITEMP-CO2", molecule="CO2", local_databases="~/.radisdb/hitemp/"
+        name="HITEMP-CO2",
+        molecule="CO2",
+        engine="default",
+        local_databases="~/.radisdb/hitemp/",
     )
-    local_files, _ = ldb.get_filenames(engine="auto")
+    local_files, _ = ldb.get_filenames()
     relevant_file, file_wmin, file_wmax = keep_only_relevant(
         local_files, wavenum_min=wmin, wavenum_max=wmax
     )
 
     assert len(relevant_file) == 1
-    assert basename(relevant_file[0]) == "CO2-02_02500-03000_HITEMP2010.h5"
+    assert basename(relevant_file[0]).startswith("CO2-02_02500-03000_HITEMP2010.")
     assert file_wmin == 2500
     assert file_wmax == 3000
 
@@ -123,7 +126,7 @@ def test_fetch_hitemp_partial_download_CO2(verbose=True, *args, **kwargs):
     assert "HITEMP-CO2" in getDatabankList()
 
     assert len(local_files) == 1
-    assert basename(local_files[0]) == "CO2-02_02500-03000_HITEMP2010.h5"
+    assert basename(local_files[0]).startswith("CO2-02_02500-03000_HITEMP2010.")
 
 
 @pytest.mark.needs_connection
@@ -161,7 +164,7 @@ def test_fetch_hitemp_all_molecules(molecule, verbose=True, *args, **kwargs):
         columns=["int", "wav"],
         verbose=verbose,
         return_local_path=True,
-        engine="pytables",
+        engine="default",
     )
 
     assert f"HITEMP-{molecule}" in getDatabankList()
@@ -170,6 +173,7 @@ def test_fetch_hitemp_all_molecules(molecule, verbose=True, *args, **kwargs):
         name=f"HITEMP-{molecule}",
         molecule=molecule,
         verbose=verbose,
+        engine="default",
         local_databases="~/.radisdb/hitemp/",
     )
     url, Nlines, _, _ = ldb.fetch_url_Nlines_wmin_wmax()
@@ -183,25 +187,39 @@ def test_partial_loading(*args, **kwargs):
 
     Also check 'vaex' engine and ``radis.config["AUTO_UPDATE_DATABASE"]``"""
 
-    wmin, wmax = 2500, 4500
+    from radis.test.utils import getTestFile
 
-    # First ensures that the database is wider than this (else test is irrelevant) :
-    df = fetch_hitemp("OH", engine="pytables")
-    assert df.wav.min() < wmin
-    assert df.wav.max() > wmax
+    df = fetch_hitemp(
+        "OH",
+        local_databases=getTestFile("."),
+        databank_name="HITEMP-OH-TEST-PARTIAL-LOADING",
+        engine="pytables",
+    )
 
     # Now fetch with partial loading
+    wmin, wmax = 2500, 4500
+    # ... First ensures that the database is wider than this (else test is irrelevant) :
+    assert df.wav.min() < wmin
+    assert df.wav.max() > wmax
+    # ... load :
     df = fetch_hitemp(
-        "OH", load_wavenum_min=wmin, load_wavenum_max=wmax, engine="pytables"
+        "OH",
+        local_databases=getTestFile("."),
+        databank_name="HITEMP-OH-TEST-PARTIAL-LOADING",
+        load_wavenum_min=wmin,
+        load_wavenum_max=wmax,
+        engine="pytables",
     )
     assert df.wav.min() >= wmin
     assert df.wav.max() <= wmax
 
     # Test with isotope:
-    wmin2 = 0
+    wmin2 = 1
     wmax2 = 300
     df = fetch_hitemp(
         "OH",
+        local_databases=getTestFile("."),
+        databank_name="HITEMP-OH-TEST-PARTIAL-LOADING",
         load_wavenum_min=wmin2,
         load_wavenum_max=wmax2,
         isotope="2",
@@ -210,6 +228,8 @@ def test_partial_loading(*args, **kwargs):
     assert df.iso.unique() == 2
     df = fetch_hitemp(
         "OH",
+        local_databases=getTestFile("."),
+        databank_name="HITEMP-OH-TEST-PARTIAL-LOADING",
         load_wavenum_min=wmin2,
         load_wavenum_max=wmax2,
         isotope="1,2",
@@ -223,21 +243,31 @@ def test_partial_loading(*args, **kwargs):
     old_config = radis.config["AUTO_UPDATE_DATABASE"]
     try:
         radis.config["AUTO_UPDATE_DATABASE"] = True
-        df = fetch_hitemp("OH", engine="vaex")
-    finally:
-        radis.config["AUTO_UPDATE_DATABASE"] = old_config
-    assert df.wav.min() < wmin
-    assert df.wav.max() > wmax
-
-    # multiple isotope selection not implemetned with vaex
-    with pytest.raises(NotImplementedError):
-        fetch_hitemp(
+        df = fetch_hitemp(
             "OH",
-            load_wavenum_min=wmin,
-            load_wavenum_max=wmax,
-            isotope="1,2,3",
+            local_databases=getTestFile("."),
+            databank_name="HITEMP-OH-TEST-PARTIAL-LOADING",
             engine="vaex",
         )
+    finally:
+        radis.config["AUTO_UPDATE_DATABASE"] = old_config
+    # assert that database is wider than future selection
+    assert df.wav.min() < wmin2
+    assert df.wav.max() > wmax2
+
+    # now test wrange & multiple isotope selection with vaex
+    df = fetch_hitemp(
+        "OH",
+        local_databases=getTestFile("."),
+        databank_name="HITEMP-OH-TEST-PARTIAL-LOADING",
+        load_wavenum_min=wmin2,
+        load_wavenum_max=wmax2,
+        isotope="2,3",
+        engine="vaex",
+    )
+    assert df.wav.min() >= wmin2
+    assert df.wav.max() <= wmax2
+    assert set(df.iso.unique()) == {2, 3}
 
 
 @pytest.mark.needs_connection
diff --git a/radis/test/io/test_query.py b/radis/test/io/test_query.py
index 6c359f931..ed5a775e5 100644
--- a/radis/test/io/test_query.py
+++ b/radis/test/io/test_query.py
@@ -99,11 +99,52 @@ def test_fetch_astroquery_cache(verbose=True, *args, **kwargs):
     return
 
 
+# ignored by pytest with argument -m "not needs_connection"
+@pytest.mark.needs_connection
+def test_fetch_hitran(*args, **kwargs):
+
+    from radis.io.hitran import fetch_hitran
+
+    df = fetch_hitran("CO")
+
+    assert len(df) == 5381
+    assert df.wav.min() == 3.40191
+    assert df.wav.max() == 14477.377142
+
+
+@pytest.mark.needs_connection
+def test_calc_hitran_spectrum(verbose=True, plot=False, *args, **kwargs):
+    """
+    Test full range & partial loading of CO spectrum
+    """
+
+    from radis import test_spectrum
+
+    # The main thing here is to make sure the two function calls work :
+    s = test_spectrum(databank=("hitran", "full"), name="full range", verbose=verbose)
+    s2 = test_spectrum(
+        databank=("hitran", "range"), name="partial range", verbose=verbose
+    )
+
+    if plot:
+        from radis import plot_diff
+
+        plot_diff(s, s2)
+
+    # assert s == s2
+    # Note: tiny differences (0.56%) probably due to how data is stored on disk
+    assert s.compare_with(s2, spectra_only="abscoeff", method="ratio", rtol=0.006)
+
+    return
+
+
 def _run_testcases(verbose=True, *args, **kwargs):
 
     test_fetch_astroquery(verbose=verbose, *args, **kwargs)
     test_fetch_astroquery_empty(verbose=verbose, *args, **kwargs)
     test_fetch_astroquery_cache(verbose=verbose, *args, **kwargs)
+    test_fetch_hitran(*args, **kwargs)
+    test_calc_hitran_spectrum(*args, **kwargs)
 
     return True
 
diff --git a/radis/test/lbl/test_base.py b/radis/test/lbl/test_base.py
index a4a5a9850..d4403e34e 100644
--- a/radis/test/lbl/test_base.py
+++ b/radis/test/lbl/test_base.py
@@ -241,7 +241,7 @@ def test_populations_CO2_hamiltonian(
         verbose=verbose,
     )
     sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
-    sf.load_databank("HITEMP-CO2-HAMIL-TEST", load_energies=True)
+    sf.load_databank("HITEMP-CO2-HAMIL-TEST", load_energies=True, load_columns="noneq")
 
     # First run a calculation at equilibrium
     s = sf.eq_spectrum(300)
@@ -338,7 +338,7 @@ def test_optically_thick_limit_1iso(verbose=True, plot=True, *args, **kwargs):
             verbose=False,
         )
         sf.warnings["NegativeEnergiesWarning"] = "ignore"
-        sf.load_databank("HITEMP-CO2-TEST")
+        sf.load_databank("HITEMP-CO2-TEST", load_columns="noneq")
         pb = ProgressBar(3, active=verbose)
         s_eq = sf.eq_spectrum(Tgas=Tgas, mole_fraction=1, name="Equilibrium")
         pb.update(1)
@@ -450,7 +450,7 @@ def test_optically_thick_limit_2iso(verbose=True, plot=True, *args, **kwargs):
             pressure=P,
             verbose=False,
         )
-        sf.load_databank("HITEMP-CO2-TEST")
+        sf.load_databank("HITEMP-CO2-TEST", load_columns="noneq")
         sf.warnings["NegativeEnergiesWarning"] = "ignore"
         #        sf.warnings['MissingSelfBroadeningWarning'] = 'ignore'
         pb = ProgressBar(3, active=verbose)
@@ -507,6 +507,8 @@ def test_optically_thick_limit_2iso(verbose=True, plot=True, *args, **kwargs):
 
 def test_get_waverange(*args, **kwargs):
 
+    assert get_waverange(1 * u.um, 10 * u.um, medium="vacuum") == (1000, 10000)
+
     # 'wunit' is none
     # ...'wmin/wmax' is none
     # ...... wavenumber is passed > should return the same wavenumbers
diff --git a/radis/test/lbl/test_broadening.py b/radis/test/lbl/test_broadening.py
index a316beb1c..4aaed7b1c 100644
--- a/radis/test/lbl/test_broadening.py
+++ b/radis/test/lbl/test_broadening.py
@@ -91,7 +91,7 @@ def test_broadening_vs_hapi(rtol=1e-2, verbose=True, plot=False, *args, **kwargs
             "HighTemperatureWarning": "ignore",
             "GaussianBroadeningWarning": "ignore",
         },
-    )  # 0.2)
+    )
     sf.load_databank(
         path=join(hapi_data_path, "CO.data"), format="hitran", parfuncfmt="hapi"
     )
@@ -318,6 +318,7 @@ def test_broadening_DLM(verbose=True, plot=False, *args, **kwargs):
         pressure=p,
         truncation=truncation,
         neighbour_lines=neighbour_lines,
+        optimization=None,
         isotope="1",
         verbose=False,
         warnings={
@@ -326,16 +327,16 @@ def test_broadening_DLM(verbose=True, plot=False, *args, **kwargs):
             "HighTemperatureWarning": "ignore",
             "GaussianBroadeningWarning": "ignore",
         },
-    )  # 0.2)
+    )
     sf.load_databank("HITRAN-CO-TEST")
 
     # Reference: calculate without DLM
-    assert sf.misc["chunksize"] is None
+    assert sf.params["optimization"] is None
     s_ref = sf.eq_spectrum(Tgas=T)
     s_ref.name = "Reference ({0:.2f}s)".format(s_ref.conditions["calculation_time"])
 
     # DLM:
-    sf.misc["chunksize"] = "DLM"
+    sf.params["optimization"] = "simple"
     sf.params.broadening_method = "convolve"
     s_dlm = sf.eq_spectrum(Tgas=T)
     s_dlm.name = "DLM ({0:.2f}s)".format(s_dlm.conditions["calculation_time"])
@@ -398,14 +399,14 @@ def test_broadening_DLM_FT(verbose=True, plot=False, *args, **kwargs):
         neighbour_lines=neighbour_lines,
         isotope="1",
         verbose=verbose,
-        chunksize="DLM",
+        optimization="simple",
         warnings={
             "MissingSelfBroadeningWarning": "ignore",
             "NegativeEnergiesWarning": "ignore",
             "HighTemperatureWarning": "ignore",
             "GaussianBroadeningWarning": "ignore",
         },
-    )  # 0.2)
+    )
     sf.load_databank("HITRAN-CO-TEST")
 
     # DLM, real space
@@ -478,11 +479,11 @@ def test_broadening_DLM_noneq(verbose=True, plot=False, *args, **kwargs):
             "HighTemperatureWarning": "ignore",
             "GaussianBroadeningWarning": "ignore",
         },
-    )  # 0.2)
-    sf.load_databank("HITRAN-CO2-TEST")
+    )
+    sf.load_databank("HITRAN-CO2-TEST", load_columns="noneq")
 
     # DLM:
-    sf.misc["chunksize"] = "DLM"
+    sf.params["optimization"] = "simple"
     s_dlm_eq = sf.eq_spectrum(Tgas=3000)
     s_dlm_eq.name = "DLM eq ({0:.2f}s)".format(s_dlm_eq.conditions["calculation_time"])
 
@@ -830,9 +831,8 @@ def test_noneq_continuum(plot=False, verbose=2, warnings=True, *args, **kwargs):
         }
     )
     sf.fetch_databank(
-        "hitran"
+        "hitran", load_columns="noneq"
     )  # uses HITRAN: not really valid at this temperature, but runs on all machines without install
-    #        sf.load_databank('HITEMP-CO2-DUNHAM')       # to take a real advantage of abscoeff continuum, should calculate with HITEMP
     sf._export_continuum = True  # activate it
 
     # Calculate one without pseudo-continuum
diff --git a/radis/test/lbl/test_calc.py b/radis/test/lbl/test_calc.py
index 67745e9b7..f4d053846 100644
--- a/radis/test/lbl/test_calc.py
+++ b/radis/test/lbl/test_calc.py
@@ -173,6 +173,10 @@ def test_calc_spectrum(verbose=True, plot=True, warnings=True, *args, **kwargs):
     #         0.11158052, 0.04583138, 0.00232573, 0.00136583])
     # Update on 05/11/21 (Radis 0.10.4) with use of tabulated instead of ab-initio
     # partition functions for Qref  (0.02% constant difference)
+    #
+    # Update on 14/11/21: after switching from HITRAN fetched by astroquery (partial range)
+    # to HITRAN fetched by HAPI (full database) > atol=1e-6 fails, because of
+    # numerical errors in how the two databases are stored. Switched to rtol=1e-3
     w_ref = np.array(
         [
             4197.60321744,
@@ -228,7 +232,7 @@ def test_calc_spectrum(verbose=True, plot=True, warnings=True, *args, **kwargs):
     # [71:-71] because of the shift introduced in 0.9.30 where convolved
     # arrays are not cropped; but np.nan values are added
     assert np.allclose(w[71:-71][::100], w_ref, atol=1e-6)
-    assert np.allclose(I[71:-71][::100], I_ref, atol=1e-6)
+    assert np.allclose(I[71:-71][::100], I_ref, rtol=1e-3)
 
     return True
 
@@ -347,6 +351,10 @@ def test_calc_spectrum_overpopulations(
     #        0.3293075 , 0.13525292, 0.00645308, 0.00366828])
     # Update on 05/11/21 (Radis 0.10.4) with use of tabulated instead of ab-initio
     # partition functions for Qref  (0.02% constant difference)
+    #
+    # Update on 14/11/21: after switching from HITRAN fetched by astroquery (partial range)
+    # to HITRAN fetched by HAPI (full database) > atol=1e-6 fails, because of
+    # numerical errors in how the two databases are stored. Switched to rtol=1e-3
 
     I_ref = np.array(
         [
@@ -378,7 +386,7 @@ def test_calc_spectrum_overpopulations(
 
     # [71:-71] because of the shift introduced in 0.9.30 where convolved
     # arrays are not cropped; but np.nan values are added
-    assert np.allclose(I[71:-71][::100], I_ref, atol=1e-6)
+    assert np.allclose(I[71:-71][::100], I_ref, rtol=1e-3)
 
     if verbose:
         printm("Test overpopulations: OK")
@@ -514,7 +522,7 @@ def test_all_calc_methods_CO2pcN(
     del database_kwargs["info"]
     database_kwargs["levelsfmt"] = "cdsd-pcN"
     # ... load the new database
-    sf.load_databank(**database_kwargs, load_energies=True)
+    sf.load_databank(**database_kwargs, load_energies=True, load_columns="noneq")
 
     # Now, define Evib:
     Q_calc = sf.parsum_calc["CO2"][1]["X"]
@@ -593,7 +601,7 @@ def test_eq_vs_noneq_isotope(verbose=True, plot=False, warnings=True, *args, **k
     sf.warnings["NegativeEnergiesWarning"] = "ignore"
     sf.warnings["HighTemperatureWarning"] = "ignore"
     sf.fetch_databank(
-        "hitran"
+        "hitran", load_columns="noneq"
     )  # uses HITRAN: not really valid at this temperature, but runs on all machines without install
     s_nq = sf.non_eq_spectrum(Tvib=Tgas, Trot=Tgas, name="Non-eq")
     s_eq = sf.eq_spectrum(Tgas=Tgas, name="Eq")
diff --git a/radis/test/lbl/test_loader.py b/radis/test/lbl/test_loader.py
index 925d0601d..8c21b19c9 100644
--- a/radis/test/lbl/test_loader.py
+++ b/radis/test/lbl/test_loader.py
@@ -57,7 +57,8 @@ def test_retrieve_from_database(
         )
         sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
         sf.init_databank(
-            "HITEMP-CO2-TEST"
+            "HITEMP-CO2-TEST",
+            load_columns="noneq",
         )  # unlike load_databank, will automatically be built when needed
         db = sf.init_database(temp_database_name, autoretrieve=True)
 
@@ -180,7 +181,7 @@ def test_custom_abundance(verbose=True, plot=False, *args, **kwargs):
         verbose=0,
     )
     sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
-    sf.load_databank("HITEMP-CO2-TEST")
+    sf.load_databank("HITEMP-CO2-TEST", load_columns="all")
 
     #%% Compute a spectrum ; set the terrestrial abundance to False to
     # force computation of linestrength from Einstein coefficients instead