Issues with elastic constants predicted by chgnet

[chtchelkatchev]

Thank you for the ChgNet program! I plan to continue using it and would like to discuss one of the issues. Specifically, I noticed that the atomate2 utility from the elastic.ipynb example produces incorrect calculations of elastic modules. When I start "https://github.com/JaGeo/TutorialAtomate2Forcefields/blob/main/atomate2_workflow_tutorial/elastic.ipynb" with the ChgNet-3.0 force field, I receive non-physical values for elastic constants. Namely for Si: {'k_voigt': -13943156.074867852, 'k_reuss': -13943156.229127234, 'k_vrh': -13943156.151997544, 'g_voigt': -8365897.858119255, ...}
When I run the program again, different incorrect results are obtained. It appears that the answer for the elastic constants is a collection of random numbers, which is not what should happen.

Using the same example-program with MACE force field I get the reasonable stable answer:
{'k_voigt': 74.7004463937039, 'k_reuss': 74.70044639370386, 'k_vrh': 74.70044639370388, 'g_voigt': 27.183833551852977, 'g_reuss': 26.23384316388329, 'g_vrh': 26.70883835786813, 'universal_anisotropy': 0.18106199347825047, 'homogeneous_poisson': 0.3402646259018905, 'y_mod': 71593822499.96439, 'trans_v': 3409.3623514873443, 'long_v': 6928.791534089603, 'snyder_ac': 61.4081610788213, 'snyder_opt': 0.21207312240521548, 'snyder_total': 61.620234201226516, 'clark_thermalcond': None, 'cahill_thermalcond': 1.02856605617899, 'debye_temperature': 417.8290081748565}

Please, help me figure out the problem. Hope, this is not a fatal Chgnet error. Program used is attached.
elastic.zip

[chtchelkatchev]

I asked my colleagues to independently test the calculation of elastic moduli in chgnet. Using the LAMMPS interface and standard LAMMPS approaches, the elastic moduli of aluminum were calculated. The answer for elastic modules turned out to be meaningless. When the force field in Lammps was changed from chgnet to "eam", the elastic modules became reasonable. This suggests that there may be serious issues with CHGNet.

[naveenmohandas]

I can also support this, I tried with the given zip file and observed similar errors. I tried with BCC Fe, it appears to give reasonable values most of the time but occasionally it gives random values which are orders of magnitude large. When I use Matcalc which was provided in the old examples of CHGNet for calculating elastic constants, it gives consistent values as far as I tried.

from chgnet.model import CHGNetCalculator, CHGNet
from pymatgen.core import Structure
from matcalc.elasticity import ElasticityCalc


calc = CHGNetCalculator(model=CHGNet.load())                                                            

Pure_bulk = Structure(                                                                                                                                                       
    lattice=[[0, 2.73, 2.73], [2.73, 0, 2.73], [2.73, 2.73, 0]],                                                                                                             
    species=["Si", "Si"],                                                                                                                                                    
    coords=[[0, 0, 0], [0.25, 0.25, 0.25]],                                                                                                                                  
) 
calc_elastic = ElasticityCalc(calculator=calc)                                                           
elastic_prop = calc_elastic.calc(Pure_bulk)                                                              
elastic_tensor = elastic_prop['elastic_tensor'] 
print("bulk modulus vrh",elastic_prop['bulk_modulus_vrh']*160.21766208) # convert to GPa from eV/Angstrom^3
# print(elastic_prop)


c_11 = elastic_tensor[0,0,0,0]*160.21766208
c_44 = elastic_tensor[1,2,1,2]*160.21766208
c_44_or = elastic_tensor[2,1,2,1] * 160.21766208
c_12 = elastic_tensor[0,0,1,1]*160.21766208
print(f"{c_11=},{c_12=}, {c_44=}, {c_44_or=}")