Repository: radis/radis (language: python)

Issue:
[Reduced Memory Use] calculating spectrum using vaex
This  PR addresses memory use for computing Spectrum, existing Pandas implementation though efficient for small databases, doesn't perform well on larger line database.
To solve this issue, vaex is used which performs exceptionally well for memory uses.
The time efficiency of vaex is also better than Pandas.

![image](https://github.com/radis/radis/assets/78527441/9e04fcbf-2ce2-474d-9289-a6a808453543)



Fixes #118 

Reduce Memory usage
Radis is quite memory hungry but there is room for improvement. 

## Still using Pandas :

Some columns do not need to be loaded: 

- [ ] Linestrengths "int" can be recomputed from Einstein coefficients at runtime, or the other way around. 
- [x] molecule `id` is not needed as we anyway compute only one molecule at a time. Same goes for isotope if only one isotope ? 
- [ ] A good starting point may be to list all the columns present in sf.df0 and sf.df1 at the end of a SpectrumFactory calculation, and remove all that are not needed. 

Also : 

- [x] do not fill a full column "Qgas" with partition functions  (an internal dictionary could be more than enough). Same for all molecular parameters (mass, etc. )  #287 

- [ ] Some column formats can be reduced (float64 > float32, int64 > uint8) with very limited impact, or formats such as "categorical" can be used. 

- [ ] the `.df0` "buffer database" can be fully removed and replaced by a .h5 file : #176 
> a very limited H5 file for each calculation could totally replace the current mechanism where the SpectrumFactory.df0 Pandas dataframe containing all lines loaded, is copied at calculation time to SpectrumFactory.df1 (where populations will be computed based on temperature, etc.). Using a .h5 file may induce a (limited ?) performance overhead, but we could also divide the memory usage by up to 2


---

## Moving away from Pandas 

Some other libraries may be more adapted to 'Big Data' (which I define here as : doesn't fit in RAM) 

- [vaex](https://towardsdatascience.com/how-to-analyse-100s-of-gbs-of-data-on-your-laptop-with-python-f83363dda94)    
- Dask

> Dask is more focused on scaling the code to compute clusters, while Vaex makes it easier to work with large datasets on a single machine. [from here](https://towardsdatascience.com/are-you-still-using-pandas-to-process-big-data-in-2021-850ab26ad919)

Let's try them ? 

---

This is discussed on [Slack](https://radis.github.io/slack-invite/) in channel [#memory_usage](https://radis-radiation.slack.com/archives/C0193Q836N4)

Anyone wanting to work on that topic is very welcome !!


---

Check memory consumption in [radis-benchmark](https://radis.github.io/radis-benchmark/)

Crash test / suggested test case : 

- read all HITEMP CH4 database in memory : 
```
31 880 412 lines,  0-13400cm-1 , 5 GB .par file,  1.9 GB h5-compressed file.
```
> 27/01/21 : can't get it to work on a 16 GB RAM laptop. 

- [🥇 Challenge : compute H2O on Radis-Lab](https://github.com/radis/radis-lab/issues/24)

--- 

🔗 _Dedicated 💬 Slack channel : https://radis-radiation.slack.com/archives/C0193Q836N4_

🔗 [_GSOC 2021 Project - Reduce Memory Usage_](https://openastronomy.org/gsoc/gsoc2021/#/projects?project=memory_usage)


Discussion / hints:
## [Codecov](https://app.codecov.io/gh/radis/radis/pull/576?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) Report
> Merging [#576](https://app.codecov.io/gh/radis/radis/pull/576?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) (a41d987) into [develop](https://app.codecov.io/gh/radis/radis/commit/0b87b6f9e0541fce9f6c9e8b7be30ffefbdd3c30?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) (0b87b6f) will **increase** coverage by `0.40%`.
> The diff coverage is `73.78%`.

:exclamation: Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please [install the Github App Integration](https://github.com/apps/codecov) for your organization. [Read more](https://about.codecov.io/blog/codecov-is-updating-its-github-integration/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis).

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop     #576      +/-   ##
===========================================
+ Coverage    73.90%   74.31%   +0.40%     
===========================================
  Files          144      144              
  Lines        19643    20009     +366     
===========================================
+ Hits         14517    14869     +352     
- Misses        5126     5140      +14     
```

</details>
Hello, 
Here is first a review of the results, I didn't check the inside of your work and the coding quality yet.

**1/** you should try to sort the following performance decrease for small spectra. On my machine, the CO spectrum simulated from 2000 to 2010 cm-1 takes 0.24 s with `vaex` and 0.04 s with `pandas`
- [ ] Improve time performance for small spectra

```
from radis import calc_spectrum
s = calc_spectrum(2000, 2010,         # cm-1
                  molecule='CO',
                  isotope='1,2,3',
                  pressure=1.01325,   # bar
                  Tgas=1000,
                  mole_fraction=0.1,
                  databank='hitran',  # or 'hitemp'
                  diluent = "air",     
                  verbose = 1,
                  engine = "vaex"
                  )
s.plot('absorbance')
```
with `pandas`:
```python
----------------------------------------
0.04s - Spectrum calculated
```
with `vaex`:
```python
----------------------------------------
0.24s - Spectrum calculated
```

**2/** When I tried to find the origin of the time difference, I ran the `calc_spectrum` function with `verbose = 3` and I got a bug. Could you please implement more tests of this new implementation and sort the bug?

- [ ] Add more tests
- [x] Debug `verbose = 3` in `calc_spectrum`
 


```python
from radis import calc_spectrum
s = calc_spectrum(2000, 2010,         # cm-1
                  molecule='CO',
                  isotope='1,2,3',
                  pressure=1.01325,   # bar
                  Tgas=1000,
                  mole_fraction=0.1,
                  databank='hitran',  # or 'hitemp'
                  diluent = "air",     
                  verbose = 3,
                  engine = "vaex"
                  )
s.plot('absorbance')
```
```python
 File "c:\users\nicolas minesi\python\radis\api\untitled0.py", line 9, in <module>
    s = calc_spectrum(2000, 2010,         # cm-1

  File "c:\users\nicolas minesi\python\radis\lbl\calc.py", line 537, in calc_spectrum
    generated_spectrum = _calc_spectrum_one_molecule(

  File "c:\users\nicolas minesi\python\radis\lbl\calc.py", line 902, in _calc_spectrum_one_molecule
    s = sf.eq_spectrum(

  File "c:\users\nicolas minesi\python\radis\lbl\factory.py", line 802, in eq_spectrum
    self._cutoff_linestrength()

  File "c:\users\nicolas minesi\python\radis\lbl\base.py", line 3265, in _cutoff_linestrength
    Nlines_cutoff / len(df.S) * 100, cutoff

TypeError: object of type 'Expression' has no len()
```
For smaller data frames, I tested some simple operations like sorting and finding the minimum on a simple data frame.
Link to Kaggle notebook https://www.kaggle.com/code/someshverma/notebook2e47fa84fe
And, it clearly shows that in general for smaller data frames performance of pandas is better.

I also tested this on the local system and similar results are obtained. So, It appears that for smaller data frame performance of Pandas is better than Vaex . Though the performance of vaex improves with larger data frames.

I will try if I can improve Vaex implementation for smaller spectra.
I already have some test cases , i will commit those and add more test cases.

I have resolved issue related to verbose in my latest PR
I pushed a minor commit to show you that some improvements can still be made (in this case, how the sum is computed). I know it's a tedious task, but it is surely valuable for RADIS (but also for you to learn). So please try to improve each line of code taking some time in the computation. On my computer, the commit I pushed the time taken to calculate 25 spectra decreased from `6.0 s `to `5.4 s`. I strongly recommend you to use `verbose=3` to track where the time is taken.

My code to see the impact of a change: 
```python
from radis import calc_spectrum
import time

t0 = time.time()
for i in range(25):
    s = calc_spectrum(2000, 2010,         # cm-1
                      molecule='CO',
                      isotope='1,2,3',
                      pressure=1.01325,   # bar
                      Tgas=1000,
                      mole_fraction=0.1,
                      databank='hitemp',  # or 'hitemp'
                      diluent = "air",     
                      verbose = 0,
                      engine = "vaex"
                      )
t1 = time.time()
print(t1 -t0)
```
Thanks for improving the sum calculations, I am working on reducing the computation time.
Last commit reduces calculation time for the code below by 10 seconds

```
from radis import calc_spectrum

import time
t0=time.time()

s = calc_spectrum(1500, 2500,         # cm-1
                      molecule='H2O',
                      isotope='1,2,3',
                      pressure=1.01325,   # bar
                      Tgas=700,           # K
                      mole_fraction=0.1,
                      wstep='auto',
                      path_length=1,      # cm
                      databank='hitemp',  # or 'hitemp', 'geisa', 'exomol'
                      engine='vaex',
                      )

s.apply_slit(0.5, 'nm')       # simulate an experimental slit

t1=time.time()

print('Time taken : '+str(t1 - t0))

```


And reduces 0.5 seconds  calculation time for the code 

```
from radis import calc_spectrum
import time
t0=time.time()
s = calc_spectrum(2000, 2010,         # cm-1
                  molecule='CO',
                  isotope='1,2,3',
                  pressure=1.01325,   # bar
                  Tgas=1000,
                  mole_fraction=0.1,
                  databank='hitran',  # or 'hitemp'
                  diluent = "air",
                  verbose = 3,
                  engine = "vaex"
                  )
t1=time.time()


print('Time taken : '+str(t1-t0))
```

I am working to further improve the results.
Very interesting last commit, we were actually talking about something similar at lunch with @erwanp.

@1someshverma, so you define a threshold when the database is either converted for the *calculations* into pandas. This threshold is 2e5 entries in the database. Please explain how you picked this threshold. A graph might be sufficient. If you used a code to determine this threshold, could you please share it?
Concerning your test cases: they are not passing on my machine!
If I run the test `test_vaex_and_pandas_dataframe_load_databank()`, then the test fails at `assert np.all(df1[column].to_numpy() == df2[column])` because the column `jl` is empty in the vaex format. Could you please check this?
I have tested this on my laptop and it passed the test case , i will test it again and sort this issue.
@minouHub After some of the improvement in the time performance for Vaex , now it perform better than better pandas after 0.2 M lines . I will share the graph in support of this soon. 
@minouHub  `test_vaex_and_pandas_dataframe_load_databank()`  is passing on my laptop and also on Travis. Can you share Python version and vaex version on your system so that I can check if it is related to different versions .
These are the updated time graph, and it shows that the total time taken for vaex starts to perform better than pandas after 9e5 (maybe I made some mistake in computing the time earlier). 
I will update the threshold with 9e5 .

As you can see compared to the previous time graph, the time efficiency is improved for Vaex both for calculations and overall.
![image](https://github.com/radis/radis/assets/78527441/548123e2-d23f-4e47-a4cb-b302927020dd)
![image](https://github.com/radis/radis/assets/78527441/22a2290d-a4f1-421f-a989-5fdbe9715817)

> @minouHub `test_vaex_and_pandas_dataframe_load_databank()` is passing on my laptop and also on Travis. Can you share Python version and vaex version on your system so that I can check if it is related to different versions .

You must be right, if Travis is passing I should be passing as well. I'll try to sort this on my side. 

> These are the updated time graph, and it shows that the total time taken for vaex starts to perform better than pandas after 9e5 (maybe I made some mistake in computing the time earlier).

Could you please the code to generate these figure so I can check if the same order of magnitude is obtained on other computer?
This is the code that I used for obtaining the time taken for calculating the spectrum 

```
from radis import calc_spectrum

import time
t0=time.time()

s = calc_spectrum(2450, 2500,         # cm-1
                      molecule='H2O',
                      isotope='1,2,3',
                      pressure=1.01325,   # bar
                      Tgas=700,           # K
                      mole_fraction=0.1,
                      wstep='auto',
                      path_length=1,      # cm
                      databank='hitemp',  # or 'hitemp', 'geisa', 'exomol'
                      engine='vaex',
                      )

s.apply_slit(0.5, 'nm')       # simulate an experimental slit

t1=time.time()

print('Time taken : '+str(t1 - t0))

```

I  started with wmin = 2450 and then decreased it in steps to 10 then increased step length.

I think it may be slightly machine dependent but the difference should be quite small.
I did some calculations to get insights on how to improve database loading & RAM use : see [radis-benchmark/manual_benchmarks/test_database_loading.ipynb](https://github.com/radis/radis-benchmark/blob/master/manual_benchmarks/test_database_loading.ipynb)
Featured in GSOC 2021 project [Reduce Memory Usage 🚀 ](https://openastronomy.org/gsoc/gsoc2021/#/projects?project=memory_usage)
>  Linestrengths "int" can be recomputed from Einstein coefficients at runtime, or the other way around.

Are there other parameters as well like these that are related to each other or it just these two?

> Some column formats can be reduced (float64 > float32, int64 > uint8) with very limited impact, or formats such as "categorical" can be used.

Adjusting dtypes and precision part of https://github.com/radis/radis-benchmark/blob/master/manual_benchmarks/test_database_loading.ipynb concluded that "adjusting precision is not considered at this point". Wanted to confirm this is indeed the case. 

> do not fill a full column "Qgas" with partition functions (an internal dictionary could be more than enough). Same for all molecular parameters (mass, etc. )

I am sorry. I couldn't get this point. Could you please elaborate?  
Hello @gagan-aryan. Many questions here: 

1. I don't think of anything else right now. List of parameters can be found in the [Hitran parser](https://github.com/radis/radis/blob/develop/radis/io/hitran.py#L64), they are all independants (except from A & int). 

2. Yes, adjusting precision did not help !

3. I realized this requires extra documentation and knowledge of the equations behind the code. I just pushed a new PR #225 to show the equations populating the levels in the docs. This will help : 
 
- in [calc_populations_eq()](https://radis.readthedocs.io/en/latest/source/radis.lbl.base.html#radis.lbl.base.BaseFactory.calc_populations_eq), "Qgas", called the partition function, is used to compute the population of the emitting/absorbing level of a line. However, although `gu` or `Eu` are unique for each line, Qgas it is the same for all levels of a same molecule & isotope. 
- In initial RADIS version, because I mostly cared for CPU performance and not Memory usage, I would add an extra column so that each line has a corresponding Qgas. Then I would just have to multiply/divide the columns. Fast, but these potentially millions of lines contain the same value for Qgas, or maybe 2 or 3 if we have different isotopes. 
- In the code Qgas is added here: https://github.com/radis/radis/blob/develop/radis/lbl/base.py#L2135 , and later used in https://github.com/radis/radis/blob/develop/radis/lbl/base.py#L2166
- Instead of doing this, we could simply read the content of a dictionary. 
New benchmark test on extra columns to remove + use of vaex
(see previous benchmark for CO : https://github.com/radis/radis/pull/306#issuecomment-896426701 ) 

Test code:
```python
import radis
from radis import SpectrumFactory
import radis
radis.config['AUTO_UPDATE_DATABASE']= True

sf = SpectrumFactory(1900, 2300,         # cm-1
                  molecule='H2O',
                  isotope='1',
                  pressure=1.01325,   # bar
                  cutoff=1e-25,
                  verbose=3,)
sf.fetch_databank('hitemp', hdf5_engine="vaex")  # or "pytables"              
sf.eq_spectrum(300,)
```

**Pytables + all 31 columns**
> Generated dataframe from C:\Users\erwan\.radisdb\hitemp\H2O-01_01750-02000_HITEMP2010.h5 in 12.56s (1438070 rows, 31 columns, 905.98 MB)
> Generated dataframe from C:\Users\erwan\.radisdb\hitemp\H2O-01_02000-02250_HITEMP2010.h5 in 28.10s (3283544 rows, 31 columns, 2068.63 MB)
> Generated dataframe from C:\Users\erwan\.radisdb\hitemp\H2O-01_02250-02500_HITEMP2010.h5 in 5.19s (608073 rows, 31 columns, 383.09 MB) columns

Total 45.5 s ; 3.2 GB

**Pytables + only 13 columns**
```python
# in loader.py / call to fetch_hitemp(), I add: 
columns=['iso', 'wav', 'int', 'A', 'airbrd', 'selbrd', 'El', 'Tdpair',
                        'Pshft', 'gp', 'gpp', 'Fu', 'jl',
                        ]
```
> Generated dataframe from C:\Users\erwan\.radisdb\hitemp\H2O-01_01750-02000_HITEMP2010.h5 in 12.61s (1438070 rows, 13 columns, 238.72 MB)
> Generated dataframe from C:\Users\erwan\.radisdb\hitemp\H2O-01_02000-02250_HITEMP2010.h5 in 26.65s (3283544 rows, 13 columns, 545.07 MB)
> Generated dataframe from C:\Users\erwan\.radisdb\hitemp\H2O-01_02250-02500_HITEMP2010.h5 in 6.17s (608073 rows, 13 columns, 100.94 MB)

Total 45.2s , 800 MB


**Vaex+ all 31 columns**
> Generated dataframe from ['C:\\Users\\erwan\\.radisdb\\hitemp\\H2O-01_01750-02000_HITEMP2010.hdf5', 'C:\\Users\\erwan\\.radisdb\\hitemp\\H2O-01_02000-02250_HITEMP2010.hdf5', 'C:\\Users\\erwan\\.radisdb\\hitemp\\H2O-01_02250-02500_HITEMP2010.hdf5'] in 20.11s (5329687 rows, 31 columns, 3315.07 MB)

Total 20s + 3.3 GB

**Vaex + 13 columns**
> LDM algorithm used. Defaulting broadening method from <radis.misc.utils.Default object at 0x0000013EE5E56DC0> to Voigt
Generated dataframe from ['C:\\Users\\erwan\\.radisdb\\hitemp\\H2O-01_01750-02000_HITEMP2010.hdf5', 'C:\\Users\\erwan\\.radisdb\\hitemp\\H2O-01_02000-02250_HITEMP2010.hdf5', 'C:\\Users\\erwan\\.radisdb\\hitemp\\H2O-01_02250-02500_HITEMP2010.hdf5'] in 7.45s (5329687 rows, 13 columns, 842.09 MB)

Total 7.5s + 800 MB 


Conclusion : **4x less memory usage** by selecting half the columns;  but same loading time in Pytables. And **6x faster** with Vaex (which, at this point, still generate a Pandas DataFrame in the end)  
@1someshverma we will close this issue when you are done with the full vaex implementation 🚀 

Fix patch:
diff --git a/examples/calculate_spectrum_using_vaex.py b/examples/calculate_spectrum_using_vaex.py
new file mode 100644
index 000000000..ccbd168b1
--- /dev/null
+++ b/examples/calculate_spectrum_using_vaex.py
@@ -0,0 +1,53 @@
+# -*- coding: utf-8 -*-
+"""
+.. _calculate_vaex:
+
+================================================
+Out-of-RAM calculations for very large databases
+================================================
+
+Calculating a spectrum using pandas can be memory-intensive, as it loads data into memory
+for the computation. This can restrict the number of lines that can be processed on a
+device, especially for larger datasets. In such cases, Vaex is recommended as the engine.
+Vaex is a table management tool similar to PANDAS, but it utilizes memory mapping, a zero-
+memory-copy policy, and lazy computations for optimal performance.
+
+This means what whatever your RAM and whatever the size of the dataset, **you won't run out of memory**.
+"""
+
+# To compute spectrum using vaex add an argument engine= "vaex" in calc_spectrum, similarly add output="vaex" in fetch_hitran method to fetch data in vaex format.
+
+from radis import calc_spectrum
+
+s = calc_spectrum(
+    1900,
+    2300,  # cm-1
+    molecule="CO",
+    isotope="1,2,3",
+    pressure=1.01325,  # bar
+    Tgas=700,  # K
+    mole_fraction=0.1,
+    path_length=1,  # cm
+    databank="hitran",  # or 'hitemp', 'geisa', 'exomol'
+    engine="vaex",
+)
+s.apply_slit(0.5, "nm")  # simulate an experimental slit
+s.plot("radiance")
+
+# To compute spectrum using pandas change engine="pandas" (default is also pandas)
+from radis import calc_spectrum
+
+s = calc_spectrum(
+    1900,
+    2300,  # cm-1
+    molecule="CO",
+    isotope="1,2,3",
+    pressure=1.01325,  # bar
+    Tgas=700,  # K
+    mole_fraction=0.1,
+    path_length=1,  # cm
+    databank="hitran",  # or 'hitemp', 'geisa', 'exomol'
+    engine="pandas",
+)
+s.apply_slit(0.5, "nm")  # simulate an experimental slit
+s.plot("radiance")
diff --git a/radis/api/cache_files.py b/radis/api/cache_files.py
index 69636a404..74919b147 100644
--- a/radis/api/cache_files.py
+++ b/radis/api/cache_files.py
@@ -113,7 +113,7 @@ def load_h5_cache_file(
 
     Returns
     -------
-    df: pandas DataFrame, or None
+    df: pandas DataFrame or Vaex Dataframe or None
         None if no cache file was found, or if it was deleted
     """
     # TODO @dev: refactor to use HDF5Manager more
@@ -537,7 +537,7 @@ def check_relevancy(
 
     for k, v in relevant_if_metadata_above.items():
         # Note : check_not_deprecated already tested the existence of each key so we are safe
-        if file_metadata[k] < v:
+        if float(file_metadata[k]) < v:
             raise IrrelevantFileWarning(
                 "Database file {0} irrelevant: {1}={2} [file metadata] < {3} [expected], not loaded".format(
                     file, k, file_metadata[k], v
@@ -545,7 +545,7 @@ def check_relevancy(
             )
     for k, v in relevant_if_metadata_below.items():
         # Note : check_not_deprecated already tested the existence of each key so we are safe
-        if file_metadata[k] > v:
+        if float(file_metadata[k]) > v:
             raise IrrelevantFileWarning(
                 "Database file {0} irrelevant ({1}={2} [file metadata] > {3} [expected]), not loaded".format(
                     file, k, file_metadata[k], v
diff --git a/radis/api/cdsdapi.py b/radis/api/cdsdapi.py
index 3deca6b93..c1ee768c9 100644
--- a/radis/api/cdsdapi.py
+++ b/radis/api/cdsdapi.py
@@ -136,6 +136,7 @@ def cdsd2df(
     load_wavenum_min=None,
     load_wavenum_max=None,
     engine="pytables",
+    output="pandas",
 ):
     """Convert a CDSD-HITEMP [1]_ or CDSD-4000 [2]_ file to a Pandas dataframe.
 
@@ -174,7 +175,7 @@ def cdsd2df(
 
     Returns
     -------
-    df: pandas Dataframe
+    df: pandas Dataframe or Vaex Dataframe
         dataframe containing all lines and parameters
 
     Notes
@@ -280,7 +281,7 @@ def cdsd2df(
 
     # %% Start reading the full file
 
-    df = parse_hitran_file(fname, columns)
+    df = parse_hitran_file(fname, columns, output=output)
 
     # Remove non numerical attributes
     if drop_non_numeric:
diff --git a/radis/api/dbmanager.py b/radis/api/dbmanager.py
index 4fa0966d0..b27664c7e 100644
--- a/radis/api/dbmanager.py
+++ b/radis/api/dbmanager.py
@@ -516,7 +516,7 @@ def add_column(self, df, key, value):
 
             if isinstance(df, vaex.dataframe.DataFrameLocal):
                 # see https://github.com/vaexio/vaex/pull/1570
-                df[key] = vaex.vconstant(float(value), length=len(df))
+                df[key] = vaex.vconstant(float(value), length=df.length_unfiltered())
             else:
                 df[key] = value
         else:
diff --git a/radis/api/hdf5.py b/radis/api/hdf5.py
index f638d8e1a..0c2049832 100644
--- a/radis/api/hdf5.py
+++ b/radis/api/hdf5.py
@@ -367,7 +367,11 @@ def load(
             raise NotImplementedError(output)
 
         if output == "vaex":
-            df = df.extract()  # return DataFrame containing only the filtered rows
+            # df = df.extract()  # Not required # return DataFrame containing only the filtered rows
+            if columns:  # load only these columns (if they exist)
+                columns = [c for c in columns if c in df.columns]
+            if columns is not None:
+                df = df[columns]
 
         return df
 
diff --git a/radis/api/hitranapi.py b/radis/api/hitranapi.py
index 7884e215c..b39cba83e 100644
--- a/radis/api/hitranapi.py
+++ b/radis/api/hitranapi.py
@@ -115,16 +115,48 @@
 }
 
 
-def cast_to_int64_with_missing_values(dg, keys):
+def cast_to_int64_with_missing_values(dg, keys, dataframe_type="pandas"):
     """replace missing values of int64 columns with -1"""
     for c in keys:
-        if dg.dtypes[c] != int64:
-            dg[c].replace(
-                r"^\s+$", -1, regex=True, inplace=True
-            )  # replace empty strings by -1, e.g. HCN
-            # Warning: -1 may be a valid non-equilibrium quantum number for some
-            # molecules, e.g. H2O, see https://github.com/radis/radis/issues/280#issuecomment-896120510
-            dg[c] = dg[c].fillna(-1).astype(int64)  # replace nans with -1
+        if dataframe_type == "pandas":
+            dtype = dg.dtypes[c]
+        elif dataframe_type == "vaex":
+            dtype = dg.data_type(c)
+
+        if dtype != int64:
+            if dataframe_type == "pandas":
+                dg[c].replace(
+                    r"^\s+$", -1, regex=True, inplace=True
+                )  # replace empty strings by -1, e.g. HCN
+                # Warning: -1 may be a valid non-equilibirum quantum number for some
+                # molecules, e.g. H2O, see https://github.com/radis/radis/issues/280#issuecomment-896120510
+                dg[c] = dg[c].fillna(-1).astype(int64)  # replace nans with -1
+            elif dataframe_type == "vaex":
+                dg[c] = dg[c].str.replace(
+                    r"^\s+$", "-1", regex=True
+                )  # replace empty strings by -1, e.g. HCN
+                # Warning: -1 may be a valid non-equilibirum quantum number for some
+                # molecules, e.g. H2O, see https://github.com/radis/radis/issues/280#issuecomment-896120510
+                # dg[c] = dg[c].apply(lambda x : int(x))
+            else:
+                raise NotImplementedError(dataframe_type)
+
+
+def cast_all_to_int64(df, keys):
+    """converts string to int64"""
+    for key in keys:
+        df[key] = df[key].apply(lambda x: int(x) if (x.strip()).isdigit() else x)
+
+
+def extract_columns(df, extracted_values, columns):
+    """extracts column from extracted_values"""
+    for column in columns:
+        df[column] = extracted_values.apply(lambda x: x.get(column))
+
+    cast_all_to_int64(df, columns)
+
+    for column in columns:
+        df[column] = df.evaluate(df[column])
 
 
 def hit2df(
@@ -135,6 +167,7 @@ def hit2df(
     load_wavenum_min=None,
     load_wavenum_max=None,
     engine="pytables",
+    output="pandas",
     parse_quanta=True,
 ):
     """Convert a HITRAN/HITEMP [1]_ file to a Pandas dataframe
@@ -166,10 +199,12 @@ def hit2df(
     parse_quanta: bool
         if ``True``, parse local & global quanta (required to identify lines
         for non-LTE calculations ; but sometimes lines are not labelled.)
+    output : str
+        output format of data as pandas Dataformat or vaex Dataformat
 
     Returns
     -------
-    df: pandas Dataframe
+    df: pandas Dataframe or Vaex Dataframe
         dataframe containing all lines and parameters
 
 
@@ -239,11 +274,11 @@ def hit2df(
             + "instead of an HITRAN file"
         ) from err
 
-    df = parse_hitran_file(fname, columns)
-
+    df = parse_hitran_file(fname, columns, output=output)
     df = post_process_hitran_data(
         df,
         molecule=mol,
+        dataframe_type=output,
         parse_quanta=parse_quanta,
     )
     # cached file mode but cached file doesn't exist yet (else we had returned)
@@ -295,6 +330,7 @@ def post_process_hitran_data(
     drop_non_numeric=True,
     parse_quanta=True,
     add_HITRAN_uncertainty_code=False,
+    dataframe_type="pandas",
 ):
     """Parsing non-equilibrium parameters in HITRAN/HITEMP [1]_ file to and return final Pandas Dataframe
 
@@ -318,10 +354,12 @@ def post_process_hitran_data(
         for non-LTE calculations ; but sometimes lines are not labelled.)
     add_HITRAN_uncertainty_code: bool
         if ``True``, a column which contains HITRAN uncertainty code is converted to integer and not dropped.
+    engine: str
+        pandas or vaex
 
     Returns
     -------
-    df: pandas Dataframe
+    df: pandas Dataframe or vaex Dataframe
         dataframe containing all lines and parameters
 
 
@@ -354,7 +392,12 @@ def post_process_hitran_data(
     elif nmol != 1:
         # Crash, give explicitly error messages
         try:
-            secondline = df.iloc[1]
+            if dataframe_type == "pandas":
+                secondline = df.iloc[1]
+            elif dataframe_type == "vaex":
+                secondline = df.values[1]
+            else:
+                raise NotImplementedError(dataframe_type)
         except IndexError:
             secondline = ""
         raise ValueError(
@@ -371,7 +414,9 @@ def post_process_hitran_data(
     if parse_quanta:
         # Add local quanta attributes, based on the HITRAN group
         try:
-            df = parse_local_quanta(df, molecule, verbose=verbose)
+            df = parse_local_quanta(
+                df, molecule, verbose=verbose, dataframe_type=dataframe_type
+            )
         except ValueError as err:
             # Empty strings (unlabelled lines) have been reported for HITEMP2010-H2O.
             # In this case, do not parse (makes non-equilibrium calculations impossible).
@@ -385,7 +430,9 @@ def post_process_hitran_data(
 
         # Add global quanta attributes, based on the HITRAN class
         try:
-            df = parse_global_quanta(df, molecule, verbose=verbose)
+            df = parse_global_quanta(
+                df, molecule, verbose=verbose, dataframe_type=dataframe_type
+            )
         except ValueError as err:
             # Empty strings (unlabelled lines) have been reported for HITEMP2010-H2O.
             # In this case, do not parse (makes non-equilibrium calculations impossible).
@@ -411,7 +458,7 @@ def post_process_hitran_data(
 # %% Hitran global quanta classes
 
 
-def _parse_HITRAN_class1(df, verbose=True):
+def _parse_HITRAN_class1(df, verbose=True, dataframe_type="pandas"):
     r"""Diatomic molecules: CO, HF, HCl, HBr, HI, N2, NO+
 
 
@@ -419,7 +466,12 @@ def _parse_HITRAN_class1(df, verbose=True):
     ----------
     df: pandas Dataframe
         lines read from a HITRAN-like database
+    dataframe_type : str
+        pandas or vaex
 
+    Returns
+    -------
+        pandas Dataframe or Vaex Dataframe
 
     Notes
     -----
@@ -436,19 +488,51 @@ def _parse_HITRAN_class1(df, verbose=True):
 
     """
 
-    # 1. Parse
-    dgu = df["globu"].astype(str).str.extract(r"[ ]{13}(?P<vu>[\d ]{2})", expand=True)
-    dgl = df["globl"].astype(str).str.extract(r"[ ]{13}(?P<vl>[\d ]{2})", expand=True)
+    if dataframe_type == "vaex":
+        # 1. Parse
+
+        extracted_values = df["globu"].str.extract_regex(
+            pattern=r"[ ]{13}(?P<vu>[\d ]{2})"
+        )
+
+        extract_columns(df, extracted_values, ["vu"])
+
+        extracted_values = df["globl"].str.extract_regex(
+            pattern=r"[ ]{13}(?P<vl>[\d ]{2})"
+        )
 
-    # 2. Convert to numeric
-    cast_to_int64_with_missing_values(dgu, ["vu"])
-    cast_to_int64_with_missing_values(dgl, ["vl"])
+        extract_columns(df, extracted_values, ["vl"])
 
-    # 3. Clean
-    del df["globu"]
-    del df["globl"]
+        # # 2. Convert to numeric
+        # cast_to_int64_with_missing_values(
+        #     df, ["vu", "vl"], dataframe_type=dataframe_type
+        # )
 
-    return pd.concat([df, dgu, dgl], axis=1)
+        # 3. Clean
+        df.drop("globu", inplace=True)
+        df.drop("globl", inplace=True)
+
+        return df
+    elif dataframe_type == "pandas":
+        # 1. Parse
+        dgu = (
+            df["globu"].astype(str).str.extract(r"[ ]{13}(?P<vu>[\d ]{2})", expand=True)
+        )
+        dgl = (
+            df["globl"].astype(str).str.extract(r"[ ]{13}(?P<vl>[\d ]{2})", expand=True)
+        )
+
+        # 2. Convert to numeric
+        cast_to_int64_with_missing_values(dgu, ["vu"], dataframe_type=dataframe_type)
+        cast_to_int64_with_missing_values(dgl, ["vl"], dataframe_type=dataframe_type)
+
+        # 3. Clean
+        del df["globu"]
+        del df["globl"]
+
+        return pd.concat([df, dgu, dgl], axis=1)
+    else:
+        raise NotImplementedError(dataframe_type)
 
 
 def _parse_HITRAN_class2(df, verbose=True):
@@ -515,7 +599,7 @@ def _parse_HITRAN_class3(df, verbose=True):
     return df
 
 
-def _parse_HITRAN_class4(df, verbose=True):
+def _parse_HITRAN_class4(df, verbose=True, dataframe_type="pandas"):
     r"""Parse linear triatomic class in HITRAN [1]_: N2O, OCS, HCN
 
     Parameters
@@ -523,6 +607,12 @@ def _parse_HITRAN_class4(df, verbose=True):
 
     df: pandas Dataframe
         lines read from a HITRAN-like database
+    dataframe_type : str
+        pandas or vaex
+
+    Returns
+    -------
+        pandas Dataframe or Vaex Dataframe
 
     Notes
     -----
@@ -540,37 +630,64 @@ def _parse_HITRAN_class4(df, verbose=True):
     .. [1] `Table 3 of Rothman et al. HITRAN 2004 <https://www.cfa.harvard.edu/hitran/Download/HITRAN04paper.pdf>`__
 
     """
+    if dataframe_type == "pandas":
+        # 1. Parse
+        dgu = (
+            df["globu"]
+            .astype(str)
+            .str.extract(
+                r"[ ]{7}(?P<v1u>[\d ]{2})(?P<v2u>[\d ]{2})(?P<l2u>[\d ]{2})(?P<v3u>[\d ]{2})",
+                expand=True,
+            )
+        )
+        dgl = (
+            df["globl"]
+            .astype(str)
+            .str.extract(
+                r"[ ]{7}(?P<v1l>[\d ]{2})(?P<v2l>[\d ]{2})(?P<l2l>[\d ]{2})(?P<v3l>[\d ]{2})",
+                expand=True,
+            )
+        )
 
-    # 1. Parse
-    dgu = (
-        df["globu"]
-        .astype(str)
-        .str.extract(
-            r"[ ]{7}(?P<v1u>[\d ]{2})(?P<v2u>[\d ]{2})(?P<l2u>[\d ]{2})(?P<v3u>[\d ]{2})",
-            expand=True,
+        # 2. Convert to numeric
+        cast_to_int64_with_missing_values(dgu, ["v1u", "v2u", "l2u", "v3u"])
+        cast_to_int64_with_missing_values(dgl, ["v1l", "v2l", "l2l", "v3l"])
+
+        # 3. Clean
+        del df["globu"]
+        del df["globl"]
+
+        return pd.concat([df, dgu, dgl], axis=1)
+    elif dataframe_type == "vaex":
+        # 1. Parse
+        extracted_values = df["globu"].str.extract_regex(
+            pattern=r"[ ]{7}(?P<v1u>[\d ]{2})(?P<v2u>[\d ]{2})(?P<l2u>[\d ]{2})(?P<v3u>[\d ]{2})"
         )
-    )
-    dgl = (
-        df["globl"]
-        .astype(str)
-        .str.extract(
-            r"[ ]{7}(?P<v1l>[\d ]{2})(?P<v2l>[\d ]{2})(?P<l2l>[\d ]{2})(?P<v3l>[\d ]{2})",
-            expand=True,
+
+        extract_columns(df, extracted_values, ["v1u", "v2u", "l2u", "v3u"])
+
+        extracted_values = df["globl"].str.extract_regex(
+            pattern=r"[ ]{7}(?P<v1l>[\d ]{2})(?P<v2l>[\d ]{2})(?P<l2l>[\d ]{2})(?P<v3l>[\d ]{2})"
         )
-    )
 
-    # 2. Convert to numeric
-    cast_to_int64_with_missing_values(dgu, ["v1u", "v2u", "l2u", "v3u"])
-    cast_to_int64_with_missing_values(dgl, ["v1l", "v2l", "l2l", "v3l"])
+        extract_columns(df, extracted_values, ["v1l", "v2l", "l2l", "v3l"])
+
+        # # 2. Convert to numeric
+        # cast_to_int64_with_missing_values(
+        #     df,
+        #     ["v1u", "v2u", "l2u", "v3u", "v1l", "v2l", "l2l", "v3l"],
+        #     dataframe_type=dataframe_type,
+        # )
 
-    # 3. Clean
-    del df["globu"]
-    del df["globl"]
+        # 3. Clean
+        df.drop("globu", inplace=True)
 
-    return pd.concat([df, dgu, dgl], axis=1)
+        return df
+    else:
+        raise NotImplementedError(dataframe_type)
 
 
-def _parse_HITRAN_class5(df, verbose=True):
+def _parse_HITRAN_class5(df, verbose=True, dataframe_type="pandas"):
     r"""Parse linear triatomic with large Fermi resonance in HITRAN [1]_: CO2
 
     Parameters
@@ -578,6 +695,12 @@ def _parse_HITRAN_class5(df, verbose=True):
 
     df: pandas Dataframe
         lines read from a HITRAN-like database
+    dataframe_type : str
+        pandas or vaex
+
+    Returns
+    -------
+        pandas Dataframe or Vaex Dataframe
 
     Notes
     -----
@@ -596,36 +719,63 @@ def _parse_HITRAN_class5(df, verbose=True):
 
     """
 
-    # 1. Parse
-    dgu = (
-        df["globu"]
-        .astype(str)
-        .str.extract(
-            r"[ ]{6}(?P<v1u>[\d ]{2})(?P<v2u>[\d ]{2})(?P<l2u>[\d ]{2})(?P<v3u>[\d ]{2})(?P<ru>\d)",
-            expand=True,
+    if dataframe_type == "pandas":
+        # 1. Parse
+        dgu = (
+            df["globu"]
+            .astype(str)
+            .str.extract(
+                r"[ ]{6}(?P<v1u>[\d ]{2})(?P<v2u>[\d ]{2})(?P<l2u>[\d ]{2})(?P<v3u>[\d ]{2})(?P<ru>\d)",
+                expand=True,
+            )
         )
-    )
-    dgl = (
-        df["globl"]
-        .astype(str)
-        .str.extract(
-            r"[ ]{6}(?P<v1l>[\d ]{2})(?P<v2l>[\d ]{2})(?P<l2l>[\d ]{2})(?P<v3l>[\d ]{2})(?P<rl>\d)",
-            expand=True,
+        dgl = (
+            df["globl"]
+            .astype(str)
+            .str.extract(
+                r"[ ]{6}(?P<v1l>[\d ]{2})(?P<v2l>[\d ]{2})(?P<l2l>[\d ]{2})(?P<v3l>[\d ]{2})(?P<rl>\d)",
+                expand=True,
+            )
+        )
+
+        # 2. Convert to numeric
+        cast_to_int64_with_missing_values(dgu, ["v1u", "v2u", "l2u", "v3u", "ru"])
+        cast_to_int64_with_missing_values(dgl, ["v1l", "v2l", "l2l", "v3l", "rl"])
+
+        # 3. Clean
+        del df["globu"]
+        del df["globl"]
+
+        return pd.concat([df, dgu, dgl], axis=1)
+    elif dataframe_type == "vaex":
+        # 1. Parse
+        extracted_values = df["globu"].str.extract_regex(
+            pattern=r"[ ]{6}(?P<v1u>[\d ]{2})(?P<v2u>[\d ]{2})(?P<l2u>[\d ]{2})(?P<v3u>[\d ]{2})(?P<ru>\d)"
+        )
+
+        extract_columns(df, extracted_values, ["v1u", "v2u", "l2u", "v3u", "ru"])
+
+        extracted_values = df["globl"].str.extract_regex(
+            pattern=r"[ ]{6}(?P<v1l>[\d ]{2})(?P<v2l>[\d ]{2})(?P<l2l>[\d ]{2})(?P<v3l>[\d ]{2})(?P<rl>\d)"
         )
-    )
 
-    # 2. Convert to numeric
-    cast_to_int64_with_missing_values(dgu, ["v1u", "v2u", "l2u", "v3u", "ru"])
-    cast_to_int64_with_missing_values(dgl, ["v1l", "v2l", "l2l", "v3l", "rl"])
+        extract_columns(df, extracted_values, ["v1l", "v2l", "l2l", "v3l", "rl"])
 
-    # 3. Clean
-    del df["globu"]
-    del df["globl"]
+        # # 2. Convert to numeric
+        # cast_to_int64_with_missing_values(
+        #     df,
+        #     ["v1u", "v2u", "l2u", "v3u", "ru", "v1l", "v2l", "l2l", "v3l", "rl"],
+        #     dataframe_type=dataframe_type,
+        # )
 
-    return pd.concat([df, dgu, dgl], axis=1)
+        # 3. Clean
+        df.drop("globu", inplace=True)
+        df.drop("globl", inplace=True)
 
+        return df
 
-def _parse_HITRAN_class6(df, verbose=True):
+
+def _parse_HITRAN_class6(df, verbose=True, dataframe_type="pandas"):
     r"""Parse non-linear triatomic in HITRAN [1]_: H2O, O3, SO2, NO2, HOCl, H2S, HO2, HOBr
 
     Parameters
@@ -633,7 +783,12 @@ def _parse_HITRAN_class6(df, verbose=True):
 
     df: pandas Dataframe
         lines read from a HITRAN-like database
+    dataframe_type : str
+        pandas or vaex
 
+    Returns
+    -------
+        pandas Dataframe or Vaex Dataframe
     Notes
     -----
 
@@ -650,45 +805,74 @@ def _parse_HITRAN_class6(df, verbose=True):
     .. [1] `Table 3 of Rothman et al. HITRAN 2004 <https://www.cfa.harvard.edu/hitran/Download/HITRAN04paper.pdf>`__
 
     """
+    if dataframe_type == "pandas":
+        # 1. Parse
+        dgu = (
+            df["globu"]
+            .astype(str)
+            .str.extract(
+                #        '[ ]{9}(?P<v1u>[\d ]{2})(?P<v2u>[\d ]{2})(?P<v3u>[\d ]{2})',
+                r"[ ]{9}(?P<v1u>[\-\d ]{2})(?P<v2u>[\-\d ]{2})(?P<v3u>[\-\d ]{2})",
+                expand=True,
+            )
+        )
+        dgl = (
+            df["globl"]
+            .astype(str)
+            .str.extract(
+                #        '[ ]{9}(?P<v1l>[\d ]{2})(?P<v2l>[\d ]{2})(?P<v3l>[\d ]{2})',
+                r"[ ]{9}(?P<v1l>[\-\d ]{2})(?P<v2l>[\-\d ]{2})(?P<v3l>[\-\d ]{2})",
+                expand=True,
+            )
+        )
+        # ... note @EP: in HITRAN H2O files, for iso=2, vibrational levels are
+        # ... somehow negative. The regex above is adapted to catch negation signs with \-
+
+        # 2. Convert to numeric
+        cast_to_int64_with_missing_values(
+            dgu,
+            [
+                "v1u",
+                "v2u",
+                "v3u",
+            ],
+        )
+        cast_to_int64_with_missing_values(dgl, ["v1l", "v2l", "v3l"])
 
-    # 1. Parse
-    dgu = (
-        df["globu"]
-        .astype(str)
-        .str.extract(
-            #        '[ ]{9}(?P<v1u>[\d ]{2})(?P<v2u>[\d ]{2})(?P<v3u>[\d ]{2})',
-            r"[ ]{9}(?P<v1u>[\-\d ]{2})(?P<v2u>[\-\d ]{2})(?P<v3u>[\-\d ]{2})",
-            expand=True,
+        # 3. Clean
+        del df["globu"]
+        del df["globl"]
+
+        return pd.concat([df, dgu, dgl], axis=1)
+    elif dataframe_type == "vaex":
+        # 1. Parse
+
+        # ... note @EP: in HITRAN H2O files, for iso=2, vibrational levels are
+        # ... somehow negative. The regex below is adapted to catch negation signs with \-
+        extracted_values = df["globu"].str.extract_regex(
+            pattern=r"[ ]{9}(?P<v1u>[\-\d ]{2})(?P<v2u>[\-\d ]{2})(?P<v3u>[\-\d ]{2})"
         )
-    )
-    dgl = (
-        df["globl"]
-        .astype(str)
-        .str.extract(
-            #        '[ ]{9}(?P<v1l>[\d ]{2})(?P<v2l>[\d ]{2})(?P<v3l>[\d ]{2})',
-            r"[ ]{9}(?P<v1l>[\-\d ]{2})(?P<v2l>[\-\d ]{2})(?P<v3l>[\-\d ]{2})",
-            expand=True,
+
+        extract_columns(df, extracted_values, ["v1u", "v2u", "v3u"])
+
+        extracted_values = df["globl"].str.extract_regex(
+            pattern=r"[ ]{9}(?P<v1l>[\-\d ]{2})(?P<v2l>[\-\d ]{2})(?P<v3l>[\-\d ]{2})"
         )
-    )
-    # ... note @EP: in HITRAN H2O files, for iso=2, vibrational levels are
-    # ... somehow negative. The regex above is adapted to catch negation signs with \-
-
-    # 2. Convert to numeric
-    cast_to_int64_with_missing_values(
-        dgu,
-        [
-            "v1u",
-            "v2u",
-            "v3u",
-        ],
-    )
-    cast_to_int64_with_missing_values(dgl, ["v1l", "v2l", "v3l"])
 
-    # 3. Clean
-    del df["globu"]
-    del df["globl"]
+        extract_columns(df, extracted_values, ["v1l", "v2l", "v3l"])
 
-    return pd.concat([df, dgu, dgl], axis=1)
+        # # 2. Convert to numeric
+        # cast_to_int64_with_missing_values(
+        #     df, ["v1l", "v2l", "v3l", "v1u", "v2u", "v3u"]
+        # )
+
+        # 3. Clean
+        df.drop("globu", inplace=True)
+        df.drop("globl", inplace=True)
+
+        return df
+    else:
+        raise NotImplementedError(dataframe_type)
 
 
 def _parse_HITRAN_class7(df, verbose=True):
@@ -816,7 +1000,7 @@ def _parse_HITRAN_class10(df, verbose=True):
 # %% HITRAN Local quanta
 
 
-def _parse_HITRAN_group1(df, verbose=True):
+def _parse_HITRAN_group1(df, verbose=True, dataframe_type="pandas"):
     r"""Parse asymmetric rotors (:py:attr:`~radis.db.classes.HITRAN_GROUP1` ):
     H2O, O3, SO2, NO2, HNO3, H2CO, HOCl, H2O2, COF2, H2S, HO2, HCOOH, ClONO2, HOBr, C2H4
 
@@ -825,6 +1009,12 @@ def _parse_HITRAN_group1(df, verbose=True):
 
     df: pandas Dataframe
         lines read from a HITRAN-like database
+    dataframe_type : str
+        pandas or vaex
+
+    Returns
+    -------
+        pandas Dataframe or Vaex Dataframe
 
 
     Notes
@@ -841,47 +1031,87 @@ def _parse_HITRAN_group1(df, verbose=True):
 
     """
 
-    # 1. Parse
+    if dataframe_type == "vaex":
+        # 1. Parse
 
-    # Ref [1] : locu
-    # --------------
-    # J'  | Ka' | Kc' | F'  | Sym'
-    # I3  | I3  | I3  | A5  | A1
-    dgu = (
-        df["locu"]
-        .astype(str)
-        .str.extract(
-            r"(?P<ju>[\d ]{3})(?P<Kau>[\-\d ]{3})(?P<Kcu>[\-\d ]{3})(?P<Fu>.{5})(?P<symu>.)",
-            expand=True,
+        # Ref [1] : locu
+        # --------------
+        # J'  | Ka' | Kc' | F'  | Sym'
+        # I3  | I3  | I3  | A5  | A1
+
+        extracted_values = df["locu"].str.extract_regex(
+            pattern=r"(?P<ju>[\d ]{3})(?P<Kau>[\-\d ]{3})(?P<Kcu>[\-\d ]{3})(?P<Fu>.{5})(?P<symu>.)"
         )
-    )
-    # Ref [1] : locl
-    # --------------
-    # J'' | Ka''| Kc''| F'' | Sym''
-    # I3  | I3  | I3  | A5  | A1
-    dgl = (
-        df["locl"]
-        .astype(str)
-        .str.extract(
-            r"(?P<jl>[\d ]{3})(?P<Kal>[\-\d ]{3})(?P<Kcl>[\-\d ]{3})(?P<Fl>.{5})(?P<syml>.)",
-            expand=True,
+
+        extract_columns(df, extracted_values, ["ju", "Kau", "Kcu", "Fu", "symu"])
+
+        # Ref [1] : locl
+        # --------------
+        # J'' | Ka''| Kc''| F'' | Sym''
+        # I3  | I3  | I3  | A5  | A1
+
+        extracted_values = df["locl"].str.extract_regex(
+            pattern=r"(?P<jl>[\d ]{3})(?P<Kal>[\-\d ]{3})(?P<Kcl>[\-\d ]{3})(?P<Fl>.{5})(?P<syml>.)"
         )
-    )
-    # ... note @EP: in HITRAN H2O files, for iso=2, the Kau, Kcu can somehow
-    # ... be negative. The regex above is adapted to catch negation signs with \-
 
-    # 2. Convert to numeric
-    cast_to_int64_with_missing_values(dgu, ["ju", "Kau", "Kcu"])
-    cast_to_int64_with_missing_values(dgl, ["jl", "Kal", "Kcl"])
+        extract_columns(df, extracted_values, ["jl", "Kal", "Kcl", "Fl", "syml"])
+
+        # ... note @EP: in HITRAN H2O files, for iso=2, the Kau, Kcu can somehow
+        # ... be negative. The regex above is adapted to catch negation signs with \-
+        # # 2. Convert to numeric
+        # cast_to_int64_with_missing_values(
+        #     df, ["ju", "Kau", "Kcu", "Fu", "symu", "jl", "Kal", "Kcl", "Fl", "syml"]
+        # )
+
+        # 3. Clean
+        df.drop("locu", inplace=True)
+        df.drop("locl", inplace=True)
+
+        return df
+    elif dataframe_type == "pandas":
+        # 1. Parse
+
+        # Ref [1] : locu
+        # --------------
+        # J'  | Ka' | Kc' | F'  | Sym'
+        # I3  | I3  | I3  | A5  | A1
+
+        dgu = (
+            df["locu"]
+            .astype(str)
+            .str.extract(
+                r"(?P<ju>[\d ]{3})(?P<Kau>[\-\d ]{3})(?P<Kcu>[\-\d ]{3})(?P<Fu>.{5})(?P<symu>.)",
+                expand=True,
+            )
+        )
+        # Ref [1] : locl
+        # --------------
+        # J'' | Ka''| Kc''| F'' | Sym''
+        # I3  | I3  | I3  | A5  | A1
+        dgl = (
+            df["locl"]
+            .astype(str)
+            .str.extract(
+                r"(?P<jl>[\d ]{3})(?P<Kal>[\-\d ]{3})(?P<Kcl>[\-\d ]{3})(?P<Fl>.{5})(?P<syml>.)",
+                expand=True,
+            )
+        )
+        # ... note @EP: in HITRAN H2O files, for iso=2, the Kau, Kcu can somehow
+        # ... be negative. The regex above is adapted to catch negation signs with \-
+        # 2. Convert to numeric
+        cast_to_int64_with_missing_values(dgu, ["ju", "Kau", "Kcu"])
+        cast_to_int64_with_missing_values(dgl, ["jl", "Kal", "Kcl"])
 
-    # 3. Clean
-    del df["locu"]
-    del df["locl"]
+        # 3. Clean
+        del df["locu"]
+        del df["locl"]
 
-    return pd.concat([df, dgu, dgl], axis=1)
+        return pd.concat([df, dgu, dgl], axis=1)
+    else:
+        raise NotImplementedError(dataframe_type)
 
 
-def _parse_HITRAN_group2(df, verbose=True):
+def _parse_HITRAN_group2(df, verbose=True, dataframe_type="pandas"):
     r"""Parse diatomic and linear molecules (:py:attr:`~radis.db.classes.HITRAN_GROUP2` ):
     CO2, N2O, CO, HF, HCl, HBr, HI, OCS, N2, HCN, C2H2, NO+
 
@@ -890,6 +1120,12 @@ def _parse_HITRAN_group2(df, verbose=True):
 
     df: pandas Dataframe
         lines read from a HITRAN-like database
+    dataframe_type : str
+        pandas or vaex
+
+    Returns
+    -------
+        pandas Dataframe or Vaex Dataframe
 
 
     Notes
@@ -912,30 +1148,61 @@ def _parse_HITRAN_group2(df, verbose=True):
     # --------------
     #     | F'  |
     # 10X | A5  |
-    dgu = df["locu"].astype(str).str.extract(r"[ ]{10}(?P<Fu>.{5})", expand=True)
-    # Ref [1] : locl
-    # --------------
-    #     | Br  | J'' | Sym''| F'' |
-    # 5X  | A1  | I3  | A1   | A5  |
-    dgl = (
-        df["locl"]
-        .astype(str)
-        .str.extract(
-            r"[ ]{5}(?P<branch>[\S]{1})(?P<jl>[\d ]{3})(?P<syml>.)(?P<Fl>.{5})",
-            expand=True,
+    if dataframe_type == "pandas":
+        dgu = df["locu"].astype(str).str.extract(r"[ ]{10}(?P<Fu>.{5})", expand=True)
+
+        # Ref [1] : locl
+        # --------------
+        #     | Br  | J'' | Sym''| F'' |
+        # 5X  | A1  | I3  | A1   | A5  |
+        dgl = (
+            df["locl"]
+            .astype(str)
+            .str.extract(
+                r"[ ]{5}(?P<branch>[\S]{1})(?P<jl>[\d ]{3})(?P<syml>.)(?P<Fl>.{5})",
+                expand=True,
+            )
         )
-    )
 
-    # 2. Convert to numeric
+        # 2. Convert to numeric
+
+        # dgl['jl'] = dgl.jl.apply(pd.to_numeric)
+        cast_to_int64_with_missing_values(dgl, ["jl"], dataframe_type=dataframe_type)
+
+        # 3. Clean
+        del df["locu"]
+        del df["locl"]
 
-    # dgl['jl'] = dgl.jl.apply(pd.to_numeric)
-    cast_to_int64_with_missing_values(dgl, ["jl"])
+        return pd.concat([df, dgu, dgl], axis=1)
+    elif dataframe_type == "vaex":
+        # 1.Parse
+        extracted_values = df["locu"].str.extract_regex(pattern=r"[ ]{10}(?P<Fu>.{5})")
+        extract_columns(df, extracted_values, ["Fu"])
 
-    # 3. Clean
-    del df["locu"]
-    del df["locl"]
+        # dgu = df["locu"].astype(str).str.extract(r"[ ]{10}(?P<Fu>.{5})", expand=True)
+        # Ref [1] : locl
+        # --------------
+        #     | Br  | J'' | Sym''| F'' |
+        # 5X  | A1  | I3  | A1   | A5  |
 
-    return pd.concat([df, dgu, dgl], axis=1)
+        extracted_values = df["locl"].str.extract_regex(
+            pattern=r"[ ]{5}(?P<branch>[\S]{1})(?P<jl>[\d ]{3})(?P<syml>.)(?P<Fl>.{5})"
+        )
+
+        extract_columns(df, extracted_values, ["branch", "jl", "syml", "Fl"])
+
+        # # 2. Convert to numeric
+        # cast_to_int64_with_missing_values(
+        #     df, ["Fu", "branch", "jl", "syml", "Fl"], dataframe_type=dataframe_type
+        # )
+
+        # 3. Clean
+        df.drop("locu", inplace=True)
+        df.drop("locl", inplace=True)
+
+        return df
+    else:
+        raise NotImplementedError(dataframe_type)
 
 
 #                # 10X included in Fu
@@ -1074,7 +1341,7 @@ def _parse_HITRAN_group6(df, verbose=True):
 # %% Reading function
 
 
-def parse_local_quanta(df, mol, verbose=True):
+def parse_local_quanta(df, mol, verbose=True, dataframe_type="pandas"):
     r"""
     Parameters
     ----------
@@ -1086,9 +1353,9 @@ def parse_local_quanta(df, mol, verbose=True):
     """
 
     if mol in HITRAN_GROUP1:
-        df = _parse_HITRAN_group1(df, verbose=verbose)
+        df = _parse_HITRAN_group1(df, verbose=verbose, dataframe_type=dataframe_type)
     elif mol in HITRAN_GROUP2:
-        df = _parse_HITRAN_group2(df, verbose=verbose)
+        df = _parse_HITRAN_group2(df, verbose=verbose, dataframe_type=dataframe_type)
     elif mol in HITRAN_GROUP3:
         df = _parse_HITRAN_group3(df, verbose=verbose)
     elif mol in HITRAN_GROUP4:
@@ -1105,7 +1372,7 @@ def parse_local_quanta(df, mol, verbose=True):
     return df
 
 
-def parse_global_quanta(df, mol, verbose=True):
+def parse_global_quanta(df, mol, verbose=True, dataframe_type="pandas"):
     r"""
 
     Parameters
@@ -1118,17 +1385,17 @@ def parse_global_quanta(df, mol, verbose=True):
     """
 
     if mol in HITRAN_CLASS1:
-        df = _parse_HITRAN_class1(df, verbose=verbose)
+        df = _parse_HITRAN_class1(df, verbose=verbose, dataframe_type=dataframe_type)
     elif mol in HITRAN_CLASS2:
         df = _parse_HITRAN_class2(df, verbose=verbose)
     elif mol in HITRAN_CLASS3:
         df = _parse_HITRAN_class3(df, verbose=verbose)
     elif mol in HITRAN_CLASS4:
-        df = _parse_HITRAN_class4(df, verbose=verbose)
+        df = _parse_HITRAN_class4(df, verbose=verbose, dataframe_type=dataframe_type)
     elif mol in HITRAN_CLASS5:
-        df = _parse_HITRAN_class5(df, verbose=verbose)
+        df = _parse_HITRAN_class5(df, verbose=verbose, dataframe_type=dataframe_type)
     elif mol in HITRAN_CLASS6:
-        df = _parse_HITRAN_class6(df, verbose=verbose)
+        df = _parse_HITRAN_class6(df, verbose=verbose, dataframe_type=dataframe_type)
     elif mol in HITRAN_CLASS7:
         df = _parse_HITRAN_class7(df, verbose=verbose)
     elif mol in HITRAN_CLASS8:
diff --git a/radis/api/tools.py b/radis/api/tools.py
index 8275a74b3..72a6b18bd 100644
--- a/radis/api/tools.py
+++ b/radis/api/tools.py
@@ -9,11 +9,12 @@
 
 import numpy as np
 import pandas as pd
+import vaex
 
 from ..misc.warning import PerformanceWarning
 
 
-def parse_hitran_file(fname, columns, count=-1):
+def parse_hitran_file(fname, columns, count=-1, output="pandas"):
     """Parse a file under HITRAN ``par`` format. Parsing is done in binary
     format with :py:func:`numpy.fromfile` so it's as fast as possible.
 
@@ -29,9 +30,12 @@ def parse_hitran_file(fname, columns, count=-1):
     count: int
         number of lines to read. If ``-1`` reads all file.
 
+    output : str
+        specifies the output type
+
     Returns
     -------
-    df: pandas DataFrame
+    df: pandas DataFrame or Vaex Dataframe
         dataframe with lines.
 
     See Also
@@ -55,7 +59,13 @@ def parse_hitran_file(fname, columns, count=-1):
     data = _read_hitran_file(fname, columns, count, linereturnformat)
 
     # Return a Pandas dataframe
-    return _ndarray2df(data, columns, linereturnformat)
+    df = _ndarray2df(data, columns, linereturnformat)
+
+    import vaex
+
+    if output == "vaex":
+        df = vaex.from_pandas(df)
+    return df
 
 
 def _get_linereturnformat(data, columns, fname=""):
@@ -209,16 +219,22 @@ def _cast_to_dtype(data, dtype):
 
 
 def drop_object_format_columns(df, verbose=True):
-    """Remove 'object' columns in a pandas DataFrame.
+    """Remove 'object' columns in a pandas DataFrame. If Vaex Dataframe, there is by construction no object in the Dataframe
 
     They are not useful to us at this time, and they slow down all
     operations (as they are converted to 'object' in pandas DataFrame).
     If you want to keep them, better convert them to some numeric values
     """
 
+    df_type = type(df)
     objects = [k for k, v in df.dtypes.items() if v == object]
-    for k in objects:
-        del df[k]
+    if df_type == pd.DataFrame:
+        df.drop(objects, axis=1)
+    elif df_type == vaex.dataframe.DataFrameLocal:  # no objects in vaex
+        pass
+    else:
+        raise NotImplementedError(df_type)
+
     if verbose >= 2 and len(objects) > 0:
         print(
             (
@@ -238,7 +254,7 @@ def replace_PQR_with_m101(df):
     Parameters
     ----------
 
-    df: pandas Dataframe
+    df: pandas Dataframe or Vaex Dataframe
         ``branch`` must be a column name.
 
 
@@ -251,8 +267,19 @@ def replace_PQR_with_m101(df):
 
     # Note: somehow pandas updates dtype automatically. We have to check
     # We just have to replace the column:
+
+    # First checking if the input type is correct
+    df_type = type(df)
+    if not (df_type == pd.DataFrame or df_type == vaex.dataframe.DataFrameLocal):
+        raise NotImplementedError(df_type)
+
+    # Then do the actual mapping
     if df.dtypes["branch"] != np.int64:
-        new_col = df["branch"].replace(["P", "Q", "R"], [-1, 0, 1])
+        mapping = {"P": -1, "Q": 0, "R": 1}
+        if df_type == pd.DataFrame:  # pandas
+            new_col = df["branch"].replace(mapping)
+        else:  # vaex
+            new_col = df["branch"].map(mapping, allow_missing=True)
         df["branch"] = new_col
 
     try:
diff --git a/radis/default_radis.json b/radis/default_radis.json
index a3bf6f5b1..9f1eecb4f 100644
--- a/radis/default_radis.json
+++ b/radis/default_radis.json
@@ -24,7 +24,7 @@
     "SPARSE_WAVERANGE": "auto"              # true,/false. sparse LDM algorithm. May be smaller on dense spectra. If "auto", a scarcity criterion is used (Nlines/Ngrids > 1)
     "DEFAULT_DOWNLOAD_PATH": "~/.radisdb"   # default path for downloading databases with databank='hitran'/'hitemp'/'exomol' . You can also specify a local path for each entry of the, "database" list.
     "RESAMPLING_TOLERANCE_THRESHOLD": 5e-3  # an error if raises if areas do not match by a value above, this threshold during resampling. See :py:meth:`~radis.spectrum.spectrum.Spectrum.resample`
-
+    "DATAFRAME_ENGINE" : "pandas"
 
     # molecular parameters
     # --------------------
diff --git a/radis/io/exomol.py b/radis/io/exomol.py
index 982ee30bb..9a1d5626c 100644
--- a/radis/io/exomol.py
+++ b/radis/io/exomol.py
@@ -117,7 +117,7 @@ def fetch_exomol(
 
     Returns
     -------
-    df: pd.DataFrame
+    df: pd.DataFrame or vaex.dataframe.DataFrameLocal
         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
@@ -295,7 +295,7 @@ def fetch_exomol(
         mdb.rename_columns(df, {"Sij0": "int"})
 
     # Add Attributes of the DataFrame
-    if output == "pandas":  # no attributes in "Jax" or "Vaex" mode
+    if output in ["pandas", "vaex"]:  # no attribtes in "Jax" or "Vaex" mode
         from radis.db.classes import HITRAN_MOLECULES
 
         attrs = {}
@@ -306,9 +306,12 @@ def fetch_exomol(
         attrs["molecule"] = molecule
         attrs["iso"] = isotope
 
-        for k, v in attrs.items():
-            df.attrs[k] = v
-
+        if output == "vaex":
+            df.attrs = {}
+            df.attrs = attrs
+        elif output == "pandas":
+            for k, v in attrs.items():
+                df.attrs[k] = v
     # Return:
     out = df
     if return_local_path or return_partition_function:
diff --git a/radis/io/geisa.py b/radis/io/geisa.py
index 26ada5939..4b2abb175 100644
--- a/radis/io/geisa.py
+++ b/radis/io/geisa.py
@@ -80,7 +80,7 @@ def fetch_geisa(
 
     Returns
     -------
-    df: pd.DataFrame
+    df: pd.DataFrame or vaex.dataframe.DataFrameLocal
         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
diff --git a/radis/io/hitran.py b/radis/io/hitran.py
index b8d01b456..02437b95b 100644
--- a/radis/io/hitran.py
+++ b/radis/io/hitran.py
@@ -90,7 +90,7 @@ def fetch_hitran(
 
     Returns
     -------
-    df: pd.DataFrame
+    df: pd.DataFrame or vaex.dataframe.DataFrameLocal
         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
diff --git a/radis/io/query.py b/radis/io/query.py
index f6560d804..80b15a069 100644
--- a/radis/io/query.py
+++ b/radis/io/query.py
@@ -45,6 +45,7 @@ def fetch_astroquery(
     cache=True,
     expected_metadata={},
     engine="pytables-fixed",
+    output="pandas",
 ):
     """Download a HITRAN line database to a Pandas DataFrame.
 
@@ -76,7 +77,8 @@ def fetch_astroquery(
         if ``cache=True``, check that the metadata in the cache file correspond
         to these attributes. Arguments ``molecule``, ``isotope``, ``wmin``, ``wmax``
         are already added by default.
-
+    output : str
+        specifies the type of returned data
     References
     ----------
     .. [1] `Astroquery <https://astroquery.readthedocs.io>`_
@@ -222,10 +224,19 @@ def fetch_astroquery(
             raise ValueError(
                 "Error while parsing HITRAN output : {}".format(response.text)
             ) from err
-        df = tbl.to_pandas()
-        df = df.rename(columns=rename_columns)
+        if output == "pandas":
+            df = tbl.to_pandas()
+            df = df.rename(columns=rename_columns)
+        elif output == "vaex":
+            import vaex
+
+            df = vaex.from_astropy_table(tbl)
+            df = df.rename(columns=rename_columns)
     else:
         df = pd.DataFrame(columns=list(rename_columns.values()))
+        import vaex
+
+        df = vaex.from_pandas(df)
 
     # Cast type to float64
     cast_type = {
diff --git a/radis/lbl/bands.py b/radis/lbl/bands.py
index 9e5b915b2..c4c6ef79c 100644
--- a/radis/lbl/bands.py
+++ b/radis/lbl/bands.py
@@ -1100,7 +1100,7 @@ def dec(obj):
     return dec
 
 
-def add_bands(df, dbformat, lvlformat, verbose=True):
+def add_bands(df, dbformat, lvlformat, dataframe_type="pandas", verbose=True):
     """Assign all transitions to a vibrational band:
 
     Add 'band', 'viblvl_l' and 'viblvl_u' attributes for each line to allow
@@ -1277,9 +1277,14 @@ def add_bands(df, dbformat, lvlformat, verbose=True):
 
             vib_lvl_name = vib_lvl_name_hitran_class1
 
-            df.loc[:, "viblvl_l"] = vib_lvl_name(df["vl"])
-            df.loc[:, "viblvl_u"] = vib_lvl_name(df["vu"])
-            df.loc[:, "band"] = df["viblvl_l"] + "->" + df["viblvl_u"]
+            if dataframe_type == "pandas":
+                df.loc[:, "viblvl_l"] = vib_lvl_name(df["vl"])
+                df.loc[:, "viblvl_u"] = vib_lvl_name(df["vu"])
+                df.loc[:, "band"] = df["viblvl_l"] + "->" + df["viblvl_u"]
+            elif dataframe_type == "vaex":
+                df["viblvl_l"] = df.vl.apply(vib_lvl_name)
+                df["viblvl_u"] = df.vu.apply(vib_lvl_name)
+                df["band"] = df["viblvl_l"] + "->" + df["viblvl_u"]
 
         else:
             raise NotImplementedError(
diff --git a/radis/lbl/base.py b/radis/lbl/base.py
index bdde68594..7248efcec 100644
--- a/radis/lbl/base.py
+++ b/radis/lbl/base.py
@@ -66,6 +66,7 @@
 
 import numpy as np
 import pandas as pd
+import vaex
 from astropy import units as u
 from numpy import exp, pi
 from psutil import virtual_memory
@@ -80,7 +81,7 @@
 except ImportError:  # if ran from here
     from radis.lbl.loader import KNOWN_LVLFORMAT, DatabankLoader, df_metadata
 
-from radis.misc.arrays import anynan
+from radis.misc.arrays import anynan, anynan_vaex
 from radis.misc.basics import all_in, is_float, transfer_metadata
 from radis.misc.debug import printdbg
 from radis.misc.log import printwarn
@@ -89,7 +90,7 @@
 from radis.misc.utils import Default
 from radis.misc.warning import OutOfBoundError
 from radis.phys.constants import c_CGS, h_CGS, hc_k
-from radis.phys.convert import cm2J, nm2cm, nm_air2cm
+from radis.phys.convert import cm2J, cm2J_vaex, nm2cm, nm_air2cm
 from radis.phys.units_astropy import convert_and_strip_units
 from radis.spectrum.utils import print_conditions
 
@@ -311,22 +312,35 @@ def assert_no_nan(self, df, column):
         from radis.misc.printer import get_print_full
 
         try:
-            assert not anynan(df[column])
+            if self.dataframe_type == "pandas":
+                assert not anynan(df[column])
+            elif self.dataframe_type == "vaex":
+                assert not anynan_vaex(df[column])
         except AssertionError as err:
-            index = np.isnan(df[column]).idxmax()
-            if self.input.molecule == "CO2":
-                fix_idea = (
-                    "If using HITEMP2010 for CO2, some lines are unlabelled and therefore cannot be used at "
-                    "equilibrium. This is a known issue of the HITEMP database and will soon be fixed in the "
-                    "edition. In the meantime you can use:\n 'sf.df0.drop(sf.df0.index[sf.df0['v1u']==-1], inplace=True)' "
-                    "where 'sf' is SpectrumFactory object"
-                )
-            raise AssertionError(
-                "{0}=NaN in line database at index {1}".format(column, index)
-                + " corresponding to Line:\n {0}".format(
-                    get_print_full(df.loc[index]) + fix_idea
-                )
-            ) from err
+            if self.dataframe_type == "pandas":
+                index = np.isnan(df[column]).idxmax()
+                if self.input.molecule == "CO2":
+                    fix_idea = (
+                        "If using HITEMP2010 for CO2, some lines are unlabelled and therefore cannot be used at "
+                        "equilibrium. This is a known issue of the HITEMP database and will soon be fixed in the "
+                        "edition. In the meantime you can use:\n 'sf.df0.drop(sf.df0.index[sf.df0['v1u']==-1], inplace=True)' "
+                        "where 'sf' is SpectrumFactory object"
+                    )
+                raise AssertionError(
+                    "{0}=NaN in line database at index {1}".format(column, index)
+                    + " corresponding to Line:\n {0}".format(
+                        get_print_full(df.loc[index]) + fix_idea
+                    )
+                ) from err
+            elif self.dataframe_type == "vaex":
+                if self.input.molecule == "CO2":
+                    fix_idea = (
+                        "If using HITEMP2010 for CO2, some lines are unlabelled and therefore cannot be used at "
+                        "equilibrium. This is a known issue of the HITEMP database and will soon be fixed in the "
+                        "edition. In the meantime you can use:\n 'sf.df0.drop(sf.df0.index[sf.df0['v1u']==-1], inplace=True)' "
+                        "where 'sf' is SpectrumFactory object"
+                    )
+                raise AssertionError("Lines Have NaN Values")
 
     def _add_EvibErot(self, df, calc_Evib_harmonic_anharmonic=False):
         """Calculate Evib & Erot in Line dataframe.
@@ -539,48 +553,70 @@ def fillEvibl(r):
             #        df['Evibl'] = df.groupby('polyl','wangl','rankl').apply(Evibl, axis=1)
             #        %timeit: 43.4s per loop
 
-            # total:  ~ 15s on 460k lines   (probably faster since neq==0.9.20)
-            #            try:
-            # ~ 6.6 s (probably faster since neq==0.9.20 (radis<1.0)
-            df = df.groupby(by=["polyu", "wangu"]).apply(fillEvibu)
-            # ~ 6.6 s (probably faster since neq==0.9.20 (radis<1.0)
-            df = df.groupby(by=["polyl", "wangl"]).apply(fillEvibl)
-            # TODO : use map(dict) version (see _add_EvibErot_CDSD_pcN)
-
-            #            except KeyError:
-            #                import traceback
-            #                traceback.print_exc()
-            #                raise KeyError("{0} -> An error (see above) occurred that usually ".format(sys.exc_info()[1]) +
-            #                               "happens when the energy level is not referenced in the database. " +
-            #                               "Check your partition function calculator, and energies " +
-            #                               "for isotope {0} (Factory.parsum_calc['CO2'][{0}]['X'].df)".format(iso))
-
-            # Another version that failed because twice slower than apply() in that case
-            # ~ keep it for information
-            #        Evibdict = energies.set_index(['p','c','N'])['Evib']
-            #        Evibdict = Evibdict.drop_duplicates()
-            #        try:
-            #            dgb = df.groupby(by=['polyu', 'wangu', 'ranku'])
-            #            for (poly, wang, rank), idx in dgb.indices.items(): # ~ 3600 items for 460k lines -> total 15s
-            #                Evib = Evibdict[(poly, wang, rank)]             # ~ 7.15 µs
-            #                df.loc[idx, 'Evibu'] = Evib                     # ~ 4.38ms
-            #
-            #            dgb = df.groupby(by=['polyl', 'wangl', 'rankl'])
-            #            for (poly, wang, rank), idx in dgb.indices.items(): # ~ 3600 items for 460k lines -> total 15s
-            #                Evib = Evibdict[(poly, wang, rank)]             # ~ 7.15 µs
-            #                df.loc[idx, 'Evibl'] = Evib                     # ~ 4.38ms
-
-            return df.loc[idx, ["Evibl", "Evibu"]]
+            if self.dataframe_type == "pandas":
+                # total:  ~ 15s on 460k lines   (probably faster since neq==0.9.20)
+                #            try:
+                # ~ 6.6 s (probably faster since neq==0.9.20 (radis<1.0)
+                df = df.groupby(by=["polyu", "wangu"]).apply(fillEvibu)
+                # ~ 6.6 s (probably faster since neq==0.9.20 (radis<1.0)
+                df = df.groupby(by=["polyl", "wangl"]).apply(fillEvibl)
+
+                # TODO : use map(dict) version (see _add_EvibErot_CDSD_pcN)
+
+                #            except KeyError:
+                #                import traceback
+                #                traceback.print_exc()
+                #                raise KeyError("{0} -> An error (see above) occured that usually ".format(sys.exc_info()[1]) +
+                #                               "happens when the energy level is not referenced in the database. " +
+                #                               "Check your partition function calculator, and energies " +
+                #                               "for isotope {0} (Factory.parsum_calc['CO2'][{0}]['X'].df)".format(iso))
+
+                # Another version that failed because twice slower than apply() in that case
+                # ~ keep it for information
+                #        Evibdict = energies.set_index(['p','c','N'])['Evib']
+                #        Evibdict = Evibdict.drop_duplicates()
+                #        try:
+                #            dgb = df.groupby(by=['polyu', 'wangu', 'ranku'])
+                #            for (poly, wang, rank), idx in dgb.indices.items(): # ~ 3600 items for 460k lines -> total 15s
+                #                Evib = Evibdict[(poly, wang, rank)]             # ~ 7.15 µs
+                #                df.loc[idx, 'Evibu'] = Evib                     # ~ 4.38ms
+                #
+                #            dgb = df.groupby(by=['polyl', 'wangl', 'rankl'])
+                #            for (poly, wang, rank), idx in dgb.indices.items(): # ~ 3600 items for 460k lines -> total 15s
+                #                Evib = Evibdict[(poly, wang, rank)]             # ~ 7.15 µs
+                #                df.loc[idx, 'Evibl'] = Evib                     # ~ 4.38ms
+
+                return df.loc[idx, ["Evibl", "Evibu"]]
+
+            elif self.dataframe_type == "vaex":
+
+                def get_evib(poly, wang, evib, df_iso):
+                    if df_iso == iso:
+                        return energies.at((poly, wang), "Evib")
+                    else:
+                        return evib
+
+                df["Evibu"] = df.apply(get_evib, [df.polyu, df.wangu, df.Evibu, df.iso])
+                df["Evibl"] = df.apply(get_evib, [df.polyl, df.wangl, df.Evibl, df.iso])
 
         #        df = df.groupby('iso').apply(lambda x: add_Evib_CDSD_pc_1iso(x, x.name))
 
-        df["Evibl"] = np.nan
-        df["Evibu"] = np.nan
-        for iso, idx in df.groupby("iso").indices.items():
-            df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_CDSD_pc_1iso(df.loc[idx], iso)
+        if self.dataframe_type == "pandas":
+            df["Evibl"] = np.nan
+            df["Evibu"] = np.nan
+            for iso, idx in df.groupby("iso").indices.items():
+                df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_CDSD_pc_1iso(
+                    df.loc[idx], iso
+                )
+
+                if radis.config["DEBUG_MODE"]:
+                    assert (df.loc[idx, "iso"] == iso).all()
+        elif self.dataframe_type == "vaex":
+            df["Evibl"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evibu"] = vaex.vconstant(np.nan, df.length_unfiltered())
 
-            if radis.config["DEBUG_MODE"]:
-                assert (df.loc[idx, "iso"] == iso).all()
+            for iso in list(df.iso.unique()):
+                get_Evib_CDSD_pc_1iso(df, iso)
 
         # Get rotational energy: better recalculate than look up the database
         # (much faster!: perf ~25s -> 765µs)
@@ -666,31 +702,53 @@ def get_Evib_CDSD_pcN_1iso(df, iso):
             energies.set_index(index, inplace=True)
             Evib_dict = dict(list(zip(energies.index, energies.Evib)))
 
-            # Add lower state energy
-            df_pcN = df.set_index(["polyl", "wangl", "rankl"])
-            df["Evibl"] = df_pcN.index.map(Evib_dict.get).values
-            # Add upper state energy
-            df_pcN = df.set_index(["polyu", "wangu", "ranku"])
-            df["Evibu"] = df_pcN.index.map(Evib_dict.get).values
+            if self.dataframe_type == "pandas":
+                # Add lower state energy
+                df_pcN = df.set_index(["polyl", "wangl", "rankl"])
+                df["Evibl"] = df_pcN.index.map(Evib_dict.get).values
+                # Add upper state energy
+                df_pcN = df.set_index(["polyu", "wangu", "ranku"])
+                df["Evibu"] = df_pcN.index.map(Evib_dict.get).values
 
-            return df.loc[:, ["Evibl", "Evibu"]]
+                return df.loc[:, ["Evibl", "Evibu"]]
+            elif self.dataframe_type == "vaex":
 
-        df["Evibl"] = np.nan
-        df["Evibu"] = np.nan
+                def get_evib(poly, wang, rank, evib, iso_df):
+                    if iso_df == iso:
+                        return Evib_dict.get(poly, wang, rank)
+                    else:
+                        return evib
 
-        # multiple-isotopes in database
-        if "iso" in df:
-            for iso, idx in df.groupby("iso").indices.items():
-                df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_CDSD_pcN_1iso(
-                    df.loc[idx], iso
+                df["Evibl"] = df.apply(
+                    get_evib, [df.polyl, df.wangl, df.rankl, df.Evibl, df.iso]
+                )
+                df["Evibu"] = df.apply(
+                    get_evib, [df.polyu, df.wangu, df.ranku, df.Evibu, df.iso]
                 )
 
-                if radis.config["DEBUG_MODE"]:
-                    assert (df.loc[idx, "iso"] == iso).all()
+        if self.dataframe_type == "pandas":
+            df["Evibl"] = np.nan
+            df["Evibu"] = np.nan
 
-        else:
-            iso = df.attrs["iso"]
-            df.loc[:, ["Evibl", "Evibu"]] = get_Evib_CDSD_pcN_1iso(df, iso)
+            # multiple-isotopes in database
+            if "iso" in df:
+                for iso, idx in df.groupby("iso").indices.items():
+                    df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_CDSD_pcN_1iso(
+                        df.loc[idx], iso
+                    )
+
+                    if radis.config["DEBUG_MODE"]:
+                        assert (df.loc[idx, "iso"] == iso).all()
+
+            else:
+                iso = df.attrs["iso"]
+                df.loc[:, ["Evibl", "Evibu"]] = get_Evib_CDSD_pcN_1iso(df, iso)
+        elif self.dataframe_type == "vaex":
+            df["Evibl"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evibu"] = vaex.vconstant(np.nan, df.length_unfiltered())
+
+            for iso in list(df.iso.unique()):
+                get_Evib_CDSD_pcN_1iso(df, iso)
 
         # Get rotational energy: better recalculate than look up the database
         # (much faster!: perf ~25s -> 765µs)
@@ -778,40 +836,62 @@ def get_Evib_CDSD_pcJN_1iso(df, iso):
             energies = energies.set_index(index, inplace=False)
             Evib_dict = dict(list(zip(energies.index, energies.Evib)))
 
-            # Add lower state energy
-            df_pcJN = df.set_index(["polyl", "wangl", "jl", "rankl"])
-            #            for i in df.index:
-            #                df.loc[i, 'Evibl'] = energies.at[i, 'Evib']
-            # the map below is crazy fast compared to above loop
-            df.loc[:, "Evibl"] = df_pcJN.index.map(Evib_dict.get).values
-            # Add upper state energy
-            df_pcJN = df.set_index(["polyu", "wangu", "ju", "ranku"])
-            #            for i in df.index:
-            #                df.loc[i, 'Evibu'] = energies.at[i, 'Evib']
-            # the map below is crazy fast compared to above loop
-            df.loc[:, "Evibu"] = df_pcJN.index.map(Evib_dict.get).values
-
-            return df.loc[:, ["Evibl", "Evibu"]]
+            if self.dataframe_type == "pandas":
+                # Add lower state energy
+                df_pcJN = df.set_index(["polyl", "wangl", "jl", "rankl"])
+                #            for i in df.index:
+                #                df.loc[i, 'Evibl'] = energies.at[i, 'Evib']
+                # the map below is crazy fast compared to above loop
+                df.loc[:, "Evibl"] = df_pcJN.index.map(Evib_dict.get).values
+                # Add upper state energy
+                df_pcJN = df.set_index(["polyu", "wangu", "ju", "ranku"])
+                #            for i in df.index:
+                #                df.loc[i, 'Evibu'] = energies.at[i, 'Evib']
+                # the map below is crazy fast compared to above loop
+                df.loc[:, "Evibu"] = df_pcJN.index.map(Evib_dict.get).values
+
+                return df.loc[:, ["Evibl", "Evibu"]]
+            elif self.dataframe_type == "vaex":
+
+                def get_evib(poly, wang, j, rank, evib, iso_df):
+                    if iso_df == iso:
+                        return Evib_dict.get(poly, wang, j, rank)
+                    else:
+                        return evib
+
+                df["Evibl"] = df.apply(
+                    get_evib, [df.polyl, df.wangl, df.jl, df.rankl, df.Evibl, df.iso]
+                )
+                df["Evibu"] = df.apply(
+                    get_evib, [df.polyu, df.wangu, df.ju, df.ranku, df.Evibu, df.iso]
+                )
 
         #        df = df.groupby('iso').apply(lambda x: get_Evib_CDSD_pcJN_1iso(x, x.name))
 
-        # slower than the following:
-        df["Evibl"] = np.nan
-        df["Evibu"] = np.nan
+        if self.dataframe_type == "pandas":
+            # slower than the following:
+            df["Evibl"] = np.nan
+            df["Evibu"] = np.nan
 
-        # multiple-isotopes in database
-        if "iso" in df:
-            for iso, idx in df.groupby("iso").indices.items():
-                df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_CDSD_pcJN_1iso(
-                    df.loc[idx], iso
-                )
+            # multiple-isotopes in database
+            if "iso" in df:
+                for iso, idx in df.groupby("iso").indices.items():
+                    df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_CDSD_pcJN_1iso(
+                        df.loc[idx], iso
+                    )
 
-                if radis.config["DEBUG_MODE"]:
-                    assert (df.loc[idx, "iso"] == iso).all()
+                    if radis.config["DEBUG_MODE"]:
+                        assert (df.loc[idx, "iso"] == iso).all()
 
-        else:
-            iso = df.attrs["iso"]
-            df.loc[:, ["Evibl", "Evibu"]] = get_Evib_CDSD_pcJN_1iso(df, iso)
+            else:
+                iso = df.attrs["iso"]
+                df.loc[:, ["Evibl", "Evibu"]] = get_Evib_CDSD_pcJN_1iso(df, iso)
+        elif self.dataframe_type == "vaex":
+            df["Evibl"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evibu"] = vaex.vconstant(np.nan, df.length_unfiltered())
+
+            for iso in list(df.iso.unique()):
+                get_Evib_CDSD_pcJN_1iso(df, iso)
 
         # Get rotational energy: better recalculate than look up the database
         # (much faster!: perf ~25s -> 765µs)
@@ -926,35 +1006,85 @@ def fillEvib123l(r):
             #        df['Evibl'] = df.groupby('polyl','wangl','rankl').apply(Evibl, axis=1)
             #        %timeit: 43.4s per loop
 
-            #            try:  # total:  ~ 15s on 460k lines
-            # ~ 6.6 s   (probably faster since neq==0.9.20) (radis<1.0)
-            df = df.groupby(by=["polyu", "wangu"]).apply(fillEvib123u)
-            # ~ 6.6 s   (probably faster since neq==0.9.20) (radis<1.0)
-            df = df.groupby(by=["polyl", "wangl"]).apply(fillEvib123l)
-            #            except KeyError:
-            #                printr("{0} -> An error (see above) occurred that usually ".format(sys.exc_info()[1]) +
-            #                       "happens when the energy level is not referenced in the database. " +
-            #                       "Check your partition function calculator, and energies " +
-            #                       "for isotope {0} (Factory.parsum_calc['CO2'][{0}]['X'].df)".format(iso))
-            #                raise
-
-            return df.loc[
-                :, ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"]
-            ]
+            if self.dataframe_type == "pandas":
+                #            try:  # total:  ~ 15s on 460k lines
+                # ~ 6.6 s   (probably faster since neq==0.9.20) (radis<1.0)
+                df = df.groupby(by=["polyu", "wangu"]).apply(fillEvib123u)
+                # ~ 6.6 s   (probably faster since neq==0.9.20) (radis<1.0)
+                df = df.groupby(by=["polyl", "wangl"]).apply(fillEvib123l)
+                #            except KeyError:
+                #                printr("{0} -> An error (see above) occured that usually ".format(sys.exc_info()[1]) +
+                #                       "happens when the energy level is not referenced in the database. " +
+                #                       "Check your partition function calculator, and energies " +
+                #                       "for isotope {0} (Factory.parsum_calc['CO2'][{0}]['X'].df)".format(iso))
+                #                raise
+
+                return df.loc[
+                    :, ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"]
+                ]
+            elif self.dataframe_type == "vaex":
+
+                def get_evib1(poly, wang, evib, df_iso):
+                    if df_iso == iso:
+                        return energies.at((poly, wang), "Evib1")
+                    else:
+                        return evib
+
+                def get_evib2(poly, wang, evib, df_iso):
+                    if df_iso == iso:
+                        return energies.at((poly, wang), "Evib2")
+                    else:
+                        return evib
+
+                def get_evib3(poly, wang, evib, df_iso):
+                    if df_iso == iso:
+                        return energies.at((poly, wang), "Evib3")
+                    else:
+                        return evib
+
+                df["Evib1u"] = df.apply(
+                    get_evib1, [df.polyl, df.wangl, df.Evib1u, df.iso]
+                )
+                df["Evib2u"] = df.apply(
+                    get_evib2, [df.polyl, df.wangl, df.Evib2u, df.iso]
+                )
+                df["Evib3u"] = df.apply(
+                    get_evib3, [df.polyl, df.wangl, df.Evib3u, df.iso]
+                )
+                df["Evib1l"] = df.apply(
+                    get_evib1, [df.polyl, df.wangl, df.Evib1l, df.iso]
+                )
+                df["Evib2l"] = df.apply(
+                    get_evib2, [df.polyl, df.wangl, df.Evib2l, df.iso]
+                )
+                df["Evib3l"] = df.apply(
+                    get_evib3, [df.polyl, df.wangl, df.Evib3l, df.iso]
+                )
 
         #        df = df.groupby('iso').apply(lambda x: get_Evib123_CDSD_pc_1iso(x, x.name))
 
-        # Slower than the version below:
-        df["Evib1l"] = np.nan
-        df["Evib2l"] = np.nan
-        df["Evib3l"] = np.nan
-        df["Evib1u"] = np.nan
-        df["Evib2u"] = np.nan
-        df["Evib3u"] = np.nan
-        for iso, idx in df.groupby("iso").indices.items():
-            df.loc[
-                idx, ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"]
-            ] = get_Evib123_CDSD_pc_1iso(df.loc[idx], iso)
+        if self.dataframe_type == "pandas":
+            # Slower than the version below:
+            df["Evib1l"] = np.nan
+            df["Evib2l"] = np.nan
+            df["Evib3l"] = np.nan
+            df["Evib1u"] = np.nan
+            df["Evib2u"] = np.nan
+            df["Evib3u"] = np.nan
+            for iso, idx in df.groupby("iso").indices.items():
+                df.loc[
+                    idx, ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"]
+                ] = get_Evib123_CDSD_pc_1iso(df.loc[idx], iso)
+        elif self.dataframe_type == "vaex":
+            df["Evib1l"] = vaex.vconstant(np.nan, df.length_unfiltered)
+            df["Evib2l"] = vaex.vconstant(np.nan, df.length_unfiltered)
+            df["Evib3l"] = vaex.vconstant(np.nan, df.length_unfiltered)
+            df["Evib1u"] = vaex.vconstant(np.nan, df.length_unfiltered)
+            df["Evib2u"] = vaex.vconstant(np.nan, df.length_unfiltered)
+            df["Evib3u"] = vaex.vconstant(np.nan, df.length_unfiltered)
+
+            for iso in list(df.iso.unique()):
+                get_Evib123_CDSD_pc_1iso(df, iso)
 
         # Add total vibrational energy too (doesnt cost much, and can plot populations in spectrum)
         df["Evibu"] = df.Evib1u + df.Evib2u + df.Evib3u
@@ -1037,36 +1167,85 @@ def get_Evib_RADIS_cls1_1iso(df, iso):
             energies.set_index(index, inplace=True)
             Evib_dict = dict(list(zip(energies.index, energies.Evib)))
 
-            # Add lower state energy
-            df_v = df.set_index(["vl"])
-            df["Evibl"] = df_v.index.map(Evib_dict.get).values
+            if self.dataframe_type == "pandas":
+                # Add lower state energy
+                df_v = df.set_index(["vl"])
+                df["Evibl"] = df_v.index.map(Evib_dict.get).values
 
-            # Add upper state energy
-            df_v = df.set_index(["vu"])
-            df["Evibu"] = df_v.index.map(Evib_dict.get).values
+                # Add upper state energy
+                df_v = df.set_index(["vu"])
+                df["Evibu"] = df_v.index.map(Evib_dict.get).values
 
-            return df.loc[:, ["Evibl", "Evibu"]]
+                return df.loc[:, ["Evibl", "Evibu"]]
+            elif self.dataframe_type == "vaex":
 
-        df["Evibl"] = np.nan
-        df["Evibu"] = np.nan
+                def get_evib(iso_df, v, Evib):
+                    if iso_df == iso:
+                        return Evib_dict.get(v)
+                    else:
+                        return Evib
 
-        # multiple-isotopes in database
-        if "iso" in df:
-            for iso, idx in df.groupby("iso").indices.items():
-                df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_RADIS_cls1_1iso(
-                    df.loc[idx], iso
-                )
+                if "iso" in df:
+                    df["Evibl"] = df.apply(
+                        get_evib, arguments=[df.iso, df.vl, df.Evibl]
+                    )
+                    df["Evibu"] = df.apply(
+                        get_evib, arguments=[df.iso, df.vu, df.Evibu]
+                    )
+                else:
+                    df["Evibl"] = df["vl"].apply(lambda x: Evib_dict.get(x))
+                    df["Evibu"] = df["vu"].apply(lambda x: Evib_dict.get(x))
+
+                return df["Evibl"], df["Evibu"]
+
+        if self.dataframe_type == "pandas":
+            df["Evibl"] = np.nan
+            df["Evibu"] = np.nan
+            # multiple-isotopes in database
+            if "iso" in df:
+                for iso, idx in df.groupby("iso").indices.items():
+                    df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_RADIS_cls1_1iso(
+                        df.loc[idx], iso
+                    )
+            else:
+                iso = df.attrs["iso"]
+                df.loc[:, ["Evibl", "Evibu"]] = get_Evib_RADIS_cls1_1iso(df, iso)
 
-        else:
-            iso = df.attrs["iso"]
-            df.loc[:, ["Evibl", "Evibu"]] = get_Evib_RADIS_cls1_1iso(df, iso)
+            # Get rotational energy: better recalculate than look up the database (much faster!)
+            df["Erotu"] = df.Eu - df.Evibu
+            df["Erotl"] = df.El - df.Evibl
 
-        # Get rotational energy: better recalculate than look up the database (much faster!)
-        df["Erotu"] = df.Eu - df.Evibu
-        df["Erotl"] = df.El - df.Evibl
+            assert np.isnan(df.Evibu).sum() == 0
+            assert np.isnan(df.Evibl).sum() == 0
+
+        elif self.dataframe_type == "vaex":
+            df["Evibl"] = vaex.vconstant(np.nan, length=df.length_unfiltered())
+            df["Evibu"] = vaex.vconstant(np.nan, length=df.length_unfiltered())
+
+            # multiple-isotopes in database
+            if "iso" in df:
+                for iso in list(df.iso.unique()):
+                    # df["Evibl"], df["Evibu"] = get_Evib_RADIS_cls1_1iso(df, iso)
+                    get_Evib_RADIS_cls1_1iso(df, iso)
+            else:
+                iso = int(self.input.isotope)
+                get_Evib_RADIS_cls1_1iso(df, iso)
+            # if "iso" in df:
+            #     for iso, idx in df.groupby("iso").indices.items():
+            #         df.loc[idx, ["Evibl", "Evibu"]] = get_Evib_RADIS_cls1_1iso(
+            #             df.loc[idx], iso
+            #         )
+
+            # else:
+            #     iso = df.attrs["iso"]
+            #     df.loc[:, ["Evibl", "Evibu"]] = get_Evib_RADIS_cls1_1iso(df, iso)
+
+            # Get rotational energy: better recalculate than look up the database (much faster!)
+            df["Erotu"] = df.Eu - df.Evibu
+            df["Erotl"] = df.El - df.Evibl
 
-        assert np.isnan(df.Evibu).sum() == 0
-        assert np.isnan(df.Evibl).sum() == 0
+            assert df.Evibu.isna().sum() == 0
+            assert df.Evibl.isna().sum() == 0
 
         self.profiler.stop(
             "fetch_energy_5", "Fetched energies for all {0} transitions".format(len(df))
@@ -1139,50 +1318,103 @@ def get_Evib123_RADIS_cls5_1iso(df, iso):
             Evib2_dict = dict(list(zip(energies.index, energies.Evib2)))
             Evib3_dict = dict(list(zip(energies.index, energies.Evib3)))
 
-            # Add lower state energy
-            df_v1v2l2v3 = df.set_index(["v1l", "v2l", "l2l", "v3l"])
-            #            for i in df.index:
-            #                df.loc[i, 'Evib1l'] = energies.at[i, 'Evib1']
-            # the map below is crazy fast compared to above loop
-            df["Evib1l"] = df_v1v2l2v3.index.map(Evib1_dict.get).values
-            df["Evib2l"] = df_v1v2l2v3.index.map(Evib2_dict.get).values
-            df["Evib3l"] = df_v1v2l2v3.index.map(Evib3_dict.get).values
-            # TODO @dev # performance: try getting all 3 values at the same time?
-            # with:  Evib123_dict = dict(list(zip(energies.index, (energies.Evib1, energies.Evib2, energies.Evib3))))
-
-            # Add upper state energy
-            df_v1v2l2v3 = df.set_index(["v1u", "v2u", "l2u", "v3u"])
-            df["Evib1u"] = df_v1v2l2v3.index.map(Evib1_dict.get).values
-            df["Evib2u"] = df_v1v2l2v3.index.map(Evib2_dict.get).values
-            df["Evib3u"] = df_v1v2l2v3.index.map(Evib3_dict.get).values
+            if self.dataframe_type == "pandas":
+                # Add lower state energy
+                df_v1v2l2v3 = df.set_index(["v1l", "v2l", "l2l", "v3l"])
+                #            for i in df.index:
+                #                df.loc[i, 'Evib1l'] = energies.at[i, 'Evib1']
+                # the map below is crazy fast compared to above loop
+                df["Evib1l"] = df_v1v2l2v3.index.map(Evib1_dict.get).values
+                df["Evib2l"] = df_v1v2l2v3.index.map(Evib2_dict.get).values
+                df["Evib3l"] = df_v1v2l2v3.index.map(Evib3_dict.get).values
+                # TODO @dev # performance: try getting all 3 values at the same time?
+                # with:  Evib123_dict = dict(list(zip(energies.index, (energies.Evib1, energies.Evib2, energies.Evib3))))
 
-            return df.loc[
-                :, ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"]
-            ]
+                # Add upper state energy
+                df_v1v2l2v3 = df.set_index(["v1u", "v2u", "l2u", "v3u"])
+                df["Evib1u"] = df_v1v2l2v3.index.map(Evib1_dict.get).values
+                df["Evib2u"] = df_v1v2l2v3.index.map(Evib2_dict.get).values
+                df["Evib3u"] = df_v1v2l2v3.index.map(Evib3_dict.get).values
+
+                return df.loc[
+                    :, ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"]
+                ]
+            elif self.dataframe_type == "vaex":
+
+                def get_evib1(v1, v2, l2, v3, evib1, iso_df):
+                    if iso_df == iso:
+                        return Evib1_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib1
+
+                def get_evib2(v1, v2, l2, v3, evib2, iso_df):
+                    if iso_df == iso:
+                        return Evib2_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib2
+
+                def get_evib3(v1, v2, l2, v3, evib3, iso_df):
+                    if iso_df == iso:
+                        return Evib3_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib3
+
+                # Add lower state energy
+                df["Evib1l"] = df.apply(
+                    get_evib1, [df.v1, df.v2, df.l2, df.v3, df.Evib1l, df.iso]
+                )
+                df["Evib2l"] = df.apply(
+                    get_evib2, [df.v1, df.v2, df.l2, df.v3, df.Evib2l, df.iso]
+                )
+                df["Evib3l"] = df.apply(
+                    get_evib3, [df.v1, df.v2, df.l2, df.v3, df.Evib3l, df.iso]
+                )
+
+                # Add upper state energy
+                df["Evib1u"] = df.apply(
+                    get_evib1, [df.v1, df.v2, df.l2, df.v3, df.Evib1l, df.iso]
+                )
+                df["Evib2u"] = df.apply(
+                    get_evib2, [df.v1, df.v2, df.l2, df.v3, df.Evib2l, df.iso]
+                )
+                df["Evib3u"] = df.apply(
+                    get_evib3, [df.v1, df.v2, df.l2, df.v3, df.Evib3l, df.iso]
+                )
 
         #        df = df.groupby('iso').apply(lambda x: get_Evib123_RADIS_cls5_1iso(x, x.name))
 
-        # Slower than the version below:
-        df["Evib1l"] = np.nan
-        df["Evib2l"] = np.nan
-        df["Evib3l"] = np.nan
-        df["Evib1u"] = np.nan
-        df["Evib2u"] = np.nan
-        df["Evib3u"] = np.nan
+        if self.dataframe_type == "pandas":
+            # Slower than the version below:
+            df["Evib1l"] = np.nan
+            df["Evib2l"] = np.nan
+            df["Evib3l"] = np.nan
+            df["Evib1u"] = np.nan
+            df["Evib2u"] = np.nan
+            df["Evib3u"] = np.nan
+            # multiple-isotopes in database
+            if "iso" in df:
+                for iso, idx in df.groupby("iso").indices.items():
+                    df.loc[
+                        idx,
+                        ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"],
+                    ] = get_Evib123_RADIS_cls5_1iso(df.loc[idx], iso)
 
-        # multiple-isotopes in database
-        if "iso" in df:
-            for iso, idx in df.groupby("iso").indices.items():
+            else:
+                iso = df.attrs["iso"]
                 df.loc[
-                    idx, ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"]
-                ] = get_Evib123_RADIS_cls5_1iso(df.loc[idx], iso)
-
-        else:
-            iso = df.attrs["iso"]
-            df.loc[
-                :,
-                ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"],
-            ] = get_Evib123_RADIS_cls5_1iso(df, iso)
+                    :,
+                    ["Evib1l", "Evib2l", "Evib3l", "Evib1u", "Evib2u", "Evib3u"],
+                ] = get_Evib123_RADIS_cls5_1iso(df, iso)
+        elif self.dataframe_type == "vaex":
+            df["Evib1l"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib2l"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib3l"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib1u"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib2u"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib3u"] = vaex.vconstant(np.nan, df.length_unfiltered())
+
+            for iso in list(df.iso.unique()):
+                get_Evib123_RADIS_cls5_1iso(df, iso)
 
         # Add total vibrational energy too (doesnt cost much, and can plot populations in spectrum)
         df["Evibu"] = df.Evib1u + df.Evib2u + df.Evib3u
@@ -1267,26 +1499,84 @@ def get_Evib123_RADIS_cls5_1iso_ah(df, iso):
             Evib3_h_dict = dict(list(zip(energies.index, energies.Evib3_h)))
             Evib3_a_dict = dict(list(zip(energies.index, energies.Evib3_a)))
 
-            # Add lower state energy
-            df_v1v2l2v3 = df.set_index(["v1l", "v2l", "l2l", "v3l"])
-            # the map below is crazy fast compared to above loop
-            df["Evib1l_h"] = df_v1v2l2v3.index.map(Evib1_h_dict.get).values
-            df["Evib1l_a"] = df_v1v2l2v3.index.map(Evib1_a_dict.get).values
-            df["Evib2l_h"] = df_v1v2l2v3.index.map(Evib2_h_dict.get).values
-            df["Evib2l_a"] = df_v1v2l2v3.index.map(Evib2_a_dict.get).values
-            df["Evib3l_h"] = df_v1v2l2v3.index.map(Evib3_h_dict.get).values
-            df["Evib3l_a"] = df_v1v2l2v3.index.map(Evib3_a_dict.get).values
-            # TODO @dev # performance: try getting all 3 values at the same time?
-            # with:  Evib123_dict = dict(list(zip(energies.index, (energies.Evib1, energies.Evib2, energies.Evib3))))
-
-            # Add upper state energy
-            df_v1v2l2v3 = df.set_index(["v1u", "v2u", "l2u", "v3u"])
-            df["Evib1u_h"] = df_v1v2l2v3.index.map(Evib1_h_dict.get).values
-            df["Evib1u_a"] = df_v1v2l2v3.index.map(Evib1_a_dict.get).values
-            df["Evib2u_h"] = df_v1v2l2v3.index.map(Evib2_h_dict.get).values
-            df["Evib2u_a"] = df_v1v2l2v3.index.map(Evib2_a_dict.get).values
-            df["Evib3u_h"] = df_v1v2l2v3.index.map(Evib3_h_dict.get).values
-            df["Evib3u_a"] = df_v1v2l2v3.index.map(Evib3_a_dict.get).values
+            if self.dataframe_type == "pandas":
+                # Add lower state energy
+                df_v1v2l2v3 = df.set_index(["v1l", "v2l", "l2l", "v3l"])
+                # the map below is crazy fast compared to above loop
+                df["Evib1l_h"] = df_v1v2l2v3.index.map(Evib1_h_dict.get).values
+                df["Evib1l_a"] = df_v1v2l2v3.index.map(Evib1_a_dict.get).values
+                df["Evib2l_h"] = df_v1v2l2v3.index.map(Evib2_h_dict.get).values
+                df["Evib2l_a"] = df_v1v2l2v3.index.map(Evib2_a_dict.get).values
+                df["Evib3l_h"] = df_v1v2l2v3.index.map(Evib3_h_dict.get).values
+                df["Evib3l_a"] = df_v1v2l2v3.index.map(Evib3_a_dict.get).values
+                # TODO @dev # performance: try getting all 3 values at the same time?
+                # with:  Evib123_dict = dict(list(zip(energies.index, (energies.Evib1, energies.Evib2, energies.Evib3))))
+
+                # Add upper state energy
+                df_v1v2l2v3 = df.set_index(["v1u", "v2u", "l2u", "v3u"])
+                df["Evib1u_h"] = df_v1v2l2v3.index.map(Evib1_h_dict.get).values
+                df["Evib1u_a"] = df_v1v2l2v3.index.map(Evib1_a_dict.get).values
+                df["Evib2u_h"] = df_v1v2l2v3.index.map(Evib2_h_dict.get).values
+                df["Evib2u_a"] = df_v1v2l2v3.index.map(Evib2_a_dict.get).values
+                df["Evib3u_h"] = df_v1v2l2v3.index.map(Evib3_h_dict.get).values
+                df["Evib3u_a"] = df_v1v2l2v3.index.map(Evib3_a_dict.get).values
+            elif self.dataframe_type == "vaex":
+
+                def get_Evib1_h(v1, v2, l2, v3, evib1_h, iso_df):
+                    if iso_df == iso:
+                        return Evib1_h_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib1_h
+
+                def get_Evib1_a(v1, v2, l2, v3, evib1_a, iso_df):
+                    if iso_df == iso:
+                        return Evib1_a_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib1_a
+
+                def get_Evib2_h(v1, v2, l2, v3, evib2_h, iso_df):
+                    if iso_df == iso:
+                        return Evib2_h_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib2_h
+
+                def get_Evib2_a(v1, v2, l2, v3, evib2_a, iso_df):
+                    if iso_df == iso:
+                        return Evib2_a_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib2_a
+
+                def get_Evib3_h(v1, v2, l2, v3, evib3_h, iso_df):
+                    if iso_df == iso:
+                        return Evib3_h_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib3_h
+
+                def get_Evib3_a(v1, v2, l2, v3, evib3_a, iso_df):
+                    if iso_df == iso:
+                        return Evib3_a_dict.get(v1, v2, l2, v3)
+                    else:
+                        return evib3_a
+
+                # Add lower state energy
+                df["Evib1l_h"] = df.apply(
+                    get_Evib1_h, [df.v1, df.v2, df.l2, df.v3, df.Evib1l_h]
+                )
+                df["Evib1l_a"] = df.apply(
+                    get_Evib1_a, [df.v1, df.v2, df.l2, df.v3, df.Evib1l_a]
+                )
+                df["Evib2l_h"] = df.apply(
+                    get_Evib2_h, [df.v1, df.v2, df.l2, df.v3, df.Evib2l_h]
+                )
+                df["Evib2l_a"] = df.apply(
+                    get_Evib2_a, [df.v1, df.v2, df.l2, df.v3, df.Evib2l_a]
+                )
+                df["Evib3l_h"] = df.apply(
+                    get_Evib3_h, [df.v1, df.v2, df.l2, df.v3, df.Evib3l_h]
+                )
+                df["Evib3l_a"] = df.apply(
+                    get_Evib3_a, [df.v1, df.v2, df.l2, df.v3, df.Evib3l_a]
+                )
 
             return df.loc[
                 :,
@@ -1306,37 +1596,54 @@ def get_Evib123_RADIS_cls5_1iso_ah(df, iso):
                 ],
             ]
 
-        # Slower than the version below:
-        df["Evib1l_h"] = np.nan
-        df["Evib1l_a"] = np.nan
-        df["Evib2l_h"] = np.nan
-        df["Evib2l_a"] = np.nan
-        df["Evib3l_h"] = np.nan
-        df["Evib3l_a"] = np.nan
-        df["Evib1u_h"] = np.nan
-        df["Evib1u_a"] = np.nan
-        df["Evib2u_h"] = np.nan
-        df["Evib2u_a"] = np.nan
-        df["Evib3u_h"] = np.nan
-        df["Evib3u_a"] = np.nan
-        for iso, idx in df.groupby("iso").indices.items():
-            df.loc[
-                idx,
-                [
-                    "Evib1l_h",
-                    "Evib1l_a",
-                    "Evib2l_h",
-                    "Evib2l_a",
-                    "Evib3l_h",
-                    "Evib3l_a",
-                    "Evib1u_h",
-                    "Evib1u_a",
-                    "Evib2u_h",
-                    "Evib2u_a",
-                    "Evib3u_h",
-                    "Evib3u_a",
-                ],
-            ] = get_Evib123_RADIS_cls5_1iso_ah(df.loc[idx], iso)
+        if self.dataframe_type == "pandas":
+            # Slower than the version below:
+            df["Evib1l_h"] = np.nan
+            df["Evib1l_a"] = np.nan
+            df["Evib2l_h"] = np.nan
+            df["Evib2l_a"] = np.nan
+            df["Evib3l_h"] = np.nan
+            df["Evib3l_a"] = np.nan
+            df["Evib1u_h"] = np.nan
+            df["Evib1u_a"] = np.nan
+            df["Evib2u_h"] = np.nan
+            df["Evib2u_a"] = np.nan
+            df["Evib3u_h"] = np.nan
+            df["Evib3u_a"] = np.nan
+            for iso, idx in df.groupby("iso").indices.items():
+                df.loc[
+                    idx,
+                    [
+                        "Evib1l_h",
+                        "Evib1l_a",
+                        "Evib2l_h",
+                        "Evib2l_a",
+                        "Evib3l_h",
+                        "Evib3l_a",
+                        "Evib1u_h",
+                        "Evib1u_a",
+                        "Evib2u_h",
+                        "Evib2u_a",
+                        "Evib3u_h",
+                        "Evib3u_a",
+                    ],
+                ] = get_Evib123_RADIS_cls5_1iso_ah(df.loc[idx], iso)
+        elif self.dataframe_type == "vaex":
+            df["Evib1l_h"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib1l_a"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib2l_h"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib2l_a"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib3l_h"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib3l_a"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib1u_h"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib1u_a"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib2u_h"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib2u_a"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib2u_h"] = vaex.vconstant(np.nan, df.length_unfiltered())
+            df["Evib2u_a"] = vaex.vconstant(np.nan, df.length_unfiltered())
+
+            for iso in list(df.iso.unique()):
+                get_Evib123_RADIS_cls5_1iso_ah(df, iso)
 
         # Add total vibrational energy too (doesnt cost much, and can plot populations in spectrum)
         df["Evibu_a"] = df.Evib1u_a + df.Evib2u_a + df.Evib3u_a
@@ -1408,13 +1715,21 @@ def _add_ju(self, df):
         #        df.loc[df.branch==1,'ju'] += 1     # branch R
 
         #        # slightly less readable but ~ 20% faster than above:
-        dgb = df.groupby("branch")
-        df["ju"] = df.jl
-        for branch, idx in dgb.indices.items():
-            if branch == -1:  # 'P':
-                df.loc[idx, "ju"] -= 1
-            if branch == 1:  #'R':
-                df.loc[idx, "ju"] += 1
+
+        if self.dataframe_type == "pandas":
+            dgb = df.groupby("branch")
+            df["ju"] = df.jl
+            for branch, idx in dgb.indices.items():
+                if branch == -1:  # 'P':
+                    df.loc[idx, "ju"] -= 1
+                if branch == 1:  #'R':
+                    df.loc[idx, "ju"] += 1
+        elif self.dataframe_type == "vaex":
+
+            def fun(branch, j):
+                return j + branch
+
+            df["ju"] = df.apply(f=fun, arguments=[df.branch, df.jl])
 
         return None
 
@@ -1512,11 +1827,20 @@ def _calc_noneq_parameters(self, vib_distribution, singleTvibmode):
 
         # ... Check no negative energies
         tol = -1e-4  # tolerance for negative energies (in cm-1)
-        if not ((df.Erotu > tol).all() and (df.Erotl > tol).all()):
-            self.warn(
-                "There are negative rotational energies in the database",
-                "NegativeEnergiesWarning",
-            )
+        if self.dataframe_type == "pandas":
+            if not ((df.Erotu > tol).all() and (df.Erotl > tol).all()):
+                self.warn(
+                    "There are negative rotational energies in the database",
+                    "NegativeEnergiesWarning",
+                )
+        elif self.dataframe_type == "vaex":
+            if not (
+                (df.Erotu > tol).sum() == len(df) and (df.Erotl > tol).sum() == len(df)
+            ):
+                self.warn(
+                    "There are negative rotational energies in the database",
+                    "NegativeEnergiesWarning",
+                )
 
         # ... Make sure degeneracies are calculated
         if not all_in(["gju", "gjl", "gvibu", "gvibl", "gu", "gl"], df):
@@ -1559,29 +1883,54 @@ def _calc_degeneracies(self, df):
         if "iso" in df:  # multiple isotopes in database
             id = df.attrs["id"]
             dgb = df.groupby(by=["iso"])
-            for (iso), idx in dgb.indices.items():
-                _gs = gs(id, iso)
-                if isinstance(_gs, tuple):
-                    # Molecules that have alternating degeneracy.
-                    if id not in [2]:  # CO2
-                        raise NotImplementedError
-                    # normally we should find whether the rovibrational level is symmetric
-                    # or asymmetric. Here we just assume it's symmetric, because for
-                    # symmetric isotopes such as CO2(626), CO2 asymmetric levels
-                    # dont exist (gs=0) and they should not be in the line database.
-                    _gs = _gs[0]
-
-                dg = df.loc[idx]
-                _gi = gi(id, iso)
-                df.loc[idx, "grotu"] = dg.gju * _gs * _gi
-                df.loc[idx, "grotl"] = dg.gjl * _gs * _gi
-
-                if radis.config["DEBUG_MODE"]:
-                    assert (df.loc[idx, "iso"] == iso).all()
+            if self.dataframe_type == "pandas":
+                for (iso), idx in dgb.indices.items():
+                    _gs = gs(id, iso)
+                    if isinstance(_gs, tuple):
+                        # Molecules that have alternating degeneracy.
+                        if id not in [2]:  # CO2
+                            raise NotImplementedError
+                        # normally we should find whether the rovibrational level is symmetric
+                        # or asymmetric. Here we just assume it's symmetric, because for
+                        # symmetric isotopes such as CO2(626), CO2 asymmetric levels
+                        # dont exist (gs=0) and they should not be in the line database.
+                        _gs = _gs[0]
+
+                    dg = df.loc[idx]
+                    _gi = gi(id, iso)
+                    df.loc[idx, "grotu"] = dg.gju * _gs * _gi
+                    df.loc[idx, "grotl"] = dg.gjl * _gs * _gi
+
+                    if radis.config["DEBUG_MODE"]:
+                        assert (df.loc[idx, "iso"] == iso).all()
+            elif self.dataframe_type == "vaex":
+
+                _gs_dic = {}
+                _gi_dic = {}
+
+                for iso in list(df.iso.unique()):
+                    _gs = gs(id, iso)
+                    if isinstance(_gs, tuple):
+                        # Molecules that have alternating degeneracy.
+                        if id not in [2]:  # CO2
+                            raise NotImplementedError
+                        # normally we should find whether the rovibrational level is symmetric
+                        # or asymmetric. Here we just assume it's symmetric, because for
+                        # symmetric isotopes such as CO2(626), CO2 asymmetric levels
+                        # dont exist (gs=0) and they should not be in the line database.
+                        _gs = _gs[0]
+                    _gs_dic[iso] = _gs
+                    _gi_dic[iso] = _gi = gi(id, iso)
+
+                def grot(gj, iso):
+                    return gj * _gs_dic[iso] * _gi_dic[iso]
+
+                df["grotu"] = df.apply(f=grot, arguments=[df.gju, df.iso])
+                df["grotl"] = df.apply(f=grot, arguments=[df.gjl, df.iso])
 
         else:
             id = df.attrs["id"]
-            isotope = df.attrs["iso"]
+            isotope = int(self.input.isotope)
             _gs = gs(id, isotope)
             if isinstance(_gs, tuple):
                 # Molecules that have alternating degeneracy.
@@ -1595,8 +1944,12 @@ def _calc_degeneracies(self, df):
 
             dg = df
             _gi = gi(id, isotope)
-            df.loc[:, "grotu"] = dg.gju * _gs * _gi
-            df.loc[:, "grotl"] = dg.gjl * _gs * _gi
+            if self.dataframe_type == "pandas":
+                df.loc[:, "grotu"] = dg.gju * _gs * _gi
+                df.loc[:, "grotl"] = dg.gjl * _gs * _gi
+            elif self.dataframe_type == "vaex":
+                df["grotu"] = dg.gju * _gs * _gi
+                df["grotl"] = dg.gjl * _gs * _gi
 
         # %%
 
@@ -1613,8 +1966,12 @@ def _calc_degeneracies(self, df):
             # Complete rovibrational assignment would not be True, for instance, for
             # bending levels of CO2 if all levels are considered degenerated
             # (with a v2+1 degeneracy)
-            df["gvibu"] = 1
-            df["gvibl"] = 1
+            if self.dataframe_type == "pandas":
+                df["gvibu"] = 1
+                df["gvibl"] = 1
+            elif self.dataframe_type == "vaex":
+                df["gvibu"] = vaex.vconstant(1, length=df.length_unfiltered())
+                df["gvibl"] = vaex.vconstant(1, length=df.length_unfiltered())
         else:
             raise NotImplementedError(
                 "vibrational degeneracy assignation for dbformat={0}".format(dbformat)
@@ -1627,15 +1984,26 @@ def _calc_degeneracies(self, df):
         # Check consistency if "gp" already existed
         # https://github.com/radis/radis/pull/514#issuecomment-1229463074
         if "gp" in df:
-            if (df["gp"] == df["gu"]).all():
-                self.warn(
-                    "'gu' was recomputed although 'gp' already in DataFrame. All values are equal",
-                    "PerformanceWarning",
-                )
-            else:
-                raise ValueError(
-                    "'gu' was recomputed although 'gp' already in DataFrame. Values are not equal ! Check calculations"
-                )
+            if self.dataframe_type == "pandas":
+                if (df["gp"] == df["gu"]).all():
+                    self.warn(
+                        "'gu' was recomputed although 'gp' already in DataFrame. All values are equal",
+                        "PerformanceWarning",
+                    )
+                else:
+                    raise ValueError(
+                        "'gu' was recomputed although 'gp' already in DataFrame. Values are not equal ! Check calculations"
+                    )
+            elif self.dataframe_type == "vaex":
+                print("hello")
+                # if (df["gp"] == df["gu"]).sum() == len(df):
+                #     self.warn(
+                #         "'gu' was recomputed although 'gp' already in DataFrame. All values are equal",
+                #     )
+                # else:
+                #     raise ValueError(
+                #         "'gu' was recomputed although 'gp' already in DataFrame. Values are not equal ! Check calculations"
+                #     )
 
         return None  # dataframe updated directly
 
@@ -1961,7 +2329,12 @@ def calc_lineshift(self):
             df["shiftwav"] = df.wav
 
         # Sorted lines is needed for sparse wavenumber range algorithm.
-        df.sort_values("shiftwav", inplace=True)
+        if self.dataframe_type == "pandas":
+            df.sort_values("shiftwav", kind="mergesort", inplace=True)
+        elif self.dataframe_type == "vaex":
+            attrs = df.attrs
+            self.df1 = df.sort("shiftwav")
+            self.df1.attrs = attrs
 
         self.profiler.stop("calc_lineshift", "Calculated lineshift")
 
@@ -2221,6 +2594,7 @@ def calc_linestrength_eq(self, Tgas):
 
         Tref = self.input.Tref
         df1 = self.df1
+        df1.attrs = self.df1.attrs
 
         if len(df1) == 0:
             return  # no lines
@@ -2442,54 +2816,116 @@ def Qneq_Qvib_Qrotu_Qrotl(
                     "PerformanceWarning",
                 )
 
-            dgb = df.groupby(by=["iso"])
-            for (iso), idx in dgb.indices.items():
-
-                # Get partition function for all lines
-                parsum = self.get_partition_function_calculator(molecule, iso, state)
+            if self.dataframe_type == "pandas":
+                dgb = df.groupby(by=["iso"])
+                for (iso), idx in dgb.indices.items():
 
-                if is_float(Tvib):
-                    Q, Qvib, dfQrot = parsum.at_noneq(
-                        Tvib,
-                        Trot,
-                        vib_distribution=vib_distribution,
-                        rot_distribution=rot_distribution,
-                        overpopulation=overpopulation,
-                        returnQvibQrot=True,
-                        update_populations=self.misc.export_populations,
-                    )
-                else:
-                    raise NotImplementedError(
-                        "Cannot return detail of Qvib, Qrot for 3-Tvib mode"
+                    # Get partition function for all lines
+                    parsum = self.get_partition_function_calculator(
+                        molecule, iso, state
                     )
 
-                # ... make sure PartitionFunction above is calculated with the same
-                # ... temperatures, rovibrational distributions and overpopulations
-                # ... as the populations of active levels (somewhere below)
-                df.loc[idx, "Qvib"] = Qvib
-                df.loc[idx, "Q"] = Q
-
-                # reindexing to get a direct access to Qrot database
-                # create the lookup dictionary
-                # dfQrot index is already 'viblvl'
-                dfQrot_dict = dict(list(zip(dfQrot.index, dfQrot.Qrot)))
+                    if is_float(Tvib):
+                        Q, Qvib, dfQrot = parsum.at_noneq(
+                            Tvib,
+                            Trot,
+                            vib_distribution=vib_distribution,
+                            rot_distribution=rot_distribution,
+                            overpopulation=overpopulation,
+                            returnQvibQrot=True,
+                            update_populations=self.misc.export_populations,
+                        )
+                    else:
+                        raise NotImplementedError(
+                            "Cannot return detail of Qvib, Qrot for 3-Tvib mode"
+                        )
+
+                    # ... make sure PartitionFunction above is calculated with the same
+                    # ... temperatures, rovibrational distributions and overpopulations
+                    # ... as the populations of active levels (somewhere below)
+                    print(df.columns)
+                    df.loc[idx, "Qvib"] = Qvib
+                    df.loc[idx, "Q"] = Q
+
+                    # reindexing to get a direct access to Qrot database
+                    # create the lookup dictionary
+                    # dfQrot index is already 'viblvl'
+                    dfQrot_dict = dict(list(zip(dfQrot.index, dfQrot.Qrot)))
+
+                    dg = df.loc[idx]
+
+                    # Add lower state Qrot
+                    dg_sorted = dg.set_index(["viblvl_l"], inplace=False)
+                    df.loc[idx, "Qrotl"] = dg_sorted.index.map(dfQrot_dict.get).values
+                    # Add upper state energy
+                    dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
+                    df.loc[idx, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
+
+                    if radis.config["DEBUG_MODE"]:
+                        assert (df.loc[idx, "iso"] == iso).all()
+            elif self.dataframe_type == "vaex":
+                df["Qvib"] = vaex.vconstant(np.nan, df.length_unfiltered())
+                df["Q"] = vaex.vconstant(np.nan, df.length_unfiltered())
+
+                Q_dict = {}
+                Qvib_dict = {}
+                dfQrot_dict = {}
+                for iso in iso_set:
+                    # Get partition function for all lines
+                    parsum = self.get_partition_function_calculator(
+                        molecule, iso, state
+                    )
 
-                dg = df.loc[idx]
+                    if is_float(Tvib):
+                        Q_dict[iso], Qvib_dict[iso], dfQrot = parsum.at_noneq(
+                            Tvib,
+                            Trot,
+                            vib_distribution=vib_distribution,
+                            rot_distribution=rot_distribution,
+                            overpopulation=overpopulation,
+                            returnQvibQrot=True,
+                            update_populations=self.misc.export_populations,
+                        )
+                    else:
+                        raise NotImplementedError(
+                            "Cannot return detail of Qvib, Qrot for 3-Tvib mode"
+                        )
+                    # reindexing to get a direct access to Qrot database
+                    # create the lookup dictionary
+                    # dfQrot index is already 'viblvl'
+                    dfQrot_dict[iso] = dict(list(zip(dfQrot.index, dfQrot.Qrot)))
+
+                def update_Qvib(iso):
+                    return Qvib_dict[iso]
+
+                def update_Q(iso):
+                    return Q_dict[iso]
+
+                df["Qvib"] = df.apply(update_Qvib, arguments=[df.iso])
+                df["Q"] = df.apply(update_Q, arguments=[df.iso])
+
+                # # Add lower state Qrot
+                # dg_sorted = dg.set_index(["viblvl_l"], inplace=False)
+                # df.loc[idx, "Qrotl"] = dg_sorted.index.map(dfQrot_dict.get).values
+                # # Add upper state energy
+                # dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
+                # df.loc[idx, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
+                def update_Qrot(viblvl, iso):
+                    return dfQrot_dict[iso].get(viblvl)
+                    # if iso_df == iso:
+                    #     return dfQrot_dict.get(val)
+                    # else:
+                    #     return val
 
                 # Add lower state Qrot
-                dg_sorted = dg.set_index(["viblvl_l"], inplace=False)
-                df.loc[idx, "Qrotl"] = dg_sorted.index.map(dfQrot_dict.get).values
+                df["Qrotl"] = df.apply(update_Qrot, arguments=[df.viblvl_l, df.iso])
                 # Add upper state energy
-                dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
-                df.loc[idx, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
-
-                if radis.config["DEBUG_MODE"]:
-                    assert (df.loc[idx, "iso"] == iso).all()
+                df["Qrotu"] = df.apply(update_Qrot, arguments=[df.viblvl_u, df.iso])
 
             Q, Qvib, Qrotu, Qrotl = df.Q, df.Qvib, df.Qrotu, df.Qrotl
 
         else:
-            iso = df.attrs["iso"]
+            iso = int(self.input.isotope)
 
             parsum = self.get_partition_function_calculator(molecule, iso, state)
 
@@ -2521,14 +2957,21 @@ def Qneq_Qvib_Qrotu_Qrotl(
             # dfQrot index is already 'viblvl'
             dfQrot_dict = dict(list(zip(dfQrot.index, dfQrot.Qrot)))
 
-            dg = df.loc[:]
+            if self.dataframe_type == "pandas":
+                dg = df.loc[:]
 
-            # Add lower state Qrot
-            dg_sorted = dg.set_index(["viblvl_l"], inplace=False)
-            df.loc[:, "Qrotl"] = dg_sorted.index.map(dfQrot_dict.get).values
-            # Add upper state energy
-            dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
-            df.loc[:, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
+                # Add lower state Qrot
+                dg_sorted = dg.set_index(["viblvl_l"], inplace=False)
+                df.loc[:, "Qrotl"] = dg_sorted.index.map(dfQrot_dict.get).values
+                # Add upper state energy
+                dg_sorted = dg.set_index(["viblvl_u"], inplace=False)
+                df.loc[:, "Qrotu"] = dg_sorted.index.map(dfQrot_dict.get).values
+            elif self.dataframe_type == "vaex":
+                dg = df
+                # Add lower state Qrot
+                df["Qrotl"] = df["viblvl_l"].apply(lambda x: dfQrot_dict.get(x))
+                # Add upper state energy
+                df["Qrotu"] = df["viblvl_u"].apply(lambda x: dfQrot_dict.get(x))
 
             Q, Qvib, Qrotu, Qrotl = Q, Qvib, df.Qrotu, df.Qrotl
 
@@ -2637,6 +3080,7 @@ def calc_populations_noneq(
                 df,
                 dbformat=self.params.dbformat,
                 lvlformat=self.params.levelsfmt,
+                dataframe_type=self.dataframe_type,
                 verbose=self.verbose,
             )
             assert "viblvl_u" in df
@@ -3049,6 +3493,7 @@ def calc_linestrength_noneq(self):
         """
 
         df = self.df1
+        df.attrs = self.df1.attrs
         Tref = self.input.Tref
 
         if len(df) == 0:
@@ -3068,7 +3513,10 @@ def calc_linestrength_noneq(self):
         nl = df.nl
         # ... remove Qref, nref, etc.
         # start from for tabulated linestrength
-        line_strength = df.int.copy()  # TODO: savememory; replace the "int" column
+        if self.dataframe_type == "pandas":
+            line_strength = df.int.copy()  # TODO: savememory; replace the "int" column
+        elif self.dataframe_type == "vaex":
+            line_strength = df.int
         if not self.molparam.terrestrial_abundances:
             raise NotImplementedError(
                 "Formula not corrected for non-terrestrial isotopic abundances"
@@ -3148,7 +3596,10 @@ def calc_emission_integral(self):
 
         A_ul = df["Aul"]  # (s-1)
 
-        DeltaE = cm2J(df.wav)  # (cm-1) -> (J)
+        if self.dataframe_type == "pandas":
+            DeltaE = cm2J(df.wav)  # (cm-1) -> (J)
+        elif self.dataframe_type == "vaex":
+            DeltaE = cm2J_vaex(df.wav)  # (cm-1) -> (J)
         Ei = n_ua * A_ul / 4 / pi * DeltaE  # (W/sr)
 
         Ei *= 1e3  # (W/sr) -> (mW/sr)
@@ -3214,26 +3665,32 @@ def _cutoff_linestrength(self, cutoff=None):
         # Estimate error being made:
         if self.warnings["LinestrengthCutoffWarning"] != "ignore":
 
-            error = df.S[b].sum() / df.S.sum() * 100
+            if self.dataframe_type == "vaex":
+                error_cutoff = df[b].sum(df[b].S) / df.sum(df.S) * 100
+            else:
+                error_cutoff = df.S[b].sum() / df.S.sum() * 100
 
             if verbose >= 2:
                 print(
                     "Discarded {0:.2f}% of lines (linestrength<{1}cm-1/(#.cm-2))".format(
-                        Nlines_cutoff / len(df.S) * 100, cutoff
+                        Nlines_cutoff / len(df) * 100, cutoff
                     )
-                    + " Estimated error: {0:.2f}%".format(error)
+                    + " Estimated error: {0:.2f}%".format(error_cutoff)
                 )
-            if error > self.misc.warning_linestrength_cutoff:
+            if error_cutoff > self.misc.warning_linestrength_cutoff:
                 self.warn(
                     "Estimated error after discarding lines is large: {0:.2f}%".format(
-                        error
+                        error_cutoff
                     )
                     + ". Consider reducing cutoff",
                     "LinestrengthCutoffWarning",
                 )
 
         try:
-            assert sum(~b) > 0
+            if self.dataframe_type == "pandas":
+                assert sum(~b) > 0
+            elif self.dataframe_type == "vaex":
+                assert (~b).sum() > 0
         except AssertionError as err:
             self.plot_linestrength_hist(cutoff=cutoff)
             raise AssertionError(
@@ -3245,7 +3702,12 @@ def _cutoff_linestrength(self, cutoff=None):
             ) from err
 
         # update df1:
-        self.df1 = pd.DataFrame(df[~b])
+        if self.dataframe_type == "pandas":
+            self.df1 = pd.DataFrame(df[~b])
+        elif self.dataframe_type == "vaex":
+            self.df1 = df[~b]
+            self.df1.attrs = df.attrs
+
         #        df.drop(b.index, inplace=True)   # performance: was not faster
         # ... @dev performance: quite long to select here, but I couldn't find a faster
         # ... alternative
@@ -3308,6 +3770,7 @@ def _reinitialize(self):
             # ... Operate on a duplicate dataframe to make it possible to do different
             # ... runs without reloading database
             self.df1 = self.df0.copy()
+            self.df1.attrs = self.df0.attrs
 
             # abundance and molar_mass should have been copied even if they are attributes
             # (only 1 molecule, 1 isotope) and not a column (line specific) in the database
@@ -3324,35 +3787,39 @@ def _reinitialize(self):
 
         self.profiler.start("memory_usage_warning", 3)
         # Check memory size
-        try:
-            # Retrieving total user RAM
-            mem = virtual_memory()
-            mem = mem.total  # total physical memory available
-            # Checking if object type column exists
-            if "O" in self.df1.dtypes.unique():
-                limit = mem / 25  # 4% of user RAM
-                self.warn(
-                    "'object' type column found in database, calculations and "
-                    + "memory usage would be faster with a numeric type. Possible "
-                    + "solution is to not use 'save_memory' and convert the columns to dtype.",
-                    "PerformanceWarning",
-                )
-            else:
-                limit = mem / 2  # 50 % of user RAM
-
-            # Note: the difference between deep=True and deep=False is around 4 times
-
-            df_size = self.df1.memory_usage(deep=False).sum()
-
-            if df_size > limit:
-                self.warn(
-                    "Line database is large: {0:.0f} Mb".format(df_size * 1e-6)
-                    + ". Consider using save_memory "
-                    + "option, if you don't need to reuse this factory to calculate new spectra",
-                    "MemoryUsageWarning",
-                )
-        except ValueError:  # had some unexplained ValueError: __sizeof__() should return >= 0
-            pass
+        if self.dataframe_type == "pandas":
+            try:
+                # Retrieving total user RAM
+                mem = virtual_memory()
+                mem = mem.total  # total physical memory available
+                # Checking if object type column exists
+                if "O" in self.df1.dtypes.unique():
+                    limit = mem / 25  # 4% of user RAM
+                    self.warn(
+                        "'object' type column found in database, calculations and "
+                        + "memory usage would be faster with a numeric type. Possible "
+                        + "solution is to not use 'save_memory' and convert the columns to dtype.",
+                        "PerformanceWarning",
+                    )
+                else:
+                    limit = mem / 2  # 50 % of user RAM
+
+                # Note: the difference between deep=True and deep=False is around 4 times
+                if self.dataframe_type == "pandas":
+                    df_size = self.df1.memory_usage(deep=False).sum()
+
+                    if df_size > limit:
+                        self.warn(
+                            "Line database is large: {0:.0f} Mb".format(df_size * 1e-6)
+                            + ". Consider using save_memory "
+                            + "option, if you don't need to reuse this factory to calculate new spectra",
+                            "MemoryUsageWarning",
+                        )
+                elif self.dataframe_type == "vaex":
+                    if self.verbose:
+                        print("Using Vaex to save memory")
+            except ValueError:  # had some unexplained ValueError: __sizeof__() should return >= 0
+                pass
 
         self.profiler.stop("memory_usage_warning", "Check Memory usage of database")
 
diff --git a/radis/lbl/broadening.py b/radis/lbl/broadening.py
index 93df61566..5c0ff3f8d 100644
--- a/radis/lbl/broadening.py
+++ b/radis/lbl/broadening.py
@@ -64,6 +64,7 @@
 from warnings import warn
 
 import numpy as np
+import vaex
 from numba import float64, jit
 from numpy import arange, exp
 from numpy import log as ln
@@ -234,7 +235,8 @@ def gaussian_FT(w_centered, hwhm):
 
     Returns
     -------
-    array
+    numpy array or vaex expression.
+    Depends upon the type of hwhm parameter, if it of vaex expression type then,returned value is vaex expression
 
     See Also
     --------
@@ -293,7 +295,8 @@ def pressure_broadening_HWHM(
 
     Returns
     -------
-    pandas Series        [shape N]
+    pandas Series or Vaex Expression        [shape N]
+        Depends on the engine used . Vaex Expression are memory efficient
         Lorentzian half-width at half-maximum (HWHM) for each line profile
 
     References
@@ -896,6 +899,7 @@ def _calc_broadening_HWHM(self):
 
         # diluent and their broadening coeff dictionary
         diluent_broadening_coeff = {}
+
         for key in diluent:
             if key == "air":
                 # no need to add broadening dictionary with air, as "airbrd" and "Tdpair" is already in the dataframe
@@ -1076,7 +1080,13 @@ def _add_voigt_broadening_HWHM(
         """
 
         # Check self broadening is here
-        if not "Tdpsel" in list(df.keys()):
+
+        if self.dataframe_type == "pandas":
+            columns = list(df.keys())
+        elif self.dataframe_type == "vaex":
+            columns = df.column_names
+
+        if not "Tdpsel" in columns:
             self.warn(
                 "Self-broadening temperature coefficient Tdpsel not given in database: used Tdpair instead",
                 "MissingSelfBroadeningWarning",
@@ -1088,7 +1098,7 @@ def _add_voigt_broadening_HWHM(
 
         molar_mass = self.get_molar_mass(df)
 
-        if not "selbrd" in list(df.keys()):
+        if not "selbrd" in columns:
             self.warn(
                 "Self-broadening coefficient selbrd not given in database: used airbrd instead",
                 "MissingSelfBroadeningWarning",
@@ -1155,7 +1165,11 @@ def _add_collisional_broadening_HWHM(
         """
 
         # Check self broadening temperature-dependance coefficient is here
-        if not "Tdpsel" in list(df.keys()):
+        if self.dataframe_type == "pandas":
+            columns = list(df.keys())
+        elif self.dataframe_type == "vaex":
+            columns = df.column_names
+        if not "Tdpsel" in columns:
             self.warn(
                 "Self-broadening temperature coefficient `Tdpsel` not given in database: used `Tdpair` instead",
                 "MissingSelfBroadeningTdepWarning",
@@ -1166,7 +1180,7 @@ def _add_collisional_broadening_HWHM(
             Tdpsel = df.Tdpsel
 
         # Check self broadening is here
-        if not "selbrd" in list(df.keys()):
+        if not "selbrd" in columns:
             self.warn(
                 "Self-broadening reference width `selbrd` not given in database: used air broadening reference width `airbrd` instead",
                 "MissingSelfBroadeningWarning",
@@ -1246,7 +1260,10 @@ def _collisional_lineshape(self, dg, wbroad_centered):
 
         # Prepare vectorized operations:
         try:  # make it a (1, N) row vector
-            gamma_lb = gamma_lb.values.reshape((1, -1))
+            if isinstance(gamma_lb, np.ndarray) or self.dataframe_type == "pandas":
+                gamma_lb = gamma_lb.values.reshape((1, -1))
+            elif self.dataframe_type == "vaex":
+                gamma_lb = gamma_lb.to_numpy().reshape((1, -1))
         except AttributeError:  # probably not a dataframe: assert there is one line only.
             assert type(gamma_lb) is np.float64
 
@@ -1308,7 +1325,10 @@ def _gaussian_lineshape(self, dg, wbroad_centered):
 
         # Prepare vectorized operations:
         try:  # make it a (1,N) row vector
-            gamma_db = gamma_db.values.reshape((1, -1))
+            if self.dataframe_type == "pandas":
+                gamma_db = gamma_db.values.reshape((1, -1))
+            elif self.dataframe_type == "vaex":
+                gamma_db = gamma_db.to_numpy().reshape((1, -1))
         except AttributeError:  # probably not a dataframe: assert there is one line only.
             assert type(gamma_db) is np.float64
 
@@ -1374,12 +1394,17 @@ def _voigt_broadening(self, dg, wbroad_centered, jit=True):
 
         # Prepare coefficients, vectorize
         # ... get Voigt HWHM
-        if not "hwhm_voigt" in list(dg.keys()):
+        if self.dataframe_type == "vaex":
+            columns = dg.column_names
+        elif self.dataframe_type == "pandas":
+            columns = list(dg.keys())
+
+        if not "hwhm_voigt" in columns:
             raise KeyError(
                 "hwhm_voigt: Calculate broadening with "
                 + "calc_voigt_broadening_HWHM first"
             )
-        if not "hwhm_lorentz" in list(dg.keys()):
+        if not "hwhm_lorentz" in columns:
             raise KeyError(
                 "hwhm_lorentz: Calculate broadening with "
                 + "calc_voigt_broadening_HWHM first"
@@ -1390,11 +1415,16 @@ def _voigt_broadening(self, dg, wbroad_centered, jit=True):
 
         # Prepare vectorized operations: make (1, N) row vectors
         try:
-            hwhm_lorentz = hwhm_lorentz.values.reshape((1, -1))
-            hwhm_voigt = hwhm_voigt.values.reshape((1, -1))
+            if self.dataframe_type == "pandas":
+                hwhm_lorentz = hwhm_lorentz.values.reshape((1, -1))
+                hwhm_voigt = hwhm_voigt.values.reshape((1, -1))
+            elif self.dataframe_type == "vaex":
+                hwhm_lorentz = hwhm_lorentz.to_numpy().reshape((1, -1))
+                hwhm_voigt = hwhm_voigt.to_numpy().reshape((1, -1))
         except AttributeError:  # probably not a dataframe: assert there is one line only.
-            assert type(hwhm_lorentz) is np.float64
-            assert type(hwhm_voigt) is np.float64
+            if self.dataframe_type == "pandas":
+                assert type(hwhm_lorentz) is np.float64
+                assert type(hwhm_voigt) is np.float64
 
         # Calculate broadening for all lines
         # ----------------------------------
@@ -1463,8 +1493,12 @@ def _calc_lineshape(self, dg):
         shifted_wavenum = dg.shiftwav
 
         try:  # make it a row vector
-            shifted_wavenum = shifted_wavenum.values.reshape((1, -1))
-            N = len(dg)
+            if self.dataframe_type == "pandas":
+                shifted_wavenum = shifted_wavenum.values.reshape((1, -1))
+                N = len(dg)
+            elif self.dataframe_type == "vaex":
+                shifted_wavenum = shifted_wavenum.to_numpy().reshape((1, -1))
+                N = len(dg)
         except AttributeError:  # probably not a dataframe: one line only.
             assert type(shifted_wavenum) is np.float64
             N = 1
@@ -1579,8 +1613,12 @@ def _init_w_axis(w_dat, log_p):
         log_pL = self.params.dxL  # LDM user params
         log_pG = self.params.dxG  # LDM user params
 
-        wL_dat = df.hwhm_lorentz.values * 2  # FWHM
-        wG_dat = df.hwhm_gauss.values * 2  # FWHM
+        if self.dataframe_type == "pandas":
+            wL_dat = df.hwhm_lorentz.values * 2.000  # FWHM
+            wG_dat = df.hwhm_gauss.values * 2.000  # FWHM
+        elif self.dataframe_type == "vaex":
+            wL_dat = df.hwhm_lorentz * 2.000  # FWHM
+            wG_dat = df.hwhm_gauss * 2.000  # FWHM
 
         wL = _init_w_axis(wL_dat, log_pL)  # FWHM
         self.NwL = len(wL)
@@ -1777,7 +1815,13 @@ def _apply_lineshape(self, broadened_param, line_profile, shifted_wavenum):
         )  # calculation vector of wavenumbers (shape W + space B on the sides)
 
         # Vectorize the chunk of lines
-        S = broadened_param.reshape((1, -1))
+        if isinstance(broadened_param, np.ndarray):
+            S = broadened_param.reshape((1, -1))
+        elif self.dataframe_type == "vaex":
+            S = broadened_param.to_numpy().reshape((1, -1))
+        else:
+            raise NotImplementedError
+
         shifted_wavenum = shifted_wavenum.reshape((1, -1))  # make it a row vector
 
         # Get truncation array
@@ -1970,9 +2014,16 @@ def _apply_lineshape_LDM(
 
         # LDM : Next calculate how the line is distributed over the 2x2x2 bins.
         ki0, ki1, tvi = self._get_indices(shifted_wavenum, wavenumber_calc)
-        li0, li1, tGi = self._get_indices(np.log(wG_dat), np.log(wG))
-        mi0, mi1, tLi = self._get_indices(np.log(wL_dat), np.log(wL))
 
+        if self.dataframe_type == "vaex":
+            # TODO : check memory usage & if there is no other implementation without numpy conversion
+            wG_dat = wG_dat.to_numpy()
+            wL_dat = wL_dat.to_numpy()
+            li0, li1, tGi = self._get_indices(np.log(wG_dat), np.log(wG))
+            mi0, mi1, tLi = self._get_indices(np.log(wL_dat), np.log(wL))
+        elif self.dataframe_type == "pandas":
+            li0, li1, tGi = self._get_indices(np.log(wG_dat), np.log(wG))
+            mi0, mi1, tLi = self._get_indices(np.log(wL_dat), np.log(wL))
         # Next assign simple weights:
         if optimization == "min-RMS":
 
@@ -2415,7 +2466,15 @@ def _broaden_lines(self, df):
                             )
                         )
 
-                    for i, (_, dg) in enumerate(df.groupby(arange(len(df)) % N)):
+                    # Cut lines in smaller bits for better memory handling
+                    if self.dataframe_type == "pandas":
+                        df = df.groupby(arange(len(df)) % N)
+                    elif self.dataframe_type == "vaex":
+                        # idx is added to as vaex don't have index column . Then df is divided into chunks
+                        df["idx"] = vaex.vrange(0, len(df)) % N
+                        df = df.groupby(df.idx)
+
+                    for i, (_, dg) in enumerate(df):
                         line_profile = self._calc_lineshape(dg)
                         (wavenumber, absorption) = self._apply_lineshape(
                             dg.S.values, line_profile, dg.shiftwav.values
@@ -2429,7 +2488,16 @@ def _broaden_lines(self, df):
                     self.reftracker.add(doi["DIT-2020"], "algorithm")
                     # Iterating over the chunks of the line database
                     # Using DIT Algorithm calculations for optimized loops
-                    for i, (_, dg) in enumerate(df.groupby(arange(len(df)) % N)):
+
+                    # Cut lines in smaller bits for better memory handlin
+                    if self.dataframe_type == "pandas":
+                        df = df.groupby(arange(len(df)) % N)
+                    elif self.dataframe_type == "vaex":
+                        # idx is added to as vaex don't have index column . Then df is divided into chunks
+                        df["idx"] = vaex.vrange(0, len(df)) % N
+                        df = df.groupby(df.idx)
+
+                    for i, (_, dg) in enumerate(df):
                         (
                             line_profile_LDM,
                             wL_i,
@@ -2760,7 +2828,7 @@ def _find_weak_lines(self, weak_rel_intensity_threshold):
         self.profiler.start("weak_lines", 2, "... classifying lines as weak or strong")
         # Get approximate spectral absorption coefficient
         rough_spectrum, S_density_on_grid, line2grid_proj_left = project_lines_on_grid(
-            df, wavenumber_calc, wstep
+            df, wavenumber_calc, wstep, self.dataframe_type
         )
 
         #     :
@@ -2873,7 +2941,7 @@ def calculate_pseudo_continuum(self, noneq=False):
             # Calculate continuum
             if noneq:
                 k_continuum, j_continuum, _, _, _ = project_lines_on_grid_noneq(
-                    df_weak_lines, wavenumber_calc, wstep
+                    df_weak_lines, wavenumber_calc, wstep, self.dataframe_type
                 )
 
                 if __debug__:
@@ -2896,7 +2964,7 @@ def calculate_pseudo_continuum(self, noneq=False):
 
             else:
                 k_continuum, _, _ = project_lines_on_grid(
-                    df_weak_lines, wavenumber_calc, wstep
+                    df_weak_lines, wavenumber_calc, wstep, self.dataframe_type
                 )
 
                 if __debug__:
@@ -2990,7 +3058,7 @@ def _add_pseudo_continuum(self, abscoeff_v, k_continuum):
         return abscoeff_v
 
 
-def project_lines_on_grid(df, wavenumber, wstep):
+def project_lines_on_grid(df, wavenumber, wstep, dataframe_type="pandas"):
     """Quickly sums all lines on wavespace grid as rectangles of HWHM
     corresponding to ``hwhm_voigt`` and a spectral absorption coefficient value so
     that linestrength is conserved.
@@ -3021,9 +3089,15 @@ def project_lines_on_grid(df, wavenumber, wstep):
         ``wavenumber``. Size ``N``
     """
 
-    shiftwav = df.shiftwav.values  # cm-1  ,   size N (number of lines)
-    S = df.S.values  # cm/#  ~   cm-1/(#.cm-2)  ,   size N
-    wv = df.hwhm_voigt.values * 2  # HWHM > FWHM
+    if dataframe_type == "pandas":
+        shiftwav = df.shiftwav.values  # cm-1  ,   size N (number of lines)
+        S = df.S.values  # cm/#  ~   cm-1/(#.cm-2)  ,   size N
+        wv = df.hwhm_voigt.values * 2  # HWHM > FWHM
+    elif dataframe_type == "vaex":
+        # TODO : check memory usage & if there is no other implementation without numpy conversion
+        shiftwav = df.shiftwav.to_numpy()  # cm-1  ,   size N (number of lines)
+        S = df.S.to_numpy()  # cm/#  ~   cm-1/(#.cm-2)  ,   size N
+        wv = df.hwhm_voigt.to_numpy() * 2  # HWHM > FWHM
 
     # ... First get closest matching line (left, and right):
     iwav_on_grid_left = np.searchsorted(wavenumber, shiftwav.T, side="left").ravel() - 1
@@ -3124,7 +3198,7 @@ def rough_sum_on_grid():
     return k_rough_spectrum, S_density_on_grid, line2grid_projection_left
 
 
-def project_lines_on_grid_noneq(df, wavenumber, wstep):
+def project_lines_on_grid_noneq(df, wavenumber, wstep, dataframe_type="pandas"):
     """Quickly sums all lines on wavespace grid as rectangles of HWHM
     corresponding to ``hwhm_voigt`` and a spectral absorption coefficient value so
     that linestrength is conserved.
@@ -3170,10 +3244,17 @@ def project_lines_on_grid_noneq(df, wavenumber, wstep):
     be recomputed from the linestrength).
     """
 
-    shiftwav = df.shiftwav.values  # cm-1  ,   size N (number of lines)
-    S = df.S.values  # cm/#  ~   cm-1/(#.cm-2)  ,   size N
-    Ei = df.Ei.values  # mW/cm3/sr
-    wv = df.hwhm_voigt.values * 2  # HWHM > FWHM
+    if dataframe_type == "pandas":
+        shiftwav = df.shiftwav.values  # cm-1  ,   size N (number of lines)
+        S = df.S.values  # cm/#  ~   cm-1/(#.cm-2)  ,   size N
+        Ei = df.Ei.values  # mW/cm3/sr
+        wv = df.hwhm_voigt.values * 2  # HWHM > FWHM
+    elif dataframe_type == "vaex":
+        # TODO : check memory usage & if there is no other implementation without numpy conversion
+        shiftwav = df.shiftwav.to_numpy()  # cm-1  ,   size N (number of lines)
+        S = df.S.to_numpy()  # cm/#  ~   cm-1/(#.cm-2)  ,   size N
+        Ei = df.Ei.to_numpy()  # mW/cm3/sr
+        wv = df.hwhm_voigt.to_numpy() * 2  # HWHM > FWHM
 
     # ... First get closest matching line (left, and right):
     iwav_on_grid_left = np.searchsorted(wavenumber, shiftwav.T, side="left").ravel() - 1
diff --git a/radis/lbl/calc.py b/radis/lbl/calc.py
index cc4b4d04e..29c6ac7fc 100644
--- a/radis/lbl/calc.py
+++ b/radis/lbl/calc.py
@@ -29,6 +29,7 @@
     from radis.lbl.factory import SpectrumFactory
     from radis.lbl.base import get_wavenumber_range
 
+from radis import config
 from radis.misc.basics import all_in
 from radis.misc.utils import Default
 from radis.spectrum.spectrum import Spectrum
@@ -290,6 +291,8 @@ def calc_spectrum(
                 s, sf = calc_spectrum(..., return_factory=True)
                 sf.df1  # see the lines calculated
                 sf.eq_spectrum(...)  #  new calculation without reloading the database
+    engine : string
+        Vaex or Pandas . Default Pandas, if engine is vaex memory performance is improved
     **kwargs: other inputs forwarded to SpectrumFactory
         For instance: ``warnings``.
         See :py:class:`~radis.lbl.factory.SpectrumFactory` documentation for more
@@ -530,6 +533,7 @@ def _check_molecules_are_consistent(
         diluent_for_this_molecule = diluents_for_molecule(
             mole_fraction, diluent, molecule
         )
+        engine = config["DATAFRAME_ENGINE"]
 
         generated_spectrum = _calc_spectrum_one_molecule(
             wavenum_min=wavenum_min,
@@ -561,6 +565,7 @@ def _check_molecules_are_consistent(
             export_lines=export_lines,
             return_factory=return_factory,
             diluent=diluent_for_this_molecule,
+            engine=engine,
             **kwargs_molecule,
         )
 
@@ -620,6 +625,7 @@ def _calc_spectrum_one_molecule(
     export_lines,
     return_factory=False,
     diluent="air",
+    engine="pandas",
     **kwargs,
 ) -> Spectrum:
     """See :py:func:`~radis.lbl.calc.calc_spectrum`
@@ -719,6 +725,7 @@ def _is_at_equilibrium():
 
     # Load databank
     # -------------
+    sf.dataframe_type = engine
 
     # Get databank
     if (
diff --git a/radis/lbl/factory.py b/radis/lbl/factory.py
index 1acaf550f..8e9ce1e4c 100644
--- a/radis/lbl/factory.py
+++ b/radis/lbl/factory.py
@@ -96,6 +96,7 @@
     from radis.lbl.bands import BandFactory
     from radis.lbl.base import get_wavenumber_range
 
+from radis import config
 from radis.misc.basics import flatten, is_float, is_range, list_if_float, round_off
 from radis.misc.utils import Default
 from radis.phys.constants import k_b
@@ -534,6 +535,10 @@ def __init__(
         self.input.Tref = convert_and_strip_units(Tref, u.K)
         self.input.pressure_mbar = convert_and_strip_units(pressure, u.bar) * 1e3
         self.input.mole_fraction = mole_fraction
+        if config["DATAFRAME_ENGINE"] in ["pandas", "vaex"]:
+            self.dataframe_engine = config["DATAFRAME_ENGINE"]
+        else:
+            raise NotImplementedError
 
         self.input.path_length = convert_and_strip_units(path_length, u.cm)
         self.input.molecule = (
diff --git a/radis/lbl/loader.py b/radis/lbl/loader.py
index 4df88185d..e9b939401 100644
--- a/radis/lbl/loader.py
+++ b/radis/lbl/loader.py
@@ -57,6 +57,7 @@
 
 import numpy as np
 import pandas as pd
+import vaex
 
 from radis import config
 from radis.api.cdsdapi import cdsd2df
@@ -1094,6 +1095,9 @@ def fetch_databank(
         if not molecule:
             raise ValueError("Please define `molecule=` so the database can be fetched")
 
+        output = config["DATAFRAME_ENGINE"]
+        self.dataframe_type = output
+
         if include_neighbouring_lines:
             wavenum_min = self.params.wavenum_min_calc
             wavenum_max = self.params.wavenum_max_calc
@@ -1105,6 +1109,7 @@ def fetch_databank(
         # and saved in output spectra information
         self.params.dbformat = dbformat
         self.misc.load_energies = load_energies
+        self.dataframe_type = output
         self.levels = levels
         if levels is not None:
             self.levelspath = ",".join([format_paths(lvl) for lvl in levels.values()])
@@ -1171,12 +1176,15 @@ def fetch_databank(
                     verbose=self.verbose,
                     return_local_path=True,
                     engine=memory_mapping_engine,
+                    output=output,
                     parallel=parallel,
                     extra_params=extra_params,
                 )
                 self.params.dbpath = ",".join(local_paths)
 
-                # ... explicitly write all isotopes based on isotopes found in the database
+                df.attrs = {}
+
+                # ... explicitely write all isotopes based on isotopes found in the database
                 if isotope == "all":
                     self.input.isotope = ",".join(
                         [str(k) for k in self._get_isotope_list(df=df)]
@@ -1216,9 +1224,15 @@ def fetch_databank(
                     )
                 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
+                    if output == "pandas":
+                        for df in frames:
+                            assert "iso" in df.columns
+                        df = pd.concat(frames, ignore_index=True)  # reindex
+                    elif output == "vaex":
+                        for df in frames:
+                            assert "iso" in df.column_names
+
+                        df = vaex.concat(frames)
                 else:
                     df = frames[0]
                 self.params.dbpath = "fetched from hitran"
@@ -1256,6 +1270,7 @@ def fetch_databank(
                 verbose=self.verbose,
                 return_local_path=True,
                 engine=memory_mapping_engine,
+                output=output,
                 parallel=parallel,
             )
             self.params.dbpath = ",".join(local_paths)
@@ -1312,6 +1327,7 @@ def fetch_databank(
                     return_local_path=True,
                     return_partition_function=True,
                     engine=memory_mapping_engine,
+                    output=output,
                 )
                 # @dev refactor : have a DatabaseClass from which we load lines and partition functions
                 if len(df) > 0:
@@ -1327,21 +1343,50 @@ def fetch_databank(
                 )
             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:
-                    if "iso" not in df.columns:
-                        assert "iso" in df.attrs
-                        df["iso"] = df.attrs["iso"]
-                # Keep attributes:
-                from radis.misc.basics import intersect
-
-                attrs = frames[0].attrs
-                for df in frames[1:]:
-                    attrs = intersect(attrs, df.attrs)
-                del attrs["iso"]  # added as a column (different for each line)
-                # Merge:
-                df = pd.concat(frames, ignore_index=True)  # reindex
-                df.attrs = attrs
-                self.params.dbpath = ",".join(local_paths)
+                if output == "pandas":
+                    for df in frames:
+                        if "iso" not in df.columns:
+                            # to save some memory, we have added "iso" as an attribute rather than a column if it is the same for all lines
+                            assert "iso" in df.attrs
+                            df["iso"] = df.attrs["iso"]
+                    # Keep attributes:
+                    from radis.misc.basics import intersect
+
+                    attrs = frames[0].attrs
+                    for df in frames[1:]:
+                        attrs = intersect(attrs, df.attrs)
+                    del attrs["iso"]  # added as a column (different for each line)
+                    # Merge:
+                    df = pd.concat(frames, ignore_index=True)  # reindex
+                    df.attrs = attrs
+                    self.params.dbpath = ",".join(local_paths)
+                elif output == "vaex":
+                    Frames = []
+                    for df in frames:
+                        if "iso" not in df.columns:
+                            assert "iso" in df.attrs
+                            import vaex
+
+                            df["iso"] = vaex.vconstant(
+                                int(df.attrs["iso"]), length=df.length_unfiltered()
+                            )
+                            Frames.append(df)
+                    # Keep attributes:
+                    from radis.misc.basics import intersect
+
+                    attrs = Frames[0].attrs
+                    for df in Frames[1:]:
+                        attrs = intersect(attrs, df.attrs)
+                    del attrs["iso"]  # added as a column (different for each line)
+                    # Merge:
+                    import vaex
+
+                    df = vaex.concat(Frames)  # reindex
+                    df.attrs = attrs
+
+                    self.params.dbpath = ",".join(local_paths)
+                else:
+                    raise NotImplementedError(output)
             else:
                 df = frames[0]
                 self.params.dbpath = local_paths[0]
@@ -1376,6 +1421,7 @@ def fetch_databank(
                 verbose=self.verbose,
                 return_local_path=True,
                 engine=memory_mapping_engine,
+                output=output,
                 parallel=parallel,
             )
             self.params.dbpath = ",".join(local_paths)
@@ -1396,18 +1442,32 @@ def fetch_databank(
             )
 
         # Always sort line database by wavenumber (required to SPARSE_WAVERANGE mode)
-        df.sort_values("wav", ignore_index=True, inplace=True)
+        if output == "pandas":
+            df.sort_values("wav", kind="mergesort", ignore_index=True, inplace=True)
+        elif output == "vaex":
+            try:
+                attrs = df.attrs
+            except:
+                attrs = {}
+            df = df.sort("wav", ascending=True)
+            df.attrs = attrs  # It is required because dataframe returned by sort_values doesn't have attrs, so I have to add it again.
+        else:
+            raise NotImplementedError(output)
 
         # %% Post-processing of the line database
         # ------------------------------------
         # (note : this is now done in 'fetch_hitemp' before saving to the disk)
         # spectroscopic quantum numbers will be needed for nonequilibrium calculations, and line survey.
         if parse_local_global_quanta and "locu" in df and source != "geisa":
-            df = parse_local_quanta(df, molecule, verbose=self.verbose)
+            df = parse_local_quanta(
+                df, molecule, verbose=self.verbose, dataframe_type=output
+            )
         if (
             parse_local_global_quanta and "globu" in df and source != "geisa"
         ):  # spectroscopic quantum numbers will be needed for nonequilibrium calculations :
-            df = parse_global_quanta(df, molecule, verbose=self.verbose)
+            df = parse_global_quanta(
+                df, molecule, verbose=self.verbose, dataframe_type=output
+            )
 
         # Remove non numerical attributes
         if drop_non_numeric:
@@ -1448,8 +1508,7 @@ def fetch_databank(
                     + "in fetch_databank"
                 )
 
-        self._remove_unnecessary_columns(df)
-
+        self._remove_unecessary_columns(df, output)
         return
 
     def load_databank(
@@ -1572,6 +1631,9 @@ def load_databank(
         # %% Check inputs
         # ---------
 
+        output = config["DATAFRAME_ENGINE"]
+        self.dataframe_type = output
+
         # use radis default for calculations non equilibrium calculations
         # if the levelsfmt is not specified in the databank
         if levelsfmt is None:
@@ -1626,6 +1688,7 @@ def load_databank(
         # %% Load Line databases
         # ----------------------
         self._reset_references()  # bibliographic references
+        self.dataframe_type = output
 
         self.df0 = self._load_databank(
             path,
@@ -1635,6 +1698,7 @@ def load_databank(
             drop_columns=drop_columns,
             load_columns=load_columns,
             include_neighbouring_lines=include_neighbouring_lines,
+            output=output,
         )
         self.misc.total_lines = len(self.df0)  # will be stored in Spectrum metadata
 
@@ -2064,7 +2128,8 @@ def _check_line_databank(self):
         #        # Reset index
         #        #    (cost ~ 1 ms but is needed if the user manually edited the database
         #        #    in between the load_database() and the calculation command
-        self.df0.reset_index(inplace=True, drop=True)  # drop: don't add old index
+        if self.dataframe_type == "pandas":
+            self.df0.reset_index(inplace=True, drop=True)  # drop: don't add old index
         # Finally commented: code may crash if users edit the database manually
         # (ex: modify broadening coefficients) and forgot to reset the index,
         # but that's for advanced users anyway. The cost (time+dont know what
@@ -2117,7 +2182,8 @@ def _load_databank(
         drop_columns,
         load_columns,
         include_neighbouring_lines=True,
-    ) -> pd.DataFrame:
+        output="pandas",
+    ):
 
         """Loads all available database files and keep the relevant one.
         Returns a Pandas dataframe.
@@ -2216,8 +2282,10 @@ def _load_databank(
 
         # subroutine load_and_concat
         # --------------------------------------
-        def load_and_concat(files):
-            """Concatenate many files in RAM
+
+        def load_and_concat(files, output="pandas"):
+            """Contatenate many files in RAM
+
             Parameters
             ----------
             files: list of str
@@ -2247,6 +2315,12 @@ def load_and_concat(files):
                                 "Missing doi for CDSD-HITEMP. Use HITEMP-2010?",
                                 "MissingReferenceWarning",
                             )
+
+                        if self.dataframe_type == "pandas":
+                            engine = "pytables"
+                        elif self.dataframe_type == "vaex":
+                            engine = "vaex"
+
                         df = cdsd2df(
                             filename,
                             version="hitemp" if dbformat == "cdsd-hitemp" else "4000",
@@ -2256,7 +2330,8 @@ def load_and_concat(files):
                             drop_non_numeric=True,
                             load_wavenum_min=wavenum_min,
                             load_wavenum_max=wavenum_max,
-                            engine="pytables",
+                            engine=engine,
+                            output=output,
                         )
                         # TODO: implement load_columns
                     elif dbformat in ["hitran", "hitemp"]:
@@ -2268,6 +2343,12 @@ def load_and_concat(files):
                             self.reftracker.add(
                                 doi["HITEMP-2010"], "line database"
                             )  # [HITEMP-2010]_
+
+                        if self.dataframe_type == "pandas":
+                            engine = "pytables"
+                        elif self.dataframe_type == "vaex":
+                            engine = "vaex"
+
                         df = hit2df(
                             filename,
                             cache=db_use_cached,
@@ -2276,7 +2357,8 @@ def load_and_concat(files):
                             drop_non_numeric=True,
                             load_wavenum_min=wavenum_min,
                             load_wavenum_max=wavenum_max,
-                            engine="pytables",
+                            engine=engine,
+                            output=output,
                         )
                     elif dbformat in ["hdf5-radisdb", "hitemp-radisdb"]:
                         if dbformat == "hitemp-radisdb":
@@ -2288,6 +2370,12 @@ def load_and_concat(files):
                                 f"Missing doi reference for database used {filename}",
                                 "MissingReferenceWarning",
                             )
+
+                        if self.dataframe_type == "pandas":
+                            engine = "guess"
+                        elif self.dataframe_type == "vaex":
+                            engine = "vaex"
+
                         df = hdf2df(
                             filename,
                             columns=columns,
@@ -2299,7 +2387,8 @@ def load_and_concat(files):
                             else None,
                             load_wavenum_min=wavenum_min,
                             load_wavenum_max=wavenum_max,
-                            engine="guess",
+                            engine=engine,
+                            output=output,
                         )
                     elif dbformat in ["exomol"]:
                         # self.reftracker.add("10.1016/j.jqsrt.2020.107228", "line database")  # [ExoMol-2020]
@@ -2322,9 +2411,13 @@ def load_and_concat(files):
                     if col in drop_columns or (
                         drop_columns == "all" and col not in drop_all_but_these
                     ):
-                        del df[col]
-                        dropped.append(col)
-                if self.verbose >= 2 and len(dropped) > 0:
+                        if output == "pandas":
+                            del df[col]
+                            dropped.append(col)
+                        elif output == "vaex":
+                            df.drop(col, inplace=True)
+                if verbose >= 2 and len(dropped) > 0:
+
                     print("Dropped columns: {0}".format(dropped))
 
                 # Crop to the wavenumber of interest
@@ -2348,18 +2441,28 @@ def load_and_concat(files):
 
             # Finally: Concatenate all
             if frames == []:
-                df = (
-                    pd.DataFrame()
-                )  # a database empty error will be raised a few lines below
+                if output == "pandas":
+                    df = (
+                        pd.DataFrame()
+                    )  # a database empty error will be raised a few lines below
+                elif output == "vaex":
+                    raise EmptyDatabaseError()
+                else:
+                    raise NotImplementedError(output)
             else:
-                df = pd.concat(frames, ignore_index=True)  # reindex
+                if output == "pandas":
+                    df = pd.concat(frames, ignore_index=True)  # reindex
+                elif output == "vaex":
+
+                    df = vaex.concat(frames)
+                    df.attrs = {}
 
             return df
 
         # end subroutine load_and_concat
         # --------------------------------------
 
-        df = load_and_concat(database)
+        df = load_and_concat(database, output=output)
 
         # Final checks
 
@@ -2411,7 +2514,11 @@ def load_and_concat(files):
             # check no isotope shows 0 line in this range. Raise an warning if it
             # happens
             for k in isotope_list:
-                if not (sum(df.iso == k) > 0):
+                if self.dataframe_type == "pandas":
+                    iso_count = sum(df.iso == k)
+                elif self.dataframe_type == "vaex":
+                    iso_count = (df.iso == k).sum()
+                if not (iso_count > 0):
                     msg = (
                         "Reference databank ({0:.2f}-{1:.2f}cm-1)".format(
                             minwavdb, maxwavdb
@@ -2457,12 +2564,12 @@ def load_and_concat(files):
                 )
             )
 
-        self._remove_unnecessary_columns(df)
+        self._remove_unecessary_columns(df, output)
 
         return df
 
-    def _remove_unnecessary_columns(self, df):
-        """Remove unnecessary columns and add values as attributes
+    def _remove_unecessary_columns(self, df, output="pandas"):
+        """Remove unecessary columns and add values as attributes
 
         Returns
         -------
@@ -2478,9 +2585,17 @@ def _remove_unnecessary_columns(self, df):
                     + "in SpectrumFactory. Use radis.calc_spectrum, which "
                     + "calculates them independently then use MergeSlabs"
                 )
+            if output == "pandas":
+                df.drop("id", axis=1, inplace=True)
+            elif output == "vaex":
+                df.drop("id", inplace=True)
+            else:
+                raise NotImplementedError(output)
 
-            df.drop("id", axis=1, inplace=True)
             df_metadata.append("id")
+
+            if "attrs" not in df.columns:
+                df.attrs = {}
             df.attrs["id"] = id_set[0]
         else:
             assert "id" in df.attrs or "molecule" in df.attrs
@@ -2489,7 +2604,13 @@ def _remove_unnecessary_columns(self, df):
             isotope_set = df.iso.unique()
 
             if len(isotope_set) == 1:
-                df.drop("iso", axis=1, inplace=True)
+                if output == "pandas":
+                    df.drop("iso", axis=1, inplace=True)
+                elif output == "vaex":
+                    df.drop("iso", inplace=True)
+                else:
+                    raise NotImplementedError(output)
+
                 df_metadata.append("iso")
                 df.attrs["iso"] = isotope_set[0]
         else:
@@ -2520,7 +2641,7 @@ def _get_isotope_list(self, molecule=None, df=None):
         if df is None:
             isotope_list = self.input.isotope.split(",")
         else:  # get all isotopes in line database
-            isotope_list = set(df.iso)
+            isotope_list = set(df.iso.unique())
 
         return [int(k) for k in isotope_list]
 
diff --git a/radis/misc/arrays.py b/radis/misc/arrays.py
index 16f24c526..2558096c5 100644
--- a/radis/misc/arrays.py
+++ b/radis/misc/arrays.py
@@ -354,6 +354,10 @@ def anynan(a):
     return np.isnan(np.dot(a, a))
 
 
+def anynan_vaex(a):
+    return not (a.countnan() == 0)
+
+
 @numba.njit
 def first_nonnan_index(a):
     """Returns index of first non-nan value in ``a``
diff --git a/radis/phys/convert.py b/radis/phys/convert.py
index ec4b3e37b..5cec3d632 100644
--- a/radis/phys/convert.py
+++ b/radis/phys/convert.py
@@ -17,6 +17,7 @@
 """
 
 import numpy as np
+import vaex
 
 from radis.phys.air import air2vacuum, vacuum2air
 from radis.phys.constants import c, eV, h, hc_k, k_b
@@ -92,6 +93,12 @@ def cm2J(E):
     return (E * 100) * (h * c)
 
 
+def cm2J_vaex(E):
+    """cm-1 to J."""
+    _assertcm_vaex(E)
+    return (E * 100) * (h * c)
+
+
 def cm2K(E):
     """cm-1 to K.
 
@@ -332,7 +339,10 @@ def atm2bar(p_atm):
 
 
 def _magn(x):
-    return np.round((np.log10(np.abs(x))))
+    if isinstance(x, vaex.expression.Expression):
+        return x.abs().log10().round()
+    else:
+        return np.round((np.log10(np.abs(x))))
 
 
 def _assertK(E):
@@ -353,6 +363,15 @@ def _assertcm(E):
             print(("Warning. Input values may not be in cm-1", E, "cm-1?"))
 
 
+def _assertcm_vaex(E):
+    if E.abs().sum():  # check E != 0 for both floats and arrays
+        try:
+            m = _magn(E)
+            assert (((1 <= m) & (m <= 5))).unique() == [True]
+        except AssertionError:
+            print(("Warning. Input values may not be in cm-1", E, "cm-1?"))
+
+
 def _asserteV(E):
     if np.sum(np.abs(E)) != 0:  # check E != 0 for both floats and arrays
         try:


Test patch:
diff --git a/radis/test/benchmark/vaex_vs_pandas_performance.py b/radis/test/benchmark/vaex_vs_pandas_performance.py
new file mode 100644
index 000000000..0316caad0
--- /dev/null
+++ b/radis/test/benchmark/vaex_vs_pandas_performance.py
@@ -0,0 +1,125 @@
+import time
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from radis import calc_spectrum, config
+from radis.misc.progress_bar import ProgressBar
+
+
+def compare_vaex_pandas_time():
+    """
+    Compares the time performance of pandas and Vaex and generates a plot. This scripts takes several minutes to run.
+    This results shoud shown that vaex and pandas provide similar performances in term if speed.
+    Returns
+    -------
+    None.
+    """
+    time_list, timeC_list, lines_list = [], [], []
+    time_list_va, timeC_list_va, lines_list_va = [], [], []
+    wmin = 1000
+    steps = 5
+    wmax_arr = np.geomspace(10, 1000, steps)
+
+    initial_engine = config[
+        "DATAFRAME_ENGINE"
+    ]  # To make sure dataframe engine not changed after running this test
+    pb = ProgressBar(N=2 * steps)
+    for i, engine in enumerate(["vaex", "pandas"]):
+        config["DATAFRAME_ENGINE"] = engine
+        for j, w_range in enumerate(wmax_arr):
+            t0 = time.time()
+            s = calc_spectrum(
+                wmin,
+                wmin + w_range,  # cm-1
+                molecule="H2O",
+                isotope="1,2,3",
+                pressure=1.01325,  # bar
+                Tgas=1000,
+                mole_fraction=0.1,
+                databank="hitemp",  # or 'hitemp'
+                wstep="auto",
+                cutoff=1e-28,
+                verbose=0,
+            )
+            t1 = time.time()
+            if engine == "vaex":
+                timeC_list_va.append(s.conditions["calculation_time"])
+                lines_list_va.append(s.conditions["lines_calculated"])
+                time_list_va.append(t1 - t0)
+                # lines_list_va.append(s.conditions['lines_calculated']+s.conditions['lines_cutoff'])
+            else:
+                timeC_list.append(s.conditions["calculation_time"])
+                lines_list.append(s.conditions["lines_calculated"])
+                time_list.append(t1 - t0)
+                # lines_list.append(s.conditions['lines_calculated']+s.conditions['lines_cutoff'])
+            pb.update(i * steps + (j + 1))
+    plt.figure()
+    plt.plot(lines_list, time_list, "k", label="pandas total")
+    plt.plot(lines_list, timeC_list, "k--", label="pandas computation")
+    plt.plot(lines_list_va, time_list_va, "r", label="vaex total")
+    plt.plot(lines_list_va, timeC_list_va, "r--", label="vaex computation")
+    plt.ylabel("Time [s]")
+    plt.xlabel("Number of lines")
+    plt.legend()
+
+    config["DATAFRAME_ENGINE"] = initial_engine
+
+
+# Compare the memory performance of Pandas and Vaex
+def compare_pandas_vs_vaex_memory():
+    """
+    Compare memory usage of `engine="vaex"` and `engine="pandas"` in calc_spectrum.
+    Expected behavior is "vaex" using much less memory. This function takes tens of seconds to run.
+    Returns
+    -------
+    None.
+    """
+
+    import tracemalloc
+
+    initial_engine = config[
+        "DATAFRAME_ENGINE"
+    ]  # To make sure dataframe engine not changed after running this test
+    for engine in ["pandas", "vaex"]:
+        config["DATAFRAME_ENGINE"] = engine
+        tracemalloc.start()
+        s = calc_spectrum(
+            1000,
+            1500,  # cm-1
+            molecule="H2O",
+            isotope="1,2,3",
+            pressure=1.01325,  # bar
+            Tgas=1000,  # K
+            mole_fraction=0.1,
+            wstep="auto",
+            databank="hitemp",  # or 'hitemp', 'geisa', 'exomol'
+            verbose=0,
+        )
+        snapshot = tracemalloc.take_snapshot()
+        memory = tracemalloc.get_traced_memory()
+        tracemalloc.stop()
+
+        # Some raw outputs
+        print("\n******** Engine = {} ***********".format(engine))
+        print(
+            "Peak, current = {:.1e}, {:.1e} for {:} lines calculated".format(
+                *memory, s.conditions["lines_calculated"]
+            )
+        )
+
+        # More sophisticated
+        print("*** List of biggest objects ***")
+        top_stats = snapshot.statistics("lineno")
+        for rank, stat in enumerate(top_stats[:3]):
+            print("#{}".format(rank + 1))
+            print(stat)
+
+        # Clear for next engine in the loop
+        tracemalloc.clear_traces()
+
+    config["DATAFRAME_ENGINE"] = initial_engine
+
+
+compare_vaex_pandas_time()
+compare_pandas_vs_vaex_memory()
diff --git a/radis/test/lbl/test_factory.py b/radis/test/lbl/test_factory.py
index bbd187404..f40063693 100644
--- a/radis/test/lbl/test_factory.py
+++ b/radis/test/lbl/test_factory.py
@@ -27,6 +27,7 @@
 import pytest
 
 import radis
+from radis import config
 from radis.lbl.factory import SpectrumFactory
 from radis.misc.printer import printm
 from radis.test.utils import setup_test_line_databases
@@ -644,6 +645,140 @@ def test_all_spectrum_using_wstep_auto(verbose=True, plot=False, *args, **kwargs
     assert sf._wstep == "auto"
 
 
+def test_vaex_and_pandas_spectrum():
+    """Compares spectrum calculated using vaex and pandas are same.
+    CONVOLUTED_QUANTITIES and dataframe df is compared to ensure the spectrum is same.
+
+    Here are things that we ensure are same for both vaex and pandas that
+    -Spectra calculated with Vaex & Pandas are the same
+    -Spectra calculated with Vaex & Pandas using Database cutoff are the same
+    Added in https://github.com/radis/radis/pull/580"""
+
+    from radis import calc_spectrum
+
+    # computing spectrum in vaex dataframe format
+    config["DATAFRAME_ENGINE"] = "vaex"
+    s_vaex, factory_s_vaex = calc_spectrum(
+        1900,
+        2300,  # cm-1
+        molecule="CO",
+        isotope="1,2,3",
+        pressure=1.01325,  # bar
+        Tgas=700,  # K
+        # mole_fraction={'CO': 0.5, 'CO2': 0.3},
+        mole_fraction=0.1,
+        path_length=1,  # cm
+        databank="hitran",  # or 'hitemp', 'geisa', 'exomol'
+        # use_cached=False,
+        return_factory=True,
+    )
+
+    # computing spectrum in pandas dataframe format
+    config["DATAFRAME_ENGINE"] = "pandas"
+    s_pandas, factory_s_pandas = calc_spectrum(
+        1900,
+        2300,  # cm-1
+        molecule="CO",
+        isotope="1,2,3",
+        pressure=1.01325,  # bar
+        Tgas=700,  # K
+        mole_fraction=0.1,
+        path_length=1,  # cm
+        databank="hitran",  # or 'hitemp', 'geisa', 'exomol'
+        # use_cached=False,
+        return_factory=True,
+    )
+
+    assert np.allclose(
+        s_vaex.get("absorbance"), s_pandas.get("absorbance"), equal_nan=True
+    )
+
+    for column in factory_s_pandas.df1.columns:
+        assert np.all(
+            factory_s_vaex.df1[column].to_numpy() == factory_s_pandas.df1[column]
+        )
+
+    #%% Additional test with other options (cutoff, ...)
+    from radis import SpectrumFactory
+
+    config["DATAFRAME_ENGINE"] = "vaex"
+    sf_vaex = SpectrumFactory(
+        wavelength_min=4200,
+        wavelength_max=4500,
+        cutoff=1e-23,
+        molecule="CO2",
+    )
+    sf_vaex.fetch_databank("hitran")
+    s_vaex = sf_vaex.eq_spectrum(Tgas=2000)  # failing on the 08/03/2023 - minouHub
+
+    config["DATAFRAME_ENGINE"] = "pandas"
+    sf_pd = SpectrumFactory(
+        wavelength_min=4200,
+        wavelength_max=4500,
+        cutoff=1e-23,
+        molecule="CO2",
+    )
+    sf_pd.fetch_databank("hitran")
+    s_pd = sf_pd.eq_spectrum(Tgas=2000)
+    for column in sf_pd.df1.columns:
+        assert np.all(sf_vaex.df1[column].to_numpy() == sf_pd.df1[column])
+
+    assert np.allclose(s_vaex.get("absorbance"), s_pd.get("absorbance"), equal_nan=True)
+
+
+#%%
+def test_vaex_and_pandas_spectrum_noneq():
+    """Compares the Spectrum calculated under non-equilibrium conditions .
+    Comparison is made between CONVOLUTED_QUANTITIES and dataframe df.
+
+    Here are things that we ensure are same for both vaex and pandas that
+    -Spectra calculated with Vaex & Pandas are the same
+    -Spectra calculated with Vaex & Pandas using Database cutoff are the same
+    Added in https://github.com/radis/radis/pull/580"""
+
+    from radis import calc_spectrum
+
+    conditions = {
+        "wmin": 2300,
+        "wmax": 2320,
+        "molecule": "CO",
+        "isotope": 1,
+        "pressure": 1.01325,
+        "mole_fraction": 0.1,
+        "wstep": "auto",
+        "path_length": 1,
+        "databank": "hitran",
+        "verbose": 3,
+        "return_factory": True,
+    }
+    # Calculating spectrum using vae
+    config["DATAFRAME_TYPE"] = "vaex"
+    s, factory_s = calc_spectrum(
+        **conditions,
+        Tgas=700,  # K
+        Tvib=710,
+        Trot=710,
+    )
+
+    # Calculating spectrum using pandas
+    config["DATAFRAME_ENGINE"] = "pandas"
+    s1, factory_s1 = calc_spectrum(
+        **conditions,
+        Tgas=700,  # K
+        Tvib=710,
+        Trot=710,
+    )
+
+    import numpy as np
+
+    # Comparing different quantities
+    assert np.allclose(s.get("absorbance"), s1.get("absorbance"), equal_nan=True)
+
+    # Comparing different columns of dataframe df
+    for column in factory_s1.df1.columns:
+        assert np.all(factory_s1.df1[column] == factory_s.df1[column].to_numpy())
+
+
 # --------------------------
 if __name__ == "__main__":
 
diff --git a/radis/test/lbl/test_loader.py b/radis/test/lbl/test_loader.py
index 696412a68..bcb35253d 100644
--- a/radis/test/lbl/test_loader.py
+++ b/radis/test/lbl/test_loader.py
@@ -11,6 +11,7 @@
 import matplotlib.pyplot as plt
 import pytest
 
+from radis import config
 from radis.lbl import SpectrumFactory
 from radis.misc.printer import printm
 from radis.test.utils import getTestFile, setup_test_line_databases
@@ -222,13 +223,188 @@ def test_custom_abundance(verbose=True, plot=False, *args, **kwargs):
     assert s.compare_with((3 * s3.take("abscoeff")), "abscoeff", rtol=0.5e-2, plot=True)
 
 
+def test_vaex_and_pandas_dataframe_fetch_databank():
+    """
+    Compares dataframes fetched and manipulated with vaex with that of pandas.
+    Every column of dataframe in vaex format is compared with corresponding column
+    of dataframe in pandas format .
+    Here expected
+    - All columns are same for dataframe in pandas and vaex format
+    - Ensure dataframes are same across all the databanks available in RADIS
+    Added in https://github.com/radis/radis/pull/580
+    """
+    import numpy as np
+    import pandas
+    import vaex
+
+    from radis.lbl import SpectrumFactory
+
+    sf = SpectrumFactory(
+        1500,
+        2000,
+        molecule="CO",
+        isotope="1,2,3",
+        pressure=1.01325,  # bar
+        path_length=1,  # cm
+    )
+
+    def compare_dataframe(df_vaex, df_pandas, columns):
+        if_same = True
+        for column in columns:
+            if_same = if_same * np.all(df_vaex[column].to_numpy() == df_pandas[column])
+
+        return if_same
+
+    # HITRAN
+    # Comparing Vaex dataframe and Pandas dataframe fetched from HITRAN
+
+    # Fetching in vaex dataframe format
+    config["DATAFRAME_ENGINE"] = "vaex"
+    sf.fetch_databank("hitran", memory_mapping_engine="vaex", load_columns="all")
+    df1 = sf.df0
+    assert isinstance(df1, vaex.dataframe.DataFrameLocal)
+
+    # Fetching in Pandas dataframe format
+    config["DATAFRAME_ENGINE"] = "pandas"
+    sf.fetch_databank("hitran", memory_mapping_engine="vaex", load_columns="all")
+    df2 = sf.df0
+    assert isinstance(df2, pandas.core.frame.DataFrame)
+
+    # Comparing the dataframe fetched from HITEMP
+    assert compare_dataframe(df1, df2, df2.columns)
+
+    # HITEMP
+    # Comparing Vaex dataframe and Pandas dataframe fetched from HITEMP
+
+    # Fetching in vaex dataframe format
+    config["DATAFRAME_ENGINE"] = "vaex"
+    sf.fetch_databank("hitemp", memory_mapping_engine="vaex", load_columns="all")
+    df1 = sf.df0
+    assert isinstance(df1, vaex.dataframe.DataFrameLocal)
+
+    # Fetching in Pandas dataframe format
+    config["DATAFRAME_ENGINE"] = "pandas"
+    sf.fetch_databank("hitemp", memory_mapping_engine="vaex", load_columns="all")
+    df2 = sf.df0
+    assert isinstance(df2, pandas.core.frame.DataFrame)
+
+    # Comparing the dataframe fetched from HITEMP
+    assert compare_dataframe(df1, df2, df2.columns)
+
+    # EXOMOL
+    # Comparing Vaex dataframe and Pandas dataframe fetched from EXOMOL
+
+    # Fetching in vaex dataframe format
+    config["DATAFRAME_ENGINE"] = "vaex"
+    sf.fetch_databank("exomol", memory_mapping_engine="vaex", load_columns="all")
+    df1 = sf.df0
+    assert isinstance(df1, vaex.dataframe.DataFrameLocal)
+
+    # Fetching in Pandas dataframe format
+    config["DATAFRAME_ENGINE"] = "pandas"
+    sf.fetch_databank("exomol", memory_mapping_engine="vaex", load_columns="all")
+    df2 = sf.df0
+    assert isinstance(df1, vaex.dataframe.DataFrameLocal)
+
+    # Comparing the dataframe fetched from EXOMOL
+    assert compare_dataframe(df1, df2, df2.columns)
+
+    # GEISA
+    # Comparing Vaex dataframe and Pandas dataframe fetched from GEISA
+
+    # Fetching in vaex dataframe format
+    config["DATAFRAME_ENGINE"] = "vaex"
+    sf.fetch_databank("geisa", memory_mapping_engine="vaex", load_columns="all")
+    df1 = sf.df0
+    assert isinstance(df1, vaex.dataframe.DataFrameLocal)
+
+    # Fetching in Pandas dataframe format
+    config["DATAFRAME_ENGINE"] = "pandas"
+    sf.fetch_databank("geisa", memory_mapping_engine="vaex", load_columns="all")
+    df2 = sf.df0
+    assert isinstance(df1, vaex.dataframe.DataFrameLocal)
+
+    # Comparing the dataframe fetched from GEISA
+    assert compare_dataframe(df1, df2, df2.columns)
+
+
+def test_vaex_and_pandas_dataframe_load_databank():
+    """
+    Compares dataframes fetched and manipulated with vaex with that of pandas.
+    Every column of dataframe in vaex format is compared with corresponding column
+    of dataframe in pandas format .
+    Added in https://github.com/radis/radis/pull/580
+    """
+    import numpy as np
+
+    from radis import SpectrumFactory
+
+    # Comparing HITRAN
+    sf = SpectrumFactory(
+        wavenum_min=2380,
+        wavenum_max=2400,
+        mole_fraction=400e-6,
+        path_length=100,  # cm
+        isotope=[1],
+    )
+
+    # Loading in vaex dataframe format
+    config["DATAFRAME_ENGINE"] = "vaex"
+    sf.load_databank("HITRAN-CO2-TEST")
+    df1 = sf.df0
+
+    # Loading in  pandas dataframe format
+    config["DATAFRAME_ENGINE"] = "pandas"
+    sf.load_databank("HITRAN-CO2-TEST")
+    df2 = sf.df0
+
+    # Comparing both the dataframes
+    for column in df2.columns:
+        assert np.all(df1[column].to_numpy() == df2[column])
+
+    # Comparing HITEMP
+    sf = SpectrumFactory(
+        2284.2,
+        2284.6,
+        wstep=0.001,  # cm-1
+        pressure=20 * 1e-3,  # bar
+        cutoff=0,
+        path_length=0.1,
+        mole_fraction=400e-6,
+        molecule="CO2",
+        isotope="1,2",
+        medium="vacuum",
+        truncation=5,
+        verbose=0,
+    )
+    sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
+
+    # Loading in vaex dataframe format
+    config["DATAFRAME_ENGINE"] = "vaex"
+    sf.load_databank("HITEMP-CO2-TEST", load_columns="all")
+    df1 = sf.df0
+
+    # Loading in pandas dataframe format
+    config["DATAFRAME_ENGINE"] = "pandas"
+    sf.load_databank("HITEMP-CO2-TEST", load_columns="all")
+    df2 = sf.df0
+
+    # Comparing both the dataframes
+    for column in df2.columns:
+        assert np.all(df1[column].to_numpy() == df2[column])
+
+
 def _run_testcases(verbose=True, plot=False):
 
     test_retrieve_from_database(plot=plot, verbose=verbose)
     test_ignore_cached_files()
     test_ignore_irrelevant_files(verbose=verbose)
     test_custom_abundance()
+    test_vaex_and_pandas_dataframe_fetch_databank()
+    test_vaex_and_pandas_dataframe_load_databank()
 
 
 if __name__ == "__main__":
     _run_testcases()
+
+# %%