BUG: strange results from elastic tensor calculation

[computron]

Describe the bug
Run the elastic tensor calculation for GaAs. Most things look OK but the Debye temperature and K_Reuss look very off, and clark_thermalcond is null.

To Reproduce
Steps to reproduce the behavior:
Run the following script:
ff_elastic.py

See screenshot for strange behavior highlighted. Expect these values to be reasonable.

Screenshots

[computron]

Actually, given the units, pretty much all the results in the previous screenshot are bad!

Important : results are catastrophically bad for primitive cells (Si and GaAs), look reasonable for conventional cells. The clark_thermalcond is still null in all cases but conventional cell at least looks reasonable for other parameters.

So questions are:

• why doesn't the workflow automatically convert to conventional cell?
• why are results so bad for primitive cells? For Si primitive cell I get a huge negative values for all derived elastic values - even worse than the GaAs example shown above. I would think the code could handle a rhombohedral cell if it is also calculating elastic tensors for monoclinic, triclinic, etc systems.

[computron]

Here is a clear set of scripts that show the difference between the conventional and primitive cell for Si. This difference exists regardless if you use a MongoStore or FileStore, so this should not be confused with some one-off issue on FileStores (few things actually are).

ff_elastic_compare.py
db_compare.py

Run the first script to generate results for Si using the conventional cell then a conversion of the same structure to primitive cell. Use the second script to query the results which shows (1) MongoStore and (2) Vastly different k_vrh (repeat for any property you desire) for the same structure.

<maggma.stores.mongolike.MongoStore object at 0x10cb34690>

----
Number of matching documents
----
2

----
Projected documents
----

Result 1
{'_id': ObjectId('69c55ce16415f5838a8e3c9e'), 'uuid': '039e82ad-8615-4f4a-b5d4-8ba1d1c1a482', 'index': 1, 'output': {'derived_properties': {'k_vrh': 79.98984931747492}}}

Result 2
{'_id': ObjectId('69c55ce96415f5838a8e3ca9'), 'uuid': '84470765-66a0-4fd5-9f07-0bc5a7d49178', 'index': 1, 'output': {'derived_properties': {'k_vrh': -27213248.98944672}}}

----
k_vrh values
----
Result 1: 79.98984931747492
Result 2: -27213248.98944672

Again, these values are for the same structure in different settings!

Hypothesizing - perhaps the elastic workflow implicitly assumes a conventional cell. If so, this needs to be enforced.

[utf]

Hi Anubhav, thanks for this. I'm not sure exactly why this is happening yet. However, just to note that if you change the force field you get vastly different results. From my testing

MACE-MP-0
  Conventional K_VRH: 74.70 GPa
  Primitive K_VRH: 74.70 GPa

MACE-MPA-0
  Conventional K_VRH: 8682963.74 GPa
  Primitive K_VRH: 5380050.21 GPa

MATPES_PBE
  Conventional K_VRH: 79.99 GPa
  Primitive K_VRH: -30242775.50 GPa

MATPES_R2SCAN
  Conventional K_VRH: -26835136.48 GPa
  Primitive K_VRH: 80.51 GPa

EDIT: The above test was run using sga.get_primitive_standard_structure() rather than si_primitive = sga.find_primitive(). I've just noticed, if I revert back to the latter, I get:

MACE-MP-0
  Conventional K_VRH: 74.70 GPa
  Primitive K_VRH: 3845433.16 GPa

[computron]

Interesting, thanks for looking into this and reproducing.

At the very least, this looks to be a very big practical concern when using MLIPs to run this workflow - even advanced ones. It's not even the case that the conventional cell always works out OK, for MATPES_R2SCAN it seems the primitive worked better.

Do you know if the same issue happens for DFT / VASP versions of this workflow? I don't have it connected otherwise I would try myself.

[utf]

There is a test for elastic constants with VASP outputs, but currently this only for the conventional cell.

I think this may be due to choice of elastic tensor fitting method and noisyness of the MLIPs. Changing the fitting method from the default of "finite_difference" to "pseudoinverse" fixes the problem for all models, e.g.:

em = ElasticMaker.from_force_field_name(
    "MATPES_R2SCAN",
    fit_elastic_tensor_kwargs={"fitting_method": "pseudoinverse"}
)
conventional_job = em.make(si_conventional)
primitive_job = em.make(si_primitive)

MACE-MP-0
  Conventional K_VRH: 74.70 GPa
  Primitive K_VRH: 74.70 GPa

MACE-MPA-0
  Conventional K_VRH: 86.40 GPa
  Primitive K_VRH: 86.41 GPa

MATPES_PBE
  Conventional K_VRH: 79.99 GPa
  Primitive K_VRH: 79.99 GPa

MATPES_R2SCAN
  Conventional K_VRH: 79.99 GPa
  Primitive K_VRH: 79.99 GPa

Secondly, I noticed that the finite_difference fitting method has a parameter called tol with a default of 1e-10 (currently it is not possible to override this setting in atomate2 but adding support for this will be trivial). Changing this default to 1e-7 (randomly picked) also stabilises the fitting for all methods.

A few things to note:

• I tried increasing the strain magnitude but this did not fix the fitting errors.
• The default elastic tensor fitting method can be controlled through the atomate settings.
• The setting has always been "finite_difference" and it was this in atomate1
• "finite_difference" is required if fitting a 3rd order tensor. The default is 2nd order. However, if pseudoinverse is requested and the order is 3, it will override and use finite_difference automatically.
• MatCalc fits the elastic tensor through least squares - maybe this will be more reliable than the current finite difference.

