Unexpected Negative Mixing Enthalpy Using a Ga-In-N MEAM Potential

LAMMPS LAMMPS General Discussion
1

Hello everyone,

I am using a published Ga-In-N MEAM potential downloaded from NIST website ( Interatomic Potentials Repository ) to calculate thermodynamic properties of liquid Ga-In alloys (at 1000 K) in LAMMPS, and I have a question about the calculated mixing enthalpy.

I am using LAMMPS (10 Sep 2025). The relevant potential setup is:

pair_style meam
pair_coeff * * library.meam Ga In N GaInN.meam Ga In

The systems were first heated to 2000 K and then cooled to 1000 K. After reaching 1000 K, I performed a 2 ns equilibration run followed by a 20 ns production run. The production trajectory was checked using a block analysis, and no systematic energy drift was observed.

The mixing enthalpy was calculated as the difference between the enthalpy of the Ga-In alloy and the mole fraction weighted enthalpies of pure Ga and pure In, with all systems evaluated at the same temperature and pressure. We performed this calculation across the full composition range and obtained negative mixing enthalpies for all compositions. This is unexpected because both experimental thermodynamic data and the values reported in the original publication of this MEAM potential are slightly positive for the Ga-In system.

We also tested the 0 K properties of the potential. The bulk modulus was estimated by applying small volume perturbations around the equilibrium structure and fitting the resulting energy-volume relation. For the 50Ga-50In composition, the calculated bulk modulus was approximately 37 GPa, compared with 48 GPa reported in the original publication.

Has anyone here used this potential for thermodynamic calculations of Ga-In alloys? Any suggestions or insights would be greatly appreciated.

2

Two things I’d check first.

  • Make sure the pure Ga and pure In reference runs read GaInN.meam too (e.g. pair_coeff * * library.meam Ga In N GaInN.meam Ga), since the docs say the parameter file overrides or complements the library’s single-element values, so a NULL or library-only reference gives different pure energies.
  • And since your 0 K bulk modulus is already off (37 vs 48 GPa), see whether the paper’s settings for flags like ialloy, mixture_ref_t, emb_lin_neg, bkgd_dyn or erose_form match what’s in the file, because LAMMPS’s defaults don’t always match DYNAMO.

Let me know if that works.

1 Like
3

Hi Joash,

Thank you for your suggestions!

I used the same MEAM potential files for both pure Ga/In and the alloys.

I found that most of the flags have corresponding equations or descriptions in the original paper, except for ialloy=2 and zbl=0.

Should I test different values for these two flags and compare the results?

4

I wouldn’t change zbl. zbl(I,J) only blends the MEAM pair term with ZBL at small atom separations, so it shouldn’t touch a liquid at 1000 K or a bulk modulus taken near equilibrium. Leaving it at 0 as the file has it is fine.

ialloy is the one that could matter. The pair_meam doc lists 0 as standard averaging of the t parameters (the LAMMPS default), 1 as the alternative averaging used in DYNAMO, and 2 as no averaging, where each atom uses its own single-element t values. That choice only changes anything when an atom has neighbours of a different element. So your pure Ga and pure In references come out the same for any ialloy, and only the alloy energy moves. That’s exactly the term that sets the sign of your mixing enthalpy.

I’d test it, but treat it as a check on sensitivity rather than as a fit. Run the same 0 K bulk modulus calculation for 50Ga-50In with ialloy = 0, 1 and 2, keeping everything else fixed. That costs minutes rather than 20 ns. If one setting brings you close to the paper’s 48 GPa, that’s a strong hint about which averaging the authors actually used. If the paper reports a 0 K formation energy for an ordered Ga-In structure, comparing that too would be an even sharper check. Then it’s worth writing to the authors, or to the NIST repository contact for that potential, to confirm which code and settings they used. The file’s ialloy = 2 may not be what their fitting code did.

1 Like
5

Thank you very much for your comment. I just tested both ialloy = 0 and ialloy = 1, and found that the 0 K bulk modulus of 50Ga50In did not change significantly, remaining around 37–38 GPa.

6

Hmm…that rules out the averaging, which helps. I think the parameter file itself is the better suspect now.

More than one version of the Ga-In cross terms is going around. The files from the Lee group’s site, which is where NIST got them in 2020, have Ec(1,2) = 4.06. NIST’s own parameter check flags Ga-In as “delta_Ec off by 10” (interactions.csv in https://github.com/lmhale99/meam2nn): 4.06 is the Ga/In average (2.705) plus 1.355 instead of 0.1355. The copy in NIST’s potentials library has Ec(1,2) = 2.8405, but it kept re(1,2) = 2.8276 and alpha(1,2) = 3.87486077 from the old file. Could you check which values your GaInN.meam has?

That alpha was worked out for Ec = 4.06. For the B1 reference, B = alpha² Ec / (9 re³). That gives 48.0 GPa with Ec = 4.06 and 33.6 GPa with 2.8405. So if your file pairs 2.8405 with that alpha, the potential can’t give the paper’s 48 GPa, and the drop is in the same direction you’re seeing. NIST’s parameter database lists alpha(1,2) = 4.6326 for Ga-In. That’s what you get if you keep 48 GPa and use the corrected Ec.

A quick way to confirm: build B1 (rock salt) GaIn with a = 2 × re(1,2) = 5.655 Å and run the same 0 K check. Energy per atom should come out close to -Ec(1,2), and B close to the formula above. Then change only alpha(1,2) to 4.6326 and repeat. If B moves to about 48 GPa, rerun one liquid composition (50:50 is enough) and see whether the mixing enthalpy moves toward the paper’s.

I’d also compare Ec(1,2), re(1,2) and alpha(1,2) with the Ga-In row of the paper’s B1 reference table, including the sign of ΔEc, since that sign decides which way the reference pushes the mixing. Whatever you find is worth sending to the NIST repository contact, since other people are downloading the same file.

8

Yes, after I set alpha(1,2) to 4.6326, the bulk modulus of B1 Ga50In50 at 0K is 48.11 GPa, I will further test if Hmix is better at 1000K with this new alpha value and update to you. And could I ask where I can find the NIST’s parameter database, is it included in the link you shared? I want to double check if other flags are the same with the MEAM I’m using. Thank you!

9

It’s in the same repo under meam2nn/parameters: element.csv has the single elements, alloy.csv the pair terms, and ternary.csv the three-body Cmin/Cmax. One correction: the Ga-In row doesn’t list alpha. It gives delta_Ec = −0.1355, re = 2.8276 and B = 0.2996 eV/ų (48.0 GPa), and 4.6326 is the alpha that B implies with Ec = 2.8405, which matches your 48.11 GPa. The row also has d+/d− (0.0735), the four Cmin values, Cmax 2.8 and rc 5.2, and the Ga-In-N row in ternary.csv has the three-body Cmin (watch the index order). LAMMPS flags like ialloy and mixture_ref_t aren’t in there, since they aren’t fitted parameters. The README calls the repo a work in progress, so where it disagrees with the paper, go with the paper.

Honestly, you should try out Verlet (https://www.verlet.app/) for these runs (full disclosure, I’m building it). You can keep the 0 K B1 check, the pure Ga and In references and the liquid runs together, each in its own folder, so your alpha = 3.87 and 4.6326 versions never overwrite each other. Then compare them side by side to see which change actually moved the mixing enthalpy. It also checks your MEAM pair_coeff element mapping before anything goes to the queue. It wouldn’t have caught the Ec/alpha mismatch itself, since that took reading the file against the paper. It’s in beta at verlet.app if you want to try it.