Exploding crystals with MACELES and MACE-OFF models

[husakm67-blip]

Hi

I work on a large test of the MACE-OFF and MACELES models on all suitable crystal structures from Cambridge Structure Database (approximately 250 000 - non disordered, no errors, no voids, only supported elements, volume less then 4000 A3 ).
The target is to check the FF can reproduce the experimentally determined crystal structure.
It works on most structured giving similar results as DFT and close to experiment...
Unfortunately some structures witch should be supported explode during optimization.
Used code is below.
Part of the BFGS geometry optimization run is below. The forces are nonsense, final geometry have nothing to do to reality.

If I will give the CSD code of such structures (e.g. ADAVAG, AFOREU,AGIQIS,AGUFOB) can anybody investigate whatever it is a general problem with the MACE technology/models ?

Used code:

from mace.calculators import MACECalculator
from ase.io import read
from ase.optimize import BFGS
from ase.filters import FrechetCellFilter
import time
import sys

#handle duplicated lib load (Win 11)
import os
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"

#instalation dependent parameters
models_path ='M:\Calculation_data\MACE\models\\'

#parameters
model="MACE-OFF24_medium.model"
precision = 'float64'
device =  'cuda'
fmax = 0.01
maxsteps = 1000

model_with_path = models_path+model

#print setup
print('Model: ', model_with_path)
print('Precision: ',precision)
print('Device: ',device)
print('Fmax: ',fmax)
print('Maximum steps: ',maxsteps)

# Read the input and output filenames
cif_path_input=sys.argv[1]
cif_path_output=sys.argv[2]

# For test
#cif_path_input = "M:\MACE_projects\molecules\DMSO.cif"
#cif_path_output = "M:\MACE_projects\molecules\DMSO_wrap_P1-out.cif"

# Read cif in P1, Warning - wraps to  <0,1> interval
atoms = read(cif_path_input)

# Set the MACE calcualtor
calculator = MACECalculator(model_paths=model_with_path,dispersion=False,default_dtype=precision, device=device)
atoms.calc = calculator

# Otimize cell
ucf = FrechetCellFilter(atoms)

# Otimization alogorithm is BFGS
opt = BFGS(ucf)

# Run optimization
start = time.time()
opt.run(fmax=fmax,steps=maxsteps)
end = time.time()

# Print the energy and elapset time
energy = atoms.get_potential_energy()
print("Final energy (eV):", energy)
print('Time (s) ',end - start)

# Save data to CIF with second wrap disabled (so we match input)
atoms.write(cif_path_output,wrap = False)

Sample result demonstrating the exploding forces:

BFGS:   29 10:41:14   -73211.648438        0.175761
BFGS:   30 10:41:14   -73211.671875        0.221249
BFGS:   31 10:41:14   -73211.695312        0.285815
BFGS:   32 10:41:14   -73211.726562        0.284400
BFGS:   33 10:41:14   -73211.757812        0.278854
BFGS:   34 10:41:14   -73211.804688        0.342102
BFGS:   35 10:41:14   -73211.867188        0.387026
BFGS:   36 10:41:15   -73211.921875        0.419696
BFGS:   37 10:41:15   -73211.984375        0.435478
BFGS:   38 10:41:15   -73212.046875        0.443961
BFGS:   39 10:41:15   -73212.156250        0.673312
BFGS:   40 10:41:15   -73212.250000        0.836853
BFGS:   41 10:41:15   -73212.382812        1.009415
BFGS:   42 10:41:15   -73212.492188        1.087747
BFGS:   43 10:41:15   -73212.703125        1.147222
BFGS:   44 10:41:15   -73212.867188        1.045183
BFGS:   45 10:41:16   -73213.140625        0.993227
BFGS:   46 10:41:16   -73213.421875        1.521351
BFGS:   47 10:41:16   -73213.601562        3.248649
BFGS:   48 10:41:16   -73214.023438        2.468689
BFGS:   49 10:41:16   -73215.273438        4.405566
BFGS:   50 10:41:16   -73216.000000        4.585445
BFGS:   51 10:41:16   -73216.734375        6.393103
BFGS:   52 10:41:16   -73217.843750        9.628826
BFGS:   53 10:41:16   -73219.914062       15.039652
BFGS:   54 10:41:17   -73222.234375       23.162329
BFGS:   55 10:41:17   -73226.539062       36.280894
BFGS:   56 10:41:17   -73222.796875      126.252441
BFGS:   57 10:41:17   -73228.070312      146.598112
BFGS:   58 10:41:17   -73232.531250       88.004586
BFGS:   59 10:41:17   -73220.484375      100.359730
BFGS:   60 10:41:17   -73222.671875      171.820705
BFGS:   61 10:41:17   -73245.593750      347.042685
BFGS:   62 10:41:17   -73287.242188      434.991750
BFGS:   63 10:41:18   -73364.046875     1043.966117

[gabor1]

Mace-off can have some holes. Can you try a few dozen steps of the FIRE optimiser, before you turn on BFGS? Also I would use LBFGS from the PreconLBFGS ASE package , it's more stable.

[husakm]

I had tested all the mentioned optimizers and nothing works. The issue are nonsence forces from the FF. Looks like:
-C-SF5 group is not handled by any MACE
-NH3+ is not handlet OK by MACELES (not at all by the OFF as expected)
I will probably come with bigger list of molecular pattern not supported well ...

[ilyes319]

Any models can have holes. You might want to try mace-omol too, might be better behaved here.

[husakm]

I see the model have holes. I am able to detect the holes. But is there anybody able/open to repair this holes ?

[gabor1]

Try the mace-omol as Ilyes suggested. It is trained on a much larger set of molecules than MACEOFF

[gabor1]

I have successfully optimised ADAVAG using the mace-omol model (https://github.com/ACEsuit/mace-foundations/releases/tag/mace_omol_0), using PreconLBFGS. No sign of any holes here. Don't forget that MACE-OFF has not seen more than one molecule at all, whereas the omol models have seen lots of multi-molecule clusters. So this is not that surprising on the whole. The omol models are a whole lot better.

[husakm]

Thank you for the ADAVAG test. This system probably failed with OFF because of the O-, NH3+ ionic bond system.
I believe the MACE OFF24 middle was well tested on crystal structure prediction off non ionic pure organic crystals:
https://pubs.rsc.org/en/content/articlelanding/2025/cp/d5cp00593k
The results were perfect close to the DFT result (I had myself replicated the calculation myself with identical results).
So it look like it see surrounding molecules ...
Question:
Can the OMOL model handle charged/ionic system with long-range interaction ?
I believe MACE without the LES extension (MACELES) can not handle ionic long range interaction ...
Does the OMOL model already handle charges based on the LES extension ?

[gabor1]

You are fundamentally not wrong, but it's a bit more nuanced. For sure, only electrostatically aware models will do very well in the long run across a wide range of soft matter and molecular tasks. We have our own electrostatic model in the works, the first version will be released shortly. But for some tasks, e.g. many organic crystals, even short range models trained on huge datasets (like omol) are doing ok. It is because the majority of the interaction is short ranged. So I expect such models to give sensible results, not explode, etc. To get the highest accuracy that is needed to discriminate correctly between polymorphs, we will need explicit electrostatics. But also a very good treatment of dispersion (better than e.g. D3)

[gabor1]

The current omol model that I linked above is short ranged, but it is trained on such a large training set that I expect it to be stable for your problem.