In terms of fixing this, options include:

• Checking whether the implementation from matcalc works for 3rd order elastic tensors and can be implemented in pymatgen directly (we always defer to pymatgen for the fitting).
• Changing the default fitting method to "pseudoinverse". However, for psuedoinverse the equilibrium stress is not used. I don't know what the knock-on effects of this will be.
• Update the value of tol in the finite difference elastic tensor fitting. Unclear what else this could impact.

[utf]

One final note - this can also be fixed by changing the default dtype of the model to float64 rather than float32. E.g.

em = ElasticMaker.from_force_field_name(
    "MATPES_R2SCAN",
    calculator_kwargs": {"default_dtype": "float64"}
)

MACE-MP-0
  Conventional K_VRH: 74.70 GPa
  Primitive K_VRH: 74.70 GPa

MACE-MPA-0
  Conventional K_VRH: 86.22 GPa
  Primitive K_VRH: 86.22 GPa

Unfortunately, it doesn't seem like you can change the MATPES models dtype easily (or I wasn't able to figure it out), so I don't know if this is useful going forward.

[computron]

Honestly find it very confusing the range of things that cause the issue to be fixed. All seem to have some relation to numerical tolerances however, I have not really wrapped my head around why all these things work and why we are so crazily sensitive to things like float64 vs float32 or the tolerance value used used in the fit procedure. This makes it difficult for me to weigh in on the appropriate solution.

Any thoughts from others?

[JaGeo]

I would assume that the dtype also affects other numerical results, potentially drastically. Enforcing float64 as a default therefore would be a safe option, also for other applications. Nevetheless, one should probably test this solution on a few hundred compounds (prim. vs conv.)

For the other solutions: one would probably need to test this on DFT results and on test compounds. I would be hesitant to change defaults here without tests

[JaGeo]

I am wondering about something else: does it make a difference if this tag here is switched on or off?

atomate2/src/atomate2/forcefields/jobs.py

Line 115 in 106e7e6

fix_symmetry: bool = False

[utf]

Great suggestion. Yes, this fixes all models. Using the default float32 and:

em = ElasticMaker.from_force_field_name("MACE-MPA-0", mlff_kwargs={"fix_symmetry": True})

I get:

MACE-MP-0
  Conventional K_VRH: 74.70 GPa
  Primitive K_VRH: 74.70 GPa

MACE-MPA-0
  Conventional K_VRH: 86.21 GPa
  Primitive K_VRH: 86.22 GPa

MATPES_PBE
  Conventional K_VRH: 79.99 GPa
  Primitive K_VRH: 79.99 GPa

MATPES_R2SCAN
  Conventional K_VRH: 80.50 GPa
  Primitive K_VRH: 80.50 GPa

It seems this is due to subtle symmetry breaking due to noise. Do we already apply this fix for the mlff phonon makers? If not perhaps we should as I imagine they will be similarly impacted.

[JaGeo]

We do not but we likely should. It is even worse for models that a not equivariant

[JaGeo]

My personal view: we should enforce float64 where possible and make the symmetry constraints the default for optimizations before complex workflows. We should discuss this in the doc strings and documentations.

@esoteric-ephemera @utf @computron what do you think with regard to defaults?

[computron]

Edit: I misread the suggestion above, leaving this comment for record-keeping but it is outdated

My philosophy is that the workflow itself should be as robust as possible. Not only does this prevent future user confusion ("I changed my model from float64 to float32, expected a ~5% error, but now my results are ridiculous negative values") but it usually also has the knock-on effect of being more robust against problems for the "strange 2%" of edge case materials that are difficult even when the user follows the instructions.

To that end:

1. Switching default float64 is probably an overall good idea, for the short term fix and also because this is likely going to just yield more stable results
2. I would not recommend auto-switching to conventional cells. There has been no clear benefit shown, this might add significant compute cost (e.g., imagine running a million elastic tensors for simple 2-atom cells, only for them to be converted to 8 atom cells), and may create the false impression that the workflow requires conventional cells (it doesn't appear to).
3. Separately and in the background, someone should launch an investigation into the fitting methods to see if it is safe to relax some of these tolerances. If the problem were something like SYMPREC being set to 1e-10 we would all agree that just loosening that tolerance would be the best solution, it may also be the case for these fitting tolerances. Ideally, I would like to see a situation where if someone did use a model with float32 precision, that the results would be in-line with expectation (a little different than float64, but not nuts).

If this sounds OK let me know, I can coordinate investigating (3) on our end also.

[computron]

I should also mention that even the historical elastic workflow had a lot of problems with strange values, negative values, etc. We had to filter out much data that was inconsistent or simply wrong in terms of elastic tensor. So further gains in workflow robustness I think are possible, I don't think the old way was necessarily perfect or even fully optimized.

[esoteric-ephemera]

It might be hard to find a set of parameters that is stable for all MLIPs, given differences in architectures, over-fitting to equilibrium energetics to favor high scores on matbench

Using things like fix_symmetry = True comes with the caveat that the symmetry is only determined up to the tolerance defined by symprec. There are many edge cases where it's hard to reliably determine the symmetry

This example is a "nice" case where the symmetry is well defined

Unfortunately, it doesn't seem like you can change the MATPES models dtype easily (or I wasn't able to figure it out), so I don't know if this is useful going forward.

matgl unfortunately defines a global variable to set the dtype, which can also affect other models

I can think of a workaround for this. I'll open a PR hopefully today to start addressing these concerns