Repository: radis/radis (language: python)

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

Add exomol support
<!-- Please be sure to check out our contributing guidelines,
https://github.com/radis/radis/blob/develop/CONTRIBUTING.md . -->

### Description
<!-- Provide a general description of what your pull request does. -->

First spectrum computed from ExoMol; borrowing some ExoJax codes from @HajimeKawahara  

```python
from radis import SpectrumFactory

sf = SpectrumFactory(
        4310,
        4320,
        molecule="H2O",
        isotope="1",
        verbose=3,
    )
sf.fetch_databank("exomol")   # recommended ExoMol database; or choose exomol_database=... to choose another
s = sf.eq_spectrum(Tgas=500)
s.plot()
```

![image](https://user-images.githubusercontent.com/16088743/126080498-5f437d32-c274-47bf-a17a-995ad7f3d4a7.png)

Still WIP,

- [x] read local ExoMol files
- [ ] convert them to RADIS's HDF5 format (Pytables #203 switching to Vaex #306 #295 )
- [x] load & convert to Radis's default Pandas DataFrame
- [x] compute missing values (Linestrengths, energies)
- [x] compute a Spectrum from ExoMol files
- [ ] add isotope abundances : [here](https://github.com/erwanp/radis/blob/add/exomol_support/radis/io/exomol.py#L391)
- [ ] add isotope number: [here](https://github.com/erwanp/radis/blob/add/exomol_support/radis/io/exomol.py#L590) and [there](https://github.com/erwanp/radis/blob/add/exomol_support/radis/io/exomol.py#L299)
- [ ] read molecular parameters from ExoMol (and update MolParams, currently = HITRAN data only. Requires #321 )
- [ ] multiple isotope support
- [ ] superline support

Interface : 
- [x] auto-download ExoMol files
- [x] auto-find ExoMol full names, and available databases
- [x] update SpectrumFactory fetch_databank
- [ ] update SpectrumFactory load_databank functions
- [ ] add explicit message to fill [KNOWN_EXOMOL_ISOTOPES_NAMES](https://github.com/radis/radis/blob/92fceb691e3695da25f8c78f3df28f487b5325d8/radis/db/classes.py#L280) if molecule not supported but exists in ExoMol
- [x] auto-find all ExoMol molecules by parsing [this page](https://exomol.com/data/molecules/).

Tests / Docs : 

- [x] add unitary tests
- [ ] add comparison with HITRAN/HITEMP
- [x] add a Gallery Example




<!-- If the pull request closes any open issues you can add this.
If you replace <Issue Number> with a number, GitHub will automatically link it.
If this pull request is unrelated to any issues, please remove
the following line. -->

Fixes #319


Discussion / hints:

New improvements in the ExoMol implementation, we are now able to compute non-HITRAN molecules (with some hardcoded parts, still). User-Interface is also as simple as possible : choose `databank='exomol'`  . See here a SiO2 spectrum from the EBJT Exomol : 

```python
from radis import calc_spectrum
s = calc_spectrum(1080, 1320,         # cm-1
                  molecule='SiO',
                  isotope='1',
                  pressure=1.01325,   # bar
                  Tgas=1000,           # K
                  mole_fraction=0.1,
                  path_length=1,      # cm
                  databank='exomol',  # or use ('exomol', "EBJT") for a specific database ("EBJT")
                  )
s.apply_slit(1, 'cm-1')       # simulate an experimental slit
s.plot('radiance')
```
![image](https://user-images.githubusercontent.com/16088743/126561745-7a35d249-78b9-48b1-8ff0-b4f17d04d329.png)

**First, good news:** nice comparison of H2CO in ExoMol and HITRAN. I will add it somewhere in the week.
![image](https://user-images.githubusercontent.com/36677294/126945708-166ed9e9-72a4-4ff0-97d8-6b83028bb56e.png)

**Now, help wanted @erwanp @gagan-aryan :**
I need to import other quantum numbers in the df0. To give you an example, the broadening of H2CO is function of J and K. This import of quantum numbers takes a huge memory (my 16 GB of RAM are full during the process). The memory overload happens here: https://github.com/radis/radis/blob/ed74864355c0ad4be21f34f06ec4cf854f6f2f1a/radis/io/exomolapi.py#L241

I plan to improve that. First:
- [ ] The number of lines should be diminished before this step. 

I really believe that I handle these data not efficiently. Otherwise, other ideas would be:
- [ ] add an argument in the fetch_database to let the user import the quantum number of his choice
- [ ] [your ideas]
This first implementation ExoMol implementation was borrowed from Exojax and used the Feather format. We'll switch to Vaex's HDF5 where your `n0 = 4 + dic_def['Landé'] + dic_def['lifetime']` definition can become a **virtual column**; i.e.; you'll use the same syntax to define the column, but the formula is evaluated only for the lines in your spectral range : no more memory problems.
This pull request **introduces 20 alerts** when merging a8b7b6518002a80188c9e870a6cd5b03b08b717f into e64b8e271715a26d13b92a91855b5eea5f96ba76 - [view on LGTM.com](https://lgtm.com/projects/g/radis/radis/rev/pr-b996f0c7c9bcca2fb14bcda7f23b77d31d150eec)

**new alerts:**

* 8 for Except block handles &#39;BaseException&#39;
* 3 for Unreachable code
* 2 for Testing equality to None
* 2 for Variable defined multiple times
* 2 for Wrong number of arguments in a call
* 1 for Mismatch in multiple assignment
* 1 for Missing call to \`\_\_init\_\_\` during object initialization
* 1 for Module is imported with &#39;import&#39; and &#39;import from&#39;
This pull request **introduces 7 alerts** and **fixes 1** when merging a85b43cb2e7f20abbbe60d5390a61dd06634ee73 into e64b8e271715a26d13b92a91855b5eea5f96ba76 - [view on LGTM.com](https://lgtm.com/projects/g/radis/radis/rev/pr-7da1373cc03bdd2534a9a94f2cc66a9cc98d9aaf)

**new alerts:**

* 6 for Except block handles &#39;BaseException&#39;
* 1 for Module is imported with &#39;import&#39; and &#39;import from&#39;

**fixed alerts:**

* 1 for Missing call to \`\_\_init\_\_\` during object initialization
Impossible to convert the trans.bz2 file of CH4 to a feather
**Code:**

```python
from radis import SpectrumFactory
molecule = "CH4"

conditions = {
    "mole_fraction": 0.01,
    "pressure": 1,
    "Tgas": 700,
}
# wmin, wmax = 2050, 2060
wmin, wmax = 2777, 2779

sf = SpectrumFactory(
        wmin=wmin,
        wmax=wmax,
        molecule=molecule,
        isotope="1",
        verbose=3,
        wstep = 2e-4,
    )
sf.fetch_databank("exomol", parfuncfmt = "exomol")
```

The file is just 1.7 GB! 

```python
MemoryError: Unable to allocate 2.85 GiB for an array with shape (2, 191003050) and data type float64
```
@gagan-aryan is currently working on writing a vaex-HDF5 builder for HITEMP #306 . We'll use the same to replace the Feather building step. 
About to merge all current implementation. There are many changes already and we don't want to diverge too much from the rest of the code.

Some things aren't final (see checkbox above) but will be finished in a later PR. 
# [Codecov](https://codecov.io/gh/radis/radis/pull/320?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) Report
> Merging [#320](https://codecov.io/gh/radis/radis/pull/320?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) (5322c49) into [develop](https://codecov.io/gh/radis/radis/commit/d91d7fdb5ccb4249335da57396b5686b7ace6398?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) (d91d7fd) will **decrease** coverage by `0.20%`.
> The diff coverage is `77.41%`.


```diff
@@             Coverage Diff             @@
##           develop     #320      +/-   ##
===========================================
- Coverage    74.86%   74.65%   -0.21%     
===========================================
  Files          135      142       +7     
  Lines        15120    15597     +477     
===========================================
+ Hits         11319    11644     +325     
- Misses        3801     3953     +152     
```

This pull request **introduces 7 alerts** and **fixes 1** when merging 5322c4973ba7ba6cf053bae74fabbfd63b450eed into d91d7fdb5ccb4249335da57396b5686b7ace6398 - [view on LGTM.com](https://lgtm.com/projects/g/radis/radis/rev/pr-3439ec018288a3975880fab812853eaa0a462f56)

**new alerts:**

* 6 for Except block handles &#39;BaseException&#39;
* 1 for Module is imported with &#39;import&#39; and &#39;import from&#39;

**fixed alerts:**

* 1 for Missing call to \`\_\_init\_\_\` during object initialization

Fix patch:
diff --git a/docs/examples/examples.rst b/docs/examples/examples.rst
index 220382572..52b002a0c 100644
--- a/docs/examples/examples.rst
+++ b/docs/examples/examples.rst
@@ -94,16 +94,16 @@ Calculate Partition Functions
 =============================
 
 By default and for equilibrium calculations, RADIS calculates Partition Functions
-using the TIPS program through [HAPI]_. These partition functions can be retrieved
+using the TIPS program ([TIPS-2020]_) through [HAPI]_. These partition functions can be retrieved
 with the :py:class:`~radis.levels.partfunc.PartFunc_Dunham` class::
 
-    from radis.levels.partfunc import PartFuncHAPI
+    from radis.levels.partfunc import PartFuncTIPS
     from radis.db.classes import get_molecule_identifier
 
     M = get_molecule_identifier('N2O')
     iso=1
 
-    Q = PartFuncHAPI(M, iso)
+    Q = PartFuncTIPS(M, iso)
     print(Q.at(T=1500))
 
 
diff --git a/docs/references/_bibliography.rst b/docs/references/_bibliography.rst
index c070e6183..70defe87a 100644
--- a/docs/references/_bibliography.rst
+++ b/docs/references/_bibliography.rst
@@ -19,7 +19,8 @@ List of bibliographic references used in this project:
 
 .. [Spectral-Synthesis-Algorithm] D.C.M. van den Bekerom, E. Pannier,
                 "A Discrete Integral Transform for Rapid Spectral Synthesis",
-                Journal of Quantitative Spectroscopy and Radiative Transfer (2021) `doi:10.1016/j.jqsrt.2020.107476 <https://www.sciencedirect.com/science/article/abs/pii/S0022407320310049>`__
+                Journal of Quantitative Spectroscopy and Radiative Transfer (2021)
+                `doi:10.1016/j.jqsrt.2020.107476 <https://www.sciencedirect.com/science/article/abs/pii/S0022407320310049>`__
 
 .. [Rothman-1998] L.S. Rothman, C.P. Rinsland, A. Goldman, S.T. Massie D.P. Edwards, J-M. Flaud,
                  A. Perrin, C. Camy-Peyret, V. Dana, J.-Y. Mandin, J. Schroeder, A. McCann,
@@ -29,6 +30,11 @@ List of bibliographic references used in this project:
                  Journal of Quantitative Spectroscopy and Radiative Transfer 60 (1998)
                  665 - 710, `doi:10.1016/S0022-4073(98)00078-8 <https://www.sciencedirect.com/science/article/abs/pii/S0022407398000788?via%3Dihub>`__
 
+.. [TIPS-2020] R.R Gamache, B. Vispoel, M. Rey, A. Nikitin, V. Tyuterev, O. Egorov, I.E Gordon, V. Boudon,
+                "Total internal partition sums for the HITRAN2020 database",
+                Journal of Quantitative Spectroscopy and Radiative Transfer 271 (2021)
+                `doi:10.1016/j.jqsrt.2021.107713 <https://www.sciencedirect.com/science/article/abs/pii/S0022407321002065?via%3Dihub>`__
+
 Line Databases
 --------------
 
diff --git a/examples/calculate_partition_functions_from_tips.py b/examples/calculate_partition_functions_from_tips.py
index 3c1c636f4..2ea7c0c4f 100644
--- a/examples/calculate_partition_functions_from_tips.py
+++ b/examples/calculate_partition_functions_from_tips.py
@@ -5,8 +5,8 @@
 =============================
 
 By default and for equilibrium calculations, RADIS calculates Partition Functions
-using the TIPS program through [HAPI]_. These partition functions can be retrieved
-with the :py:class:`~radis.levels.partfunc.PartFuncHAPI` class::
+using the TIPS program ([TIPS-2020]_) through [HAPI]_. These partition functions can be retrieved
+with the :py:class:`~radis.levels.partfunc.PartFuncTIPS` class::
 
 See Also
 --------
diff --git a/examples/plot_custom_abundances.py b/examples/plot_custom_abundances.py
new file mode 100644
index 000000000..caf3e7d37
--- /dev/null
+++ b/examples/plot_custom_abundances.py
@@ -0,0 +1,53 @@
+# -*- coding: utf-8 -*-
+"""
+=====================
+Use Custom Abundances
+=====================
+
+Custom isotopologue abundances can be defined, to model non-terrestrial
+atmospheres.
+
+Abundances are read and set directly in the :py:class:`~radis.lbl.factory.SpectrumFactory`
+using the :py:meth:`~radis.lbl.loader.DatabankLoader.get_abundance`
+and :py:meth:`~radis.lbl.loader.DatabankLoader.set_abundance` methods.
+
+Below, we compute a CO2 spectrum with custom abundances for the two main
+terrestrial isotopes (12C-16O2 ; 13C-16O2)
+
+See Also
+--------
+:py:class:`~radis.db.molparam.MolParams`
+
+"""
+
+from radis.test.utils import setup_test_line_databases
+
+setup_test_line_databases()  # creates "HITEMP-CO2-TEST" for this example
+
+from radis import SpectrumFactory
+
+sf = SpectrumFactory(
+    2284.2,
+    2284.6,
+    wstep=0.001,  # cm-1
+    pressure=20 * 1e-3,  # bar
+    mole_fraction=400e-6,
+    molecule="CO2",
+    isotope="1,2",
+    verbose=False,
+)
+sf.load_databank("HITEMP-CO2-TEST")
+#%%To explicitely identify the isotopes we can use the molparam attribute the Factory
+print(sf.molparam.get("CO2", 1, "isotope_name"))  # >> (12C)(16O)2
+print(sf.molparam.get("CO2", 2, "isotope_name"))  # >> (13C)(16O)2
+
+
+#%%Print the default abundance of the CO2 isotopes, compute a spectrum
+print("Abundance of CO2[1,2]", sf.get_abundance("CO2", [1, 2]))
+sf.eq_spectrum(2000).plot("abscoeff")
+
+#%% Set the abundance of CO2(626) to 0.8; and the abundance of CO2(636) to 0.2 (arbitrary):
+sf.set_abundance("CO2", [1, 2], [0.8, 0.2])
+print("New abundance of CO2[1,2]", sf.get_abundance("CO2", [1, 2]))
+
+sf.eq_spectrum(2000).plot("abscoeff", nfig="same")
diff --git a/examples/plot_linestrengths.py b/examples/plot_linestrengths.py
index d5e10ea71..74018a842 100644
--- a/examples/plot_linestrengths.py
+++ b/examples/plot_linestrengths.py
@@ -25,7 +25,7 @@
 
 from radis.db.classes import get_molecule, get_molecule_identifier
 from radis.io.hitemp import fetch_hitemp
-from radis.levels.partfunc import PartFuncHAPI
+from radis.levels.partfunc import PartFuncTIPS
 from radis.phys.constants import hc_k
 
 
@@ -33,11 +33,11 @@ def get_Qgas(molecule, iso, T):
 
     M = get_molecule_identifier(molecule)
 
-    Q = PartFuncHAPI(M, iso)
+    Q = PartFuncTIPS(M, iso)
     return Q.at(T=T)
 
 
-def calc_linestrength_eq(df, Tref, Tgas):
+def scale_linestrength_eq(df, Tref, Tgas):
 
     print("Scaling equilibrium linestrength")
 
@@ -78,7 +78,6 @@ def _calc_Q(molecule, iso, T_ref, T_gas):
     Tref = 296
     df = fetch_hitemp(
         molecule="CO",
-        databank_name="HITEMP",
         isotope="1, 2, 3",
         load_wavenum_min=2000,
         load_wavenum_max=2250,
@@ -86,7 +85,7 @@ def _calc_Q(molecule, iso, T_ref, T_gas):
 
     Tgas = 450
 
-    df = calc_linestrength_eq(df, Tref, Tgas)
+    df = scale_linestrength_eq(df, Tref, Tgas)
     plt.bar(df["wav"], df["S"])
     plt.xlabel("Wavenumbers in cm-1")
     plt.ylabel("Linestrengths in cm-1/(molecules/cm-2)")
diff --git a/radis/db/molparam.py b/radis/db/molparam.py
index bf94e34fe..fe907e6a7 100644
--- a/radis/db/molparam.py
+++ b/radis/db/molparam.py
@@ -145,33 +145,38 @@
 # even if HITRAN eventually changes the conventions and labels.
 
 
-class MolParams:
-    def __init__(self, file=None):
-        """Easy access to molecular parameters taken from HITRAN molparam.txt.
+class MolParams(object):
+    """Easy access to molecular parameters taken from HITRAN molparam.txt.
 
-        Parameters
-        ----------
-        file: str
-            if None the one in RADIS is taken
+    Parameters
+    ----------
+    file: str
+        if None the one in RADIS is taken
 
-        Examples
-        --------
-        Get earth abundance of CO2, isotope 1::
+    Examples
+    --------
+    Get earth abundance of CO2, isotope 1::
 
-            molpar = Molparams()
-            molpar.get(2, 1, 'abundance')         # 2 for CO2, 1 for isotope 1
+        molpar = Molparams()
+        molpar.get(2, 1, 'abundance')         # 2 for CO2, 1 for isotope 1
 
-        Note
-        ----
-        Isotope number was derived manually assuming the isonames were ordered in the database
-        The isotope name (ex: CO2 626) is kept for comparison if ever needed
+    .. minigallery:: radis.db.molparam.MolParams
+    :add-heading
 
-        References
-        ----------
+    Note
+    ----
+    Isotope number was derived manually assuming the isonames were ordered in the database
+    The isotope name (ex: CO2 626) is kept for comparison if ever needed
 
-        http://hitran.org/media/molparam.txt
-        """
+    References
+    ----------
+
+    http://hitran.org/media/molparam.txt
+    """
 
+    __slots__ = ["df", "terrestrial_abundances"]
+
+    def __init__(self, file=None, terrestrial_abundances=True):
         if file is None:
             file = getFile("molparam.txt")
 
@@ -181,6 +186,7 @@ def __init__(self, file=None):
         df["isotope_name_exomol"] = _add_exomol_name(df)
 
         self.df = df
+        self.terrestrial_abundances = terrestrial_abundances
 
         # ------
         try:  # Add hints (Python >3.6 only)
@@ -197,9 +203,22 @@ def get(self, M, I, key):
             molecule id
         I: int
             molecule isotope #
-        key: ``'abundance'``, ``'mol_mass'``, ``'isotope_name_exomol'``,
+        key: ``'abundance'``, ``'mol_mass'``, ``'isotope_name'``, ``'isotope_name_exomol'``
             parameter
 
+        Examples
+        --------
+        Get explicit name of an isotope::
+
+            from radis.db.molparam import MolParams
+            mp = MolParams()
+            print(mp.get("CO2", 1, "isotope_name"))   # >> (12C)(16O)2
+            print(mp.get("CO2", 2, "isotope_name"))   # >> (13C)(16O)2
+
+
+        .. minigallery:: radis.db.molparam.MolParams.get
+            :add-heading
+
         """
         try:
             float(M)
diff --git a/radis/lbl/bands.py b/radis/lbl/bands.py
index 20e36d11c..f583d3fb0 100644
--- a/radis/lbl/bands.py
+++ b/radis/lbl/bands.py
@@ -542,8 +542,8 @@ def non_eq_bands(
         # (Note: Emission Integral is non canonical quantity, equivalent to
         #  Linestrength for absorption)
         self.calc_populations_noneq(Tvib, Trot)
-        self._calc_linestrength_noneq()
-        self._calc_emission_integral()
+        self.calc_linestrength_noneq()
+        self.calc_emission_integral()
 
         # ----------------------------------------------------------------------
         # Cutoff linestrength
diff --git a/radis/lbl/base.py b/radis/lbl/base.py
index 922eced7f..a25521ed6 100644
--- a/radis/lbl/base.py
+++ b/radis/lbl/base.py
@@ -44,8 +44,8 @@
 - :py:meth:`radis.lbl.base.BaseFactory.calc_linestrength_eq`
 - :py:meth:`radis.lbl.base.BaseFactory.calc_populations_eq`
 - :py:meth:`radis.lbl.base.BaseFactory.calc_populations_noneq`
-- :py:meth:`radis.lbl.base.BaseFactory._calc_linestrength_noneq`
-- :py:meth:`radis.lbl.base.BaseFactory._calc_emission_integral`
+- :py:meth:`radis.lbl.base.BaseFactory.calc_linestrength_noneq`
+- :py:meth:`radis.lbl.base.BaseFactory.calc_emission_integral`
 - :py:meth:`radis.lbl.base.BaseFactory._cutoff_linestrength`
 
 Most methods are written in inherited class with the following inheritance scheme:
@@ -74,7 +74,6 @@
 
 # TODO: rename in get_molecule_name
 from radis.db.classes import get_molecule, get_molecule_identifier
-from radis.db.molparam import MolParams
 
 try:  # Proper import
     from .loader import KNOWN_LVLFORMAT, DatabankLoader, df_metadata
@@ -1598,8 +1597,8 @@ def _calc_degeneracies(self, df):
 
         return None  # dataframe updated directly
 
-    def get_abundance(self, df):
-        """Returns the isotopic abundance
+    def get_lines_abundance(self, df):
+        """Returns the isotopic abundance of each line in `df`
 
         Parameters
         ----------
@@ -1610,7 +1609,7 @@ def get_abundance(self, df):
         float or dict: The abundance of all the isotopes in the dataframe
         """
 
-        molpar = MolParams()
+        molpar = self.molparam
 
         if "id" in df.columns:
             id_set = df.id.unique()
@@ -1641,6 +1640,47 @@ def get_abundance(self, df):
             iso = df.attrs["iso"]
             return molpar.get(df.attrs["id"], iso, "abundance")
 
+    def get_molar_mass(self, df):
+        """Returns molar mass.
+
+        Parameters
+        ----------
+        df: dataframe
+
+        Returns
+        -------
+        The molar mass of all the isotopes in the dataframe
+        """
+        molpar = self.molparam
+
+        if "id" in df.columns:
+            raise NotImplementedError(">1 molecule")
+        elif "id" in df.attrs:
+            id = df.attrs["id"]
+        else:
+            # HARDCODED molar mass; for WIP ExoMol implementation, until MolParams
+            # is an attribute and can be updated with definitions from ExoMol.
+            # https://github.com/radis/radis/issues/321
+            HARCODED_MOLAR_MASS = {"SiO": {1: 43.971842}}
+            try:
+                return HARCODED_MOLAR_MASS[df.attrs["molecule"]][df.attrs["iso"]]
+            except KeyError:
+                raise NotImplementedError
+
+        if "iso" in df.columns:
+            iso_set = df.iso.unique()
+            molar_mass_dict = {}
+            for iso in iso_set:
+                molar_mass_dict[iso] = molpar.get(id, iso, "mol_mass")
+            molar_mass = df["iso"].map(molar_mass_dict)
+        else:
+            iso = df.attrs["iso"]
+            molar_mass = molpar.get(id, iso, "mol_mass")
+
+        return molar_mass
+
+        #
+
     def calc_weighted_trans_moment(self):
         """Calculate weighted transition-moment squared :math:`R_s^2` (in ``Debye^2``)
 
@@ -1725,7 +1765,11 @@ def calc_weighted_trans_moment(self):
         # Get moment
 
         # get abundance
-        abundance = self.get_abundance(df)
+        abundance = self.get_lines_abundance(df)
+        if not self.molparam.terrestrial_abundances:
+            raise NotImplementedError(
+                "Formula not corrected for non-terrestrial isotopic abundances"
+            )
 
         gl = df.gl
         El = df.El
@@ -1733,7 +1777,9 @@ def calc_weighted_trans_moment(self):
         Ia = abundance
         h = h_CGS  # erg.s
         c = c_CGS
-        S = df.int  # reference linestrength
+        S = (
+            df.int
+        )  # reference linestrength   ( computed with terrestrial isotopic abundances)
 
         weighted_trans_moment_sq = (
             (3 * h * c / 8 / pi ** 3)
@@ -1886,8 +1932,8 @@ def calc_einstein_coefficients(self):
     # calc_linestrength_eq
     # calc_populations_eq
     # calc_populations_noneq
-    # _calc_linestrength_noneq
-    # _calc_emission_integral
+    # calc_linestrength_noneq
+    # calc_emission_integral
     # _cutoff_linestrength
 
     # XXX =====================================================================
@@ -1961,6 +2007,11 @@ def calc_linestrength_eq(self, Tgas):
         Updates linestrength in df1. Cutoff criteria is applied afterwards.
 
         .. minigallery:: radis.lbl.base.BaseFactory.calc_linestrength_eq
+            :add-heading:
+
+        See Also
+        --------
+        :py:func:`~radis.lbl.base.linestrength_from_Einstein`
         """
 
         Tref = self.input.Tref
@@ -2025,6 +2076,51 @@ def _calc_Q(molecule, iso, state):
             molecule = get_molecule(df1.attrs["id"])
         state = self.input.state
 
+        # Partition functions
+        def Qgas(Tgas):
+            """Aggregate the values of Qgas
+
+            Parameters
+            ----------
+            Tgas: float (K)
+                gas temperature
+
+            Returns
+            -------
+            float or dict: Returns Qgas as a dictionary with isotope values as its keys
+
+            """
+            if "iso" in df1:
+                iso_set = df1.iso.unique()
+                if len(iso_set) == 1:
+                    self.warn(
+                        "There shouldn't be a Column 'iso' with a unique value",
+                        "PerformanceWarning",
+                    )
+                    iso = int(iso_set)
+                    parsum = self.get_partition_function_interpolator(
+                        molecule, iso, state
+                    )
+                    # TODO : use _calc_Q instead ? (which switches to partition function calculator if possible?)
+                    Q = parsum.at(Tgas)
+                    df1.attrs["Q"] = Q
+                    return Q
+                else:
+                    Qgas_dict = {}
+                    for iso in iso_set:
+                        parsum = self.get_partition_function_interpolator(
+                            molecule, iso, state
+                        )
+                        Qgas_dict[iso] = parsum.at(Tgas)
+                    return df1["iso"].map(Qgas_dict)
+
+            else:  # "iso" not in df:
+                iso = df1.attrs["iso"]
+                parsum = self.get_partition_function_interpolator(molecule, iso, state)
+                Q = parsum.at(Tgas)
+                df1.attrs["Q"] = Q
+                return Q
+
         def Qref_Qgas_ratio():
 
             if "iso" in df1:
@@ -2054,23 +2150,37 @@ def Qref_Qgas_ratio():
         # %% Calculate line strength at desired temperature
         # -------------------------------------------------
 
-        # This calculation is based on equation (A11) in Rothman 1998: "JQSRT, vol.
-        # 60, No. 5, pp. 665-710"
-        # An alternative strategy would be to calculate the linestrength from the
-        # Einstein A coefficient and the populations (see Klarenaar 2017 Eqn. 12)
-
-        # correct for Partition Function
-        df1["S"] = (
-            df1.int
-            * Qref_Qgas_ratio()
-            *
-            # ratio of Boltzman populations
-            exp(-hc_k * df1.El * (1 / Tgas - 1 / Tref))
-            *
-            # effect of stimulated emission
-            (1 - exp(-hc_k * df1.wav / Tgas))
-            / (1 - exp(-hc_k * df1.wav / Tref))
-        )  # [cm-1/(molecules/cm-2)]
+        if self.molparam.terrestrial_abundances:
+
+            # This calculation is based on equation (A11) in Rothman 1998: "JQSRT, vol.
+            # 60, No. 5, pp. 665-710"
+
+            # correct for Partition Function
+            df1["S"] = (
+                df1.int
+                * Qref_Qgas_ratio()
+                *
+                # ratio of Boltzman populations
+                exp(-hc_k * df1.El * (1 / Tgas - 1 / Tref))
+                *
+                # effect of stimulated emission
+                (1 - exp(-hc_k * df1.wav / Tgas))
+                / (1 - exp(-hc_k * df1.wav / Tref))
+            )  # [cm-1/(molecules/cm-2)]
+
+        else:
+            # An alternative strategy is to calculate the linestrength from the
+            # Einstein A coefficient and the populations (see Klarenaar 2017 Eqn. 12)
+
+            if not "gu" in df1:
+                if not "ju" in df1:
+                    self._add_ju(df1)
+                self._calc_degeneracies(df1)
+
+            Ia = self.get_lines_abundance(df1)
+            df1["S"] = linestrength_from_Einstein(
+                df1.A, df1.gu, df1.El, Ia, df1.wav, Qgas(Tgas), Tgas
+            )
 
         assert "S" in self.df1
 
@@ -2826,7 +2936,7 @@ def get_populations(self, levels="vib"):
 
         # To get isotopic abundance
         # placeholder # TODO: replace with attributes of Isotope>ElectronicState objects
-        molpar = MolParams()
+        molpar = self.molparam
 
         pops = {}
         # Loop over molecules, isotopes, electronic states
@@ -2880,11 +2990,11 @@ def get_populations(self, levels="vib"):
         # exported in a Spectrum but still connected to the Factory
         return pops
 
-    def _calc_linestrength_noneq(self):
-        """
+    def calc_linestrength_noneq(self):
+        """Calculate linestrengths at non-LTE
+
         Parameters
         ----------
-
         Pre-requisite:
 
             lower state population `nl` has already been calculated by
@@ -2893,7 +3003,6 @@ def _calc_linestrength_noneq(self):
 
         Returns
         -------
-
         None
             Linestrength `S` added in self.df
 
@@ -2912,7 +3021,9 @@ def _calc_linestrength_noneq(self):
         See Also
         --------
 
-        :meth:`~radis.lbl.base.BaseFactory.calc_populations_noneq`
+        :py:meth:`~radis.lbl.base.BaseFactory.calc_populations_noneq`,
+        :py:meth:`~radis.lbl.base.BaseFactory.calc_emission_integral`,
+        :py:func:`~radis.lbl.base.linestrength_from_Einstein`
 
         """
 
@@ -2934,63 +3045,62 @@ def _calc_linestrength_noneq(self):
 
         if "id" in df:
             id_set = df.id.unique()
-            id = id_set[0]
-
-        else:
-            id = df.attrs["id"]
-
-        iso_set = self._get_isotope_list()  # df1.iso.unique()
-
-        Qref_dict = {}
-
-        # TODO for multi-molecule code: add above line in the loop
-        if (("id" in df and len(id_set) == 1) or ("id" not in df and id)) and len(
-            iso_set
-        ) == 1:
-
-            # Shortcut if only 1 isotope. We attribute molar_mass & abundance
-            # as attributes of the line database, instead of columns. Much
-            # faster!
-
-            molecule = get_molecule(id)
-            state = self.input.state
-            parsum = self._get_parsum(molecule, iso_set[0], state)  # partition function
-            Qref_dict[iso_set[0]] = parsum.at(Tref, update_populations=False)
-
-        else:
+            if len(id_set) > 0:
+                raise NotImplementedError("> 1 molecules in same DataFrame")
+            else:
+                self.warn(
+                    "There shouldn't be a Column 'id' with a unique value",
+                    "PerformanceWarning",
+                )
+                df.attrs["id"] = int(id_set)
+        molecule = get_molecule(df.attrs["id"])
+        state = self.input.state
 
-            #            # normal method
-            #            # still much faster than the groupby().apply() method (see radis<=0.9.19)
-            #            # (tested + see https://stackoverflow.com/questions/44954514/efficient-way-to-conditionally-populate-elements-in-a-pandas-groupby-object-pos)
-            #
-            # partition function
-            if "id" in df:
-                dgb = df.groupby(by=["id", "iso"])
-                for (id, iso), idx in dgb.indices.items():
-                    molecule = get_molecule(id)
-                    state = self.input.state
-                    parsum = self._get_parsum(molecule, iso, state)
-                    Qref_dict[iso] = parsum.at(Tref, update_populations=False)
+        # Partition functions
+        def Qgas(Tref):
+            """Aggregate the values of Qgas
 
-                    if radis.DEBUG_MODE:
-                        assert (df.loc[idx, "id"] == id).all()
-                        assert (df.loc[idx, "iso"] == iso).all()
+            Parameters
+            ----------
+            Tref: float (K)
+                reference gas temperature of database
 
-            else:
-                dgb = df.groupby(by=["iso"])
-                id = df.attrs["id"]
-                molecule = get_molecule(id)
+            Returns
+            -------
+            float or dict: Returns Qgas as a dictionary with isotope values as its keys
 
-                for (iso), idx in dgb.indices.items():
-                    state = self.input.state
-                    parsum = self._get_parsum(molecule, iso, state)
-                    Qref_dict[iso] = parsum.at(Tref, update_populations=False)
+            """
+            self.profiler.start("corrected_population_se", 3)
+            if "iso" in df:
+                iso_set = df.iso.unique()
+                if len(iso_set) == 1:
+                    self.warn(
+                        "There shouldn't be a Column 'iso' with a unique value",
+                        "PerformanceWarning",
+                    )
+                    iso = int(iso_set)
+                    parsum = self.get_partition_function_interpolator(
+                        molecule, iso, state
+                    )
+                    Q = parsum.at(Tref, update_populations=False)
+                    df.attrs["Q"] = Q
+                else:
+                    Qref_dict = {}
+                    for iso in iso_set:
+                        parsum = self.get_partition_function_interpolator(
+                            molecule, iso, state
+                        )
+                        Qref_dict[iso] = parsum.at(Tref, update_populations=False)
+                    Q = df["iso"].map(Qref_dict)
 
-                    if radis.DEBUG_MODE:
-                        assert (df.loc[idx, "iso"] == iso).all()
+            else:  # "iso" not in df:
+                iso = df.attrs["iso"]
+                parsum = self.get_partition_function_interpolator(molecule, iso, state)
+                Q = parsum.at(Tref, update_populations=False)
+                df.attrs["Q"] = Q
+            self.profiler.stop("map_part_func", "map partition functions")
+            return Q
 
-        self.profiler.stop("map_part_func", "map partition functions")
-        self.profiler.start("corrected_population_se", 3)
         # Correct linestrength
 
         # ... populations without abundance dependance (already in linestrength)
@@ -2998,17 +3108,14 @@ def _calc_linestrength_noneq(self):
         nl = df.nl
         # ... remove Qref, nref, etc.
         # start from for tabulated linestrength
-        line_strength = df.int.copy()
+        line_strength = df.int.copy()  # TODO: savememory; replace the "int" column
+        if not self.molparam.terrestrial_abundances:
+            raise NotImplementedError(
+                "Formula not corrected for non-terrestrial isotopic abundances"
+            )
 
         # ... correct for lower state population
-        if "iso" in df:
-            line_strength /= (
-                df.gl * exp(-hc_k * df.El / Tref) / df["iso"].map(Qref_dict)
-            )
-        else:
-            assert len(Qref_dict) == 1
-            Qref = list(Qref_dict.values())[0]
-            line_strength /= df.gl * exp(-hc_k * df.El / Tref) / Qref
+        line_strength /= df.gl * exp(-hc_k * df.El / Tref) / Qgas(Tref)
         line_strength *= nl
 
         # ... correct effect of stimulated emission
@@ -3027,21 +3134,26 @@ def _calc_linestrength_noneq(self):
         return  # df1 automatically updated
 
     # %%
-    def _calc_emission_integral(self):
-        """Calculate Emission Integral.
+    def calc_emission_integral(self):
+        r"""Calculate Emission Integral.
+
+        .. math::
+            Ei=\frac{n_u A_{ul}}{4} \pi \Delta E_{ul}
 
         Emission Integral is a non usual quantity introduced here to have an
         equivalent of Linestrength in emission calculation
 
         Returns
         -------
-
         None
             Emission integral `Ei` added in self.df
 
+        Notes
+        -----
 
         emission_integral: (mW/sr)
-            emission integral is defined as
+            emission integral is defined as::
+
             Ei = n_u * A_ul / 4π * DeltaE_ul
                   :      :     :      :
                 (#/#)  (s-1) (sr)    (mJ)
@@ -3051,6 +3163,10 @@ def _calc_emission_integral(self):
                 ϵ(λ)   =      Ei  *   Phi(λ)  * ntot   * path_length
                  :             :         :       :          :
             mW/cm2/sr/nm    (mW/sr)   (1/nm)   (cm-3)      (cm)
+
+        See Also
+        --------
+        :py:meth:`~radis.lbl.base.BaseFactory.calc_linestrength_noneq`
         """
 
         df = self.df1
@@ -3070,7 +3186,7 @@ def _calc_emission_integral(self):
         # adim. (#/#) (multiplied by n_tot later)
         n_u = df["nu"]
         # correct for abundance
-        n_ua = n_u * self.get_abundance(df)
+        n_ua = n_u * self.get_lines_abundance(df)
 
         A_ul = df["Aul"]  # (s-1)
 
@@ -3564,7 +3680,7 @@ def get_waverange(
 
 
 def linestrength_from_Einstein(A, gu, El, Ia, nu, Q, T):
-    """Calculate linestrength at temperature ``T`` from Einstein coefficients.
+    r"""Calculate linestrength at temperature ``T`` from Einstein coefficients.
 
     Parameters
     ----------
@@ -3592,19 +3708,21 @@ def linestrength_from_Einstein(A, gu, El, Ia, nu, Q, T):
     ----------
 
     .. math::
-        S=\frac{A gu \operatorname{exp}\left(\frac{-hc_k El}{T}\right) \left(1-\operatorname{exp}\left(\frac{-hc_k \nu}{T}\right)\right)}{8\pi c_{CGS} {\nu}^2 Q}
+        S(T) =\frac{1}{8\pi c_{CGS} {n_u}^2} A \frac{I_a g_u \operatorname{exp}\left(\frac{-c_2 E_l}{T}\right)}{Q(T)} \left(1-\operatorname{exp}\left(\frac{-c_2 n_u}{T}\right)\right)
 
-    Combine Eq.(A.5), (A.8) in [Rothman-1998]_
+    Combine Eq.(A.5), (A.9) in [Rothman-1998]_
 
+    See Also
+    --------
+    :py:meth:`~radis.lbl.base.BaseFactory.calc_linestrength_eq`
 
     """
 
     return (
-        A
-        * gu
-        * np.exp(-hc_k * El / T)
+        (1 / (8 * np.pi * c_CGS * nu ** 2))
+        * A
+        * ((Ia * gu * np.exp(-hc_k * El / T)) / Q)
         * (1 - np.exp(-hc_k * nu / T))
-        / (8 * np.pi * c_CGS * nu ** 2 * Q)
     )
 
 
diff --git a/radis/lbl/broadening.py b/radis/lbl/broadening.py
index 5017e7eac..65546d9b0 100644
--- a/radis/lbl/broadening.py
+++ b/radis/lbl/broadening.py
@@ -69,7 +69,6 @@
 from numpy import pi, sin, sqrt, trapz, zeros_like
 from scipy.signal import oaconvolve
 
-from radis.db.molparam import MolParams
 from radis.lbl.base import BaseFactory
 from radis.misc.arrays import add_at, numpy_add_at
 from radis.misc.basics import is_float
@@ -956,7 +955,7 @@ def _add_voigt_broadening_HWHM(self, df, pressure_atm, mole_fraction, Tgas, Tref
         else:
             Tdpsel = df.Tdpsel
 
-        molar_mass = get_molar_mass(df)
+        molar_mass = self.get_molar_mass(df)
 
         # Calculate broadening FWHM
         wv, wl, wg = voigt_broadening_HWHM(
@@ -1053,7 +1052,7 @@ def _add_doppler_broadening_HWHM(self, df, Tgas):
             - ``hwhm_gauss``
         """
 
-        molar_mass = get_molar_mass(df)
+        molar_mass = self.get_molar_mass(df)
 
         # Calculate broadening HWHM
         wg = doppler_broadening_HWHM(df.wav, molar_mass, Tgas)
@@ -2850,48 +2849,6 @@ def rough_sum_on_grid():
     )
 
 
-def get_molar_mass(df):
-    """Returns molar mass.
-
-    Parameters
-    ----------
-    df: dataframe
-
-    Returns
-    -------
-    The molar mass of all the isotopes in the dataframe
-    """
-    molpar = MolParams()
-
-    if "id" in df.columns:
-        raise NotImplementedError(">1 molecule")
-    elif "id" in df.attrs:
-        id = df.attrs["id"]
-    else:
-        # HARDCODED molar mass; for WIP ExoMol implementation, until MolParams
-        # is an attribute and can be updated with definitions from ExoMol.
-        # https://github.com/radis/radis/issues/321
-        HARCODED_MOLAR_MASS = {"SiO": {1: 43.971842}}
-        try:
-            return HARCODED_MOLAR_MASS[df.attrs["molecule"]][df.attrs["iso"]]
-        except KeyError:
-            raise NotImplementedError
-
-    if "iso" in df.columns:
-        iso_set = df.iso.unique()
-        molar_mass_dict = {}
-        for iso in iso_set:
-            molar_mass_dict[iso] = molpar.get(id, iso, "mol_mass")
-        molar_mass = df["iso"].map(molar_mass_dict)
-    else:
-        iso = df.attrs["iso"]
-        molar_mass = molpar.get(id, iso, "mol_mass")
-
-    return molar_mass
-
-    #
-
-
 if __name__ == "__main__":
 
     from radis.test.lbl.test_broadening import _run_testcases
diff --git a/radis/lbl/factory.py b/radis/lbl/factory.py
index 9fdf8e5f5..ffc4d7ba0 100644
--- a/radis/lbl/factory.py
+++ b/radis/lbl/factory.py
@@ -88,7 +88,6 @@
 from radis import version
 from radis.db import MOLECULES_LIST_EQUILIBRIUM, MOLECULES_LIST_NONEQUILIBRIUM
 from radis.db.classes import get_molecule, get_molecule_identifier
-from radis.db.molparam import MolParams
 
 try:  # Proper import
     from .bands import BandFactory
@@ -946,7 +945,7 @@ def eq_spectrum_gpu(
 
         ### GET ISOTOPE ABUNDANCE & MOLECULAR MASS ###
 
-        molpar = MolParams()
+        molpar = self.molparam
 
         try:
             if "id" in self.df0:
@@ -1337,8 +1336,8 @@ def non_eq_spectrum(
                 overpopulation=overpopulation,
             )
 
-        self._calc_linestrength_noneq()
-        self._calc_emission_integral()
+        self.calc_linestrength_noneq()
+        self.calc_emission_integral()
 
         # ----------------------------------------------------------------------
         # Cutoff linestrength
@@ -1761,7 +1760,7 @@ def optically_thin_power(
         else:
             self.calc_populations_eq(Tgas)
             self.df1["Aul"] = self.df1.A  # update einstein coefficients
-        self._calc_emission_integral()
+        self.calc_emission_integral()
 
         #        # ----------------------------------------------------------------------
         #        # Cutoff linestrength  (note that cuting linestrength doesnt make this
diff --git a/radis/lbl/loader.py b/radis/lbl/loader.py
index 5844a4cb4..6d177a1a7 100644
--- a/radis/lbl/loader.py
+++ b/radis/lbl/loader.py
@@ -59,6 +59,7 @@
 
 from radis.db.classes import get_molecule
 from radis.db.molecules import getMolecule
+from radis.db.molparam import MolParams
 from radis.io.cache_files import cache_file_name
 from radis.io.cdsd import cdsd2df
 from radis.io.exomol import fetch_exomol
@@ -69,7 +70,7 @@
 from radis.io.tools import drop_object_format_columns, replace_PQR_with_m101
 from radis.levels.partfunc import (
     PartFunc_Dunham,
-    PartFuncHAPI,
+    PartFuncTIPS,
     RovibParFuncCalculator,
     RovibParFuncTabulator,
 )
@@ -657,6 +658,12 @@ def __init__(self):
 
         self._autoretrieveignoreconditions = []  # HACK. See _retrieve_from_database
 
+        # Molecular parameters
+        self.molparam = MolParams()
+        """MolParam: contains information about molar mass; isotopic abundance.
+
+        See :py:class:`~radis.db.molparam.MolParams`"""
+
     def _reset_profiler(self, verbose):
         """Reset :py:class:`~radis.misc.profiler.Profiler`
 
@@ -908,8 +915,10 @@ def fetch_databank(
         if exomol_database != None:
             assert source == "exomol"
         if [parfuncfmt, source].count("exomol") == 1:
-            raise NotImplementedError(
-                "ExoMol partition functions must be used with ExoMol database - and vice-versa"
+            self.warn(
+                f"Using lines from {source} but partition functions from {parfuncfmt}"
+                + "for consistency we recommend using lines and partition functions from the same database",
+                "AccuracyWarning",
             )
 
         # Get inputs
@@ -2179,7 +2188,7 @@ def _build_partition_function_interpolator(
             # Use TIPS-2017 through HAPI (HITRAN Python interface, integrated in RADIS)
             # no tabulated partition functions defined. Only non-eq spectra can
             # be calculated if energies are also given
-            parsum = PartFuncHAPI(
+            parsum = PartFuncTIPS(
                 M=molecule, I=isotope, path=parfunc, verbose=self.verbose
             )
         elif parfuncfmt == "cdsd":  # Use tabulated CDSD partition functions
@@ -2268,6 +2277,103 @@ def _build_partition_function_calculator(
 
         return parsum
 
+    def get_abundance(self, molecule, isotope):
+        """Get isotopic abundance
+
+        Parameters
+        ----------
+        molecule: str
+        isotope: int, or list
+            isotope number, sorted in terrestrial abundance
+
+        Examples
+        --------
+        Use it from SpectrumFactory::
+
+            sf.get_abundance("H2O", 1)
+            sf.get_abundance("CH4", [1,2,3])
+
+        .. minigallery:: radis.lbl.loader.DatabankLoader.get_abundance
+            :add-heading:
+
+        See Also
+        --------
+        :py:meth:`~radis.lbl.loader.DatabankLoader.set_abundance`
+        """
+
+        if isinstance(molecule, str):
+            from radis.db.classes import get_molecule_identifier
+
+            molecule = get_molecule_identifier(molecule)
+
+        if isinstance(isotope, int):
+            return self.molparam.df.loc[(molecule, isotope)].abundance
+        elif isinstance(isotope, list):
+            return np.array(
+                [self.molparam.df.loc[(molecule, iso)].abundance for iso in isotope]
+            )
+        else:
+            raise ValueError(isotope)
+
+    def set_abundance(self, molecule, isotope, abundance):
+        """Set isotopic abundance
+
+        Parameters
+        ----------
+        molecule: str
+        isotope: int, or list
+            isotope number, sorted in terrestrial abundance
+        abundance: float, or list
+
+        Examples
+        --------
+
+            from radis import SpectrumFactory
+
+            sf = SpectrumFactory(
+                2284.2,
+                2284.6,
+                wstep=0.001,  # cm-1
+                pressure=20 * 1e-3,  # bar
+                mole_fraction=400e-6,
+                molecule="CO2",
+                isotope="1,2",
+                verbose=False
+            )
+            sf.load_databank("HITEMP-CO2-TEST")
+            print("Abundance of CO2[1,2]", sf.get_abundance("CO2", [1, 2]))
+            sf.eq_spectrum(2000).plot("abscoeff")
+
+            #%% Set the abundance of CO2(626) to 0.8; and the abundance of CO2(636) to 0.2 (arbitrary):
+            sf.set_abundance("CO2", [1, 2], [0.8, 0.2])
+            print("New abundance of CO2[1,2]", sf.get_abundance("CO2", [1, 2]))
+
+            sf.eq_spectrum(2000).plot("abscoeff", nfig="same")
+
+        .. minigallery:: radis.lbl.loader.DatabankLoader.set_abundance
+            :add-heading
+
+        See Also
+        --------
+        :py:meth:`~radis.lbl.loader.DatabankLoader.get_abundance`
+
+        """
+
+        if isinstance(molecule, str):
+            from radis.db.classes import get_molecule_identifier
+
+            molecule = get_molecule_identifier(molecule)
+
+        self.molparam.terrestrial_abundances = False
+
+        if isinstance(isotope, int):
+            self.molparam.df.loc[(molecule, isotope), "abundance"] = abundance
+        elif isinstance(isotope, list):
+            assert len(isotope) == len(abundance)
+            self.molparam.df.loc[(molecule, isotope), "abundance"] = abundance
+        else:
+            raise ValueError(isotope)
+
     def get_partition_function_interpolator(self, molecule, isotope, elec_state):
         """Retrieve Partition Function Interpolator.
 
diff --git a/radis/levels/partfunc.py b/radis/levels/partfunc.py
index 5d4541354..447efbc59 100644
--- a/radis/levels/partfunc.py
+++ b/radis/levels/partfunc.py
@@ -10,7 +10,7 @@
 and require a list of energies
 
 Tabulators are more specific, and require a list of precalculated
-partition functions at different temperature. PartFuncHAPI
+partition functions at different temperature. PartFuncTIPS
 uses the HITRAN hapi.py [1]_ module to interpolate Q for HITRAN species
 
 
@@ -19,7 +19,7 @@
 
 Partition functions:
 
-- :class:`~radis.levels.partfunc.PartFuncHAPI`
+- :class:`~radis.levels.partfunc.PartFuncTIPS`
 - :class:`~radis.levels.partfunc.PartFunc_Dunham`
 
 Which inherit from:
@@ -143,8 +143,7 @@ def at(self, T, **kwargs):
 
         See Also
         --------
-
-        :py:class:`~radis.levels.partfunc.PartFuncHAPI`
+        :py:class:`~radis.levels.partfunc.PartFuncTIPS`
         """
 
         # For compatibility with the syntax of RovibParFuncCalculator.at()
@@ -1104,8 +1103,9 @@ def _at(self, T):
         return np.interp(T, self.T_range, self.Q_range)
 
 
-class PartFuncHAPI(RovibParFuncTabulator):
-    """Return partition function using interpolation of tabulated values.
+class PartFuncTIPS(RovibParFuncTabulator):
+    """Return partition function using interpolation of tabulated values
+    using the TIPS program [TIPS-2020]_
 
     Parameters
     ----------
@@ -1120,28 +1120,26 @@ class PartFuncHAPI(RovibParFuncTabulator):
     --------
     ::
 
-        from radis.levels.partfunc import PartFuncHAPI
+        from radis.levels.partfunc import PartFuncTIPS
         from radis.db.classes import get_molecule_identifier
 
         M = get_molecule_identifier('N2O')
         iso=1
 
-        Q = PartFuncHAPI(M, iso)
+        Q = PartFuncTIPS(M, iso)
         print(Q.at(T=1500))
 
     See :ref:`online examples <label_examples_partition_functions>` for more.
 
     References
     ----------
-    Partition function are calculated with HAPI [1]_ (Hitran Python Interface) using
+    Partition function are retrieved from [TIPS-2020]_ through [HAPI]_  (Hitran Python Interface) using
     partitionSum(M,I,T)
-
-    .. [1] `HAPI: The HITRAN Application Programming Interface <http://hitran.org/hapi>`_
     """
 
     def __init__(self, M, I, path=None, verbose=True):
 
-        super(PartFuncHAPI, self).__init__()
+        super(PartFuncTIPS, self).__init__()
 
         self.verbose = verbose
 
@@ -1166,7 +1164,7 @@ def __init__(self, M, I, path=None, verbose=True):
             M = get_molecule_identifier(M)
         if type(M) is not int:
             raise TypeError("Molecule id must be int: got {0} ({1})".format(M, type(M)))
-        if type(I) is not int:
+        if type(I) not in [int, np.int64]:
             raise TypeError(
                 "Isotope number must be int: got {0} ({1})".format(I, type(I))
             )
@@ -1199,7 +1197,7 @@ def _at(self, T):
         try:
             return self.partitionSum(self.M, self.I, T)
         except Exception as err:
-            if "TIPS2017" in str(err):
+            if "TIPS2017" in str(err) or "TIPS2020" in str(err):
                 # TIPS 2017 OutOfBoundError
                 from radis.db import MOLECULES_LIST_NONEQUILIBRIUM
 
@@ -1235,6 +1233,9 @@ def _get_cachefile_name(ElecState):
     return filename
 
 
+PartFuncHAPI = PartFuncTIPS  # old name; for compatibility
+
+
 class PartFunc_Dunham(RovibParFuncCalculator):
     """Calculate partition functions from spectroscopic constants, if
     molecule data is available in RADIS. If using the
diff --git a/radis/spectrum/spectrum.py b/radis/spectrum/spectrum.py
index 1a0d8c61b..622c3c038 100644
--- a/radis/spectrum/spectrum.py
+++ b/radis/spectrum/spectrum.py
@@ -3444,7 +3444,7 @@ def check_wavespace(w):
 
     def __eq__(self, other):
         """Override the default Equals behavior."""
-        return self.compare_with(other, verbose=False, plot=False)
+        return self.compare_with(other, spectra_only=True, verbose=False, plot=False)
 
     def __ne__(self, other):
         """Define a non-equality test."""


Test patch:
diff --git a/radis/test/lbl/test_base.py b/radis/test/lbl/test_base.py
index d24f04b97..fd8c4b6bf 100644
--- a/radis/test/lbl/test_base.py
+++ b/radis/test/lbl/test_base.py
@@ -20,6 +20,70 @@
 from radis.test.utils import setup_test_line_databases
 
 
+def test_linestrength_calculations(*args, **kwargs):
+    """Compare and validate linestrength calculations from
+    :py:func:`radis.lbl.base.linestrength_from_Einstein`,
+    :py:meth:`radis.lbl.base.BaseFactory.calc_linestrength_eq`,
+    :py:meth:`radis.lbl.base.BaseFactory.calc_linestrength_non_eq`,
+
+    """
+
+    # Test linestrength calculations
+    from radis.io.hitran import hit2df
+    from radis.test.utils import getTestFile
+
+    df = hit2df(getTestFile("hitran_co_3iso_2000_2300cm.par"))
+
+    from radis.db.classes import get_molecule_identifier
+    from radis.db.molparam import MolParams
+    from radis.lbl.base import linestrength_from_Einstein
+    from radis.levels.partfunc import PartFuncTIPS
+
+    molpar = MolParams()
+    M = get_molecule_identifier("CO")
+
+    Q_arr = df.iso.map({iso: PartFuncTIPS(M, iso).at(296) for iso in df.iso.unique()})
+    Ia_arr = df.iso.map(
+        {iso: molpar.get(M, iso, "abundance") for iso in df.iso.unique()}
+    )
+
+    # Reference linestrengths recomputed from Einstein coefficients match
+    # within 0.1%
+    S = linestrength_from_Einstein(df.A, df.gp, df.El, Ia_arr, df.wav, Q_arr, 296)
+    assert np.allclose(df.int, S, rtol=1e-3, atol=0)
+
+    #%% Now, at 1000 K
+    Q_arr = df.iso.map({iso: PartFuncTIPS(M, iso).at(1000) for iso in df.iso.unique()})
+    S_1000 = linestrength_from_Einstein(df.A, df.gp, df.El, Ia_arr, df.wav, Q_arr, 1000)
+
+    from radis.lbl.factory import SpectrumFactory
+
+    sf = SpectrumFactory(2000, 3000)
+    sf.load_databank(
+        path=getTestFile("hitran_co_3iso_2000_2300cm.par"),
+        format="hitran",
+        parfuncfmt="hapi",
+    )
+
+    # TODO : write an example of all the calculation steps in SpectrumFactory
+    # sf._reinitialize()  # creates scaled dataframe df1 from df0  # TODO: make a public function.
+    # sf.calc_linestrength_eq(1000)
+
+    sf.eq_spectrum(1000)
+    # Linestrengths computed from Einstein coefficients by linestrength_from_Einstein
+    # match the one from the SpectrumFactory within 0.1%
+    assert np.allclose(sf.df1.S, S_1000, rtol=1e-3, atol=0)
+
+    #%% Using non-eq calculations
+    S_eq = sf.df1.S.copy()
+
+    sf.non_eq_spectrum(1000, 1000)
+    # Linestrengths computed with equilibrium routine (scaling reference) match
+    # the ones from the nonequilibrium routine (from Einstein coefficients)
+    # within 0.1%
+    assert np.allclose(sf.df1.S, S_eq, rtol=1e-3, atol=0)
+
+
 @pytest.mark.fast
 def test_export_populations(plot=True, verbose=True, warnings=True, *args, **kwargs):
     """Check populations calculated in the nonequilibrium module are exported in Spectrum.
@@ -655,6 +719,7 @@ def test_get_waverange(*args, **kwargs):
 
 def _run_testcases(verbose=True, plot=True):
 
+    test_linestrength_calculations()
     test_export_populations(plot=plot, verbose=verbose)
     test_export_rovib_fractions(plot=plot, verbose=verbose)
     test_populations_CO2_hamiltonian(plot=plot, verbose=verbose)
diff --git a/radis/test/lbl/test_loader.py b/radis/test/lbl/test_loader.py
index 56b1fdde2..beef83113 100644
--- a/radis/test/lbl/test_loader.py
+++ b/radis/test/lbl/test_loader.py
@@ -50,7 +50,6 @@ def test_retrieve_from_database(
             mole_fraction=400e-6,
             molecule="CO2",
             isotope=[1],
-            db_use_cached=True,
             medium="vacuum",
             broadening_max_width=10,
             export_populations="rovib",
@@ -154,11 +153,80 @@ def test_ignore_irrelevant_files(*args, **kwargs):
         sf3.load_databank(path=test_file, format="cdsd-hitemp", parfuncfmt="hapi")
 
 
+@pytest.mark.fast
+def test_custom_abundance(verbose=True, plot=False, *args, **kwargs):
+    """Test calculatinos with custom abundances
+
+    set the terrestrial abundance to False to
+    force computation of linestrength from Einstein coefficients instead
+    of rescaling parameters
+
+    :py:meth:`~radis.lbl.base.DatabankLoader.get_abundance`,
+    :py:meth:`~radis.lbl.base.DatabankLoader.set_abundance`
+    """
+
+    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",
+        broadening_max_width=10,
+        verbose=0,
+    )
+    sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
+    sf.load_databank("HITEMP-CO2-TEST")
+
+    #%% Compute a spectrum ; set the terrestrial abundance to False to
+    # force computation of linestrength from Einstein coefficients instead
+    # of rescaling parameters
+    s = sf.eq_spectrum(2000)
+    sf.molparam.terrestrial_abundances = False
+    s2 = sf.eq_spectrum(2000)
+    assert (
+        s != s2
+    )  # compare spectra; we should see small differences due to different methods
+    assert s.compare_with(s2, "abscoeff", rtol=0.5e-2, plot=False)
+    # assert np.allclose(2 * s2.get("abscoeff")[1], s.get("abscoeff")[1])
+
+    if plot:
+        from radis import plot_diff
+
+        plot_diff(s, s2)
+
+    #%% Now restart by changing the abundance itself ; as a user would do
+
+    # ... check abundance get/set functions work well
+    abundance12_0 = sf.get_abundance("CO2", [1, 2])
+
+    if verbose:
+        print("Abundance: CO2[1]", sf.get_abundance("CO2", 1))
+        print("Abundance: CO2[1,2]", sf.get_abundance("CO2", [1, 2]))
+
+    sf.set_abundance("CO2", [1, 2], abundance12_0 / 3)
+    assert all(sf.get_abundance("CO2", [1, 2]) == abundance12_0 / 3)
+    sf.set_abundance("CO2", 1, abundance12_0[0] / 3)
+    sf.set_abundance("CO2", 2, abundance12_0[1] / 3)
+    assert all(sf.get_abundance("CO2", [1, 2]) == abundance12_0 / 3)
+
+    # Below we compute a spectrum with the abundance divided by 3; and then
+    # check that the resulting absorption coefficient is 3x smaller
+
+    s3 = sf.eq_spectrum(2000)
+    assert s.compare_with((3 * s3.take("abscoeff")), "abscoeff", rtol=0.5e-2, plot=True)
+
+
 def _run_testcases(verbose=True, plot=False):
 
     test_ignore_cached_files()
     test_retrieve_from_database(plot=plot, verbose=verbose)
     test_ignore_irrelevant_files(verbose=verbose)
+    test_custom_abundance()
 
 
 if __name__ == "__main__":
diff --git a/radis/test/levels/test_partfunc.py b/radis/test/levels/test_partfunc.py
index d8e71a0e4..c3ec8e14d 100644
--- a/radis/test/levels/test_partfunc.py
+++ b/radis/test/levels/test_partfunc.py
@@ -29,7 +29,7 @@
 
 from radis import SpectrumFactory
 from radis.db.molecules import Molecules
-from radis.levels.partfunc import PartFunc_Dunham, PartFuncHAPI
+from radis.levels.partfunc import PartFunc_Dunham, PartFuncTIPS
 from radis.levels.partfunc_cdsd import PartFuncCO2_CDSDcalc, PartFuncCO2_CDSDtab
 from radis.misc.printer import printm
 from radis.misc.warning import DeprecatedFileWarning
@@ -222,7 +222,7 @@ def test_calculatedQ_match_HAPI_CO(
 
     #    if plot: db.plot_states()
 
-    hapi = PartFuncHAPI(M=5, I=1)  # CO  # isotope
+    hapi = PartFuncTIPS(M=5, I=1)  # CO  # isotope
 
     us = []
     hap = []
@@ -288,7 +288,7 @@ def test_calculatedQ_match_HAPI(plot=False, verbose=True, *args, **kwargs):
 
         from radis.db.classes import get_molecule_identifier
 
-        hapi = PartFuncHAPI(M=get_molecule_identifier(molecule), I=iso)
+        hapi = PartFuncTIPS(M=get_molecule_identifier(molecule), I=iso)
 
         Q_radis = db.at(T)
         Q_hapi = hapi.at(T)