Unable to replicate bulk modulus of aluminum

pymatgen
1

Hello all, I am trying to calculate the bulk modulus of aluminum using pymatgen and a machine learned interatomic potential (MLIP) called CHGNet trained on the Materials Project database. I have my code below:

from chgnet.model.model import CHGNet
from pymatgen.core import Structure
from pymatgen.analysis.elasticity.strain import Deformation, DeformedStructureSet, Strain
from pymatgen.analysis.elasticity.stress import Stress
from pymatgen.analysis.elasticity import ElasticTensor

chgnet = CHGNet.load()
structure = Structure.from_file('cif_files/Al.cif')
prediction = chgnet.predict_structure(structure)

eq_stress = Stress(prediction['s'])

deformed_structure = DeformedStructureSet(structure)

strain_list = [i.green_lagrange_strain for i in deformed_structure.deformations]
stress_list = []

for i in range(len(deformed_structure)):
    prediction = chgnet.predict_structure(deformed_structure[i])
    stress_list.append(Stress(prediction['s']))

elastic_tensor = ElasticTensor.from_independent_strains(strain_list, stress_list, eq_stress)

print(elastic_tensor.property_dict)

The output I get is:

{'k_voigt': 9.845064083735359, 'k_reuss': 9.845063833881039, 'k_vrh': 9.845063958808199, 'g_voigt': 8.0023897488912, 'g_reuss': 4.5904549461157655, 'g_vrh': 6.2964223475034835, 'universal_anisotropy': 3.716336252207812, 'homogeneous_poisson': 0.23641618649298335, 'y_mod': 15569997014.898907}

I know the bulk modulus should be around 76GPa, but what I get from this is around 10GPa. What am I missing? Any help would be greatly appreciated! Many thanks

2

To me, this just sounds like CHGNet is not good at predicting the bulk modulus/other elastic properties of Aluminum. The values on the MP website for Aluminum are from DFT calculations, not from CHGNet.

3

I see. Yes I understand that the values on MP website are from DFT, but apparently CHGNet has been trained on values from MP, and I would assume this would have been a straightforward prediction. Just to be sure, nothing in the pymatgen code looks obviously incorrect?

4

You should probably relax the strained structure at fixed volume, rather than computing single points. See code snippet below, which gives me 36.2 GPa for the bulk modulus of fcc Al (mp-134).

To @matthewkuner’s point: it’s also quite possible that CHGNet doesn’t give an accurate bulk modulus.


from __future__ import annotations

from chgnet.model import CHGNet, StructOptimizer
import matplotlib.pyplot as plt
from mp_api.client import MPRester
import numpy as np
from pymatgen.analysis.eos import EOS
from pymatgen.core import Structure
from pymatgen.transformations.standard_transformations import (
    DeformStructureTransformation,
)


model = CHGNet.load()

def CHGNet_EOS(
    structure : Structure,
    linear_strain: tuple[float, float] = (-0.05, 0.05),
    number_of_frames: int = 6,
    eos_name :str = "birch_murnaghan"
):
    Nsites = structure.num_sites

    opt = StructOptimizer(model = model)
    relax0 = opt.relax(structure)
    s0 = relax0["final_structure"]
    e0 = model.predict_structure(s0)["e"]

    linear_strains = np.linspace(
        linear_strain[0], linear_strain[1], number_of_frames
    )
    deformations = [(np.identity(3) * (1 + eps)).tolist() for eps in linear_strains]

    EV_data = np.zeros((number_of_frames+1,2))
    EV_data[0] = [s0.volume, e0]
    for ideform in range(number_of_frames):
        deformation = DeformStructureTransformation(deformation=deformations[ideform])
        sdeform = deformation.apply_transformation(s0)
        srelax = opt.relax(sdeform,relax_cell=False)["final_structure"]
        erelax = model.predict_structure(srelax)["e"]
        EV_data[1+ideform] = [srelax.volume, erelax]

    EV_data = EV_data[np.argsort(EV_data[:,0]),:]
    EV_data[:,1] *= Nsites

    fig, ax = plt.subplots(figsize=(6,4))

    ax.scatter(EV_data[:,0],EV_data[:,1],color="darkblue")
    
    eos = EOS(eos_name=eos_name).fit(EV_data[:,0],EV_data[:,1])
    print(eos_name,{**eos.results, "b0 GPa": float(eos.b0_GPa)})
    vols = np.linspace(0.95*EV_data[:,0].min(),1.05*EV_data[:,0].max(),2000)
    interp = eos.func(vols)
    ax.plot(vols,interp,linestyle=':',label=eos_name,color="darkblue")

    ax.legend(fontsize=10)
        
    ax.set_xlim(vols[0],vols[-1])
    ax.set_ylim(
        min(EV_data[:,1].min(),interp.min()),
        max(EV_data[:,1].max(),interp.max())
    )
    ax.set_xlabel("Volume ($\\mathrm{\\AA{}}^3$)",fontsize=12)
    ax.set_ylabel("$E(V)$ (eV)",fontsize=12)

    plt.show()

if __name__ == "__main__":

    from os import path

    MPID = "mp-134"

    mp_struct_file = f"./POSCAR_{MPID}"
    if not path.isfile(mp_struct_file):
        with MPRester() as mpr:
            structure = mpr.summary.search(
                fields=["structure"],
                material_ids=[MPID]
            )[0].structure
        structure.to(mp_struct_file)
    else:
        structure = Structure.from_file(mp_struct_file)

    CHGNet_EOS(structure)
1 Like
5

I see, thank you for sharing the code snippet with me. When I run it using the cif file for aluminum that I have, I actually get a bulk modulus of -137 GPa! But besides that, for the calculation of elastic constants (bulk, shear, Poisson’s ratio, Young’s modulus), would you recommend I just use pymatgen as a standalone package? I am not very familiar with getting the stress from pymatgen using the code I have in the question, but I also wanted something that was faster.

Edit: I get 36.2 GPa as well, although if I explicitly download the cif file from MP and use that as the structure, I get -137 GPa. Will edit again once I know what this is about.

6

Confirming that I get -137 GPa for the bulk modulus when I load in the mp-134 CIF with Structure.from_file. Note however, if you pass primitive=True to from_file, or load the structure with pymatgen.io.cif.CifParser, you will again get 36.2 GPa.

This shouldn’t depend on whether the conventional or primitive cell is used, I’ll need to look into why that’s happening.

Pymatgen is intended to be modular, so no need to use it as a standalone. In this case, I think we’re seeing that the stresses from CHGNet are not sufficiently accurate for computation of the elastic tensor. Tagging in @janosh, who helped train CHGNet.

1 Like
7

Okay, that makes sense. Also, thank you for the clarification about the -137 GPa value. It would be good to use a dedicated DFT package to compute the stresses then and feed them into pymatgen. But I’m looking forward to seeing if this can be resolved within CHGNet itself as it would speed up things considerably