pauliacomi / pyGAPS

A framework for processing adsorption data and isotherm fitting
MIT License
63 stars 24 forks source link

Issue in Dubinin calculations #43

Closed sblanky closed 1 year ago

sblanky commented 1 year ago

Describe the bug I think Dubinin volume calculations are incorrect; they are overestimated due to the order of the calculations.

Current output

To Reproduce See dr_da_plots:

line 223

    pressure, loading = get_iso_loading_and_pressure_ordered(
        isotherm, branch, {
            "loading_basis": "molar",
            "loading_unit": "mol"
        }, {"pressure_mode": "relative"}
    )

and line 403

def log_v_adj(loading, molar_mass, liquid_density):
    """Base 10 logarithm of volumetric uptake."""
    return numpy.log10(loading * molar_mass / liquid_density)

I think that this is trying to combine the calculation of limiting micropore capacity and limiting micropore volume with the log(V) variable. This results in incorrect (negative) values for log(V).

System (please complete the following information):

Additional context I am working on a patch for my own use which first takes the loading in cm3

    pressure, loading = get_iso_loading_and_pressure_ordered(
        isotherm, branch, {
            "loading_basis": "volume_gas",
            "loading_unit": "cm3"
        }, {"pressure_mode": "relative"}
    )

Then takes just takes the logarithm of the loading;

def log_v_adj(loading, molar_mass, liquid_density):
    """Base 10 logarithm of volumetric uptake."""
    return numpy.log10(loading)

This means that the limiting micropore capacity is calculated instead of micropore volume

This is converted to micropore volume using the density conversion factor 0.0015500 (for N2 at 77 K) that micromeritics uses. I think this is just the conversion factor from gas to liquid density. The results are much more reasonable than above (and they match MicroActive results).

my results

I may have just totally missed something somewhere and be using the code wrong though!

sblanky commented 1 year ago

I have included a rough fix that works for my data (see above commit on my fork). the loadings read in are already in cm3/g stp, so I don't know if this will work on other units/bases. If you look over this and think it makes sense, then let me know and I'll formally submit a pull request etc.

If not, and I just don't understand how this code is meant to work then ignore me!

pauliacomi commented 1 year ago

Interesting difference. A couple of questions to clarify

If the above are correct, the changes made should not have any impact on the final number, so there may be something else at play...

pauliacomi commented 1 year ago

PS: I don't remember why I did the calcs in base 10, it just makes it more complicated.

sblanky commented 1 year ago

Your definitions are correct. And yes, this magic number is just the ratio between gas density/liquid density. Unfortunately, CoolProp doesn't have density of nitrogen at 298 K. On running;

adsorbate.gas_density(298)

The following error is produced;

Thermodynamic backend failed with error Temperature to QT_flash [298 K] must be in range [63.051 K, 126.192 K]. Attempting to read parameters dictionary...

The aif file that the isotherm is imported from has relative pressure units and loading units of 'cm3/g stp' - don't know if that reveals anything?

pauliacomi commented 1 year ago

CoolProp won't work on supercritical gases, true. If it's STP, density should be in 22.4 cm3/mmol / 28k mmol/g (N2 molar mass).

So definitions are right. Note that the volume_gas loading option will give you the isotherm volume of gas at isotherm temperature, so at 77K. In the changes

microp_volume = microp_capacity * gas_density / liquid_density

In terms of units and concepts

[cm3/g, gas at 77k] * [cm3/g, gas at 303k] / [cm3/g, liquid at 77K]

So it's not apples to apples.

sblanky commented 1 year ago

Grand, so I tested this on a file, and it turns out that I'd just been creating the point isotherm, and then the resultant aif file incorrectly.

Essentially;

file = pd.read_csv('isotherm.csv', names=['pressure', 'loading], header=0)
 isotherm = pg.PointIsotherm(
        pressure=file['pressure'],
        loading=file['loading'],
        material='example',
        adsorbate='N2',
        temperature=77,
        temperature_unit='K',
        pressure_unit='bar',
        pressure_mode='relative',
        loading_basis='volume_gas',
        loading_unit='cm3',
        material_basis='mass',
        material_unit='g',
    )
isotherm.to_aif('isotherm.aif')

So the preamble in the resulting aif file looks like;

...
_pygaps_pressure_mode 'relative'
_pygaps_pressure_unit 'None'
_pygaps_material_basis 'mass'
_pygaps_material_unit 'g'
_pygaps_loading_basis 'volume_gas'
_pygaps_loading_unit 'cm3'
_pygaps_temperature_unit 'K'

Changing to

...
_pygaps_loading_basis 'molar'
_pygaps_loading_unit 'cm3(STP)'
...

Fixes this, and running verbose gives much more reasonable pore volume;

isotherm = pgp.isotherm_from_aif('isotherm.aif')
pgc.dr_da_plots.da_plot(isotherm, p_limits=[0,0.1], verbose=True)

Returns

Exponent is: 3
Micropore volume is: 0.773 cm3/g
Effective adsorption potential is 6.35 kJ/mol

I wonder if it's worth warning when the user attempts to create an isotherm with the specific units and bases I have shown here? To me it was the most intuitive way to do so.

Anywho, closing this issue!

pauliacomi commented 1 year ago

Good stuff! If it was up to me, I would remove cm3[stp] completely...