pymatgen icon indicating copy to clipboard operation
pymatgen copied to clipboard

Unit for Young's modulus

Open anyangml opened this issue 7 months ago • 10 comments

Python version

3.11.5

Pymatgen version

2025.5.28

Operating system version

No response

Current behavior

I think the unit for Young's modulus is wrong.

If I understand correctly, in the following calculation, the units for shear modulus and bulk modulus are assumed to be GPa. https://github.com/materialsproject/pymatgen/blob/20afc862d34f3b65e30ff50391c9dedc82252c93/src/pymatgen/analysis/elasticity/elastic.py#L199-L205

However the units are in eV/A^3

https://github.com/materialsproject/pymatgen/blob/20afc862d34f3b65e30ff50391c9dedc82252c93/src/pymatgen/analysis/elasticity/elastic.py#L188-L198

Expected Behavior

a unit conversion from eV/A^3 to GPa shall be applied before calculating Young's modulus

Minimal example


Relevant files to reproduce this bug

No response

anyangml avatar Jun 12 '25 03:06 anyangml

Thanks for reporting this. I haven't had time to check this but wanted to forward it to @tschaume and @esoteric-ephemera . Maybe this is also relevant to MP data.

JaGeo avatar Jul 03 '25 09:07 JaGeo

Good catch - the units in the docstr overall look inconsistent. I'll try to take a look as well

esoteric-ephemera avatar Jul 03 '25 16:07 esoteric-ephemera

I also came across this and opened a duplicate #4532.

As seen above the core seems to be that y_mod (and other properties like long_v) are assuming the tensor input is GPa, even though it can be anything, and also the doc strings ask for/state eV/Ang^3.

The test case for this code is a nice example since it inputs what appears to be a GPa tensor, and so the units of y_mod are then correct.

I think to resolve either:

  • y_mod should return without unit conversion (to be in "ev/Ang^3") to be consistent with the other properties of ElasticTensor.

But this will impact downstream properties perhaps. E.g. synder_ac uses long_v which I think has the same issue

def long_v(...)
    """
    ...
    Returns:
        float: longitudinal sound velocity (in SI units)
    """
    ...
    return (1e9 * (self.k_vrh + 4 / 3 * self.g_vrh) / mass_density) ** 0.5
  • Or y_mod should convert correctly to SI units as indicated by the docstring. But this then depends on the units used to construct ElasticTensor (which could be anything), which begs the question of handling units more carefully in ElasticTensor. But it could just be about changing the docstring to GPa (since that is the test case).

Also I think ASE uses eV/Ang^3 by default, is that where this docstring convention has come from?

harveydevereux avatar Oct 31 '25 08:10 harveydevereux

I came across this today. Let's say k_vrh and g_vrh are in eV/A^3 This implies, without the 9.0e9 conversion factor, y_mod is also in eV/A^3 SI unit is Pa (Pascal or J/m^3). The conversion is 1 eV/A^3 = 1.6e11 Pa and not 9.0e9!

I think this should be fixed ASAP.

This also becomes a problem when ElasticityCalc is used from matcalc, which has pymatgen as its dependency to calculate Young's modulus from the ElasticTensor.

Update: Why does the commit 47dde4f1828f68da410138005c5523167e398d66 in PR #4436 contain: return 9.0e9 * self.k_vrh * self.g_vrh * self.eV_A3_to_GPa / (3 * self.k_vrh + self.g_vrh) instead of return 1.0e9 * self.k_vrh * self.g_vrh * self.eV_A3_to_GPa / (3 * self.k_vrh + self.g_vrh)?

sreeharshab avatar Nov 04 '25 00:11 sreeharshab

@sreeharshab The factor of 9 comes from the definition of the Young's modulus and not a unit, see Eq. 21 of this paper: $E_\mathrm{young} = 9KG/(3K + G)$

The presumption of eV/Å3 units comes probably from VASP calculations where this would be the derived unit without any conversion

Overall, I'm not in favor of guessing what units have actually been input by the user, and just returning consistent units. There are tools in pymatgen, ase, and scipy to permit unit conversion. Thoughts @JaGeo or @mkhorton ?

esoteric-ephemera avatar Nov 07 '25 00:11 esoteric-ephemera

@esoteric-ephemera Thank you very much for the reply and for sharing the paper. I agree with you about returning consistent units (everything in eV/Å^3 to not break backward compatibility for other properties ).

sreeharshab avatar Nov 07 '25 02:11 sreeharshab

@esoteric-ephemera i think there are many such assumptions in pymatgen. In the phonon part, we always assume THz as the frequency unit.

I think it would be worth to think about best practices for unit validation. In any case, the documentation should say what unit should go in and what comes out.

JaGeo avatar Nov 07 '25 06:11 JaGeo

I agree on all the comments. Let's document the units properly as a first step. @esoteric-ephemera can you review the PR to fix the doc on the units?

@JaGeo If you encounter any such missing documentation, you can perhaps help submit a quick PR to fix them? Thanks.

As for best practices, all solutions have some issues, especially given that units are not even consistent between different DFT codes. One way is of course ot use something like pint to handle this, but there is a simplicity to dealing with python natives like float and numpy arrays that we will lose if we go with pint.

shyuep avatar Nov 07 '25 16:11 shyuep

Sure I'll look at #4436 today. Doing unit conversion with pint + the float with unit class in pymatgen might be more of a long term solution

esoteric-ephemera avatar Nov 07 '25 16:11 esoteric-ephemera

@shyuep of course! I was more referring to the fact that we assume such units without checking rather than to cases where the units are wrongly documented.

However, we had fixed a previous issue where the unit was wrongly documented in the doc strings in the phonon part recently.

JaGeo avatar Nov 07 '25 17:11 JaGeo