Difficulty reproducing Henriksson 2013 Fe–Cr–C potential benchmarks with LAMMPS tersoff/zbl

LAMMPS
1

Hello,

I am trying to implement the full Fe–Cr–C bond-order potential developed by Henriksson et al.:

K. O. E. Henriksson, C. Björkas, K. Nordlund,
Atomistic simulations of stainless steels: a many-body potential for the Fe–Cr–C system,
J. Phys.: Condens. Matter 25 (2013) 445401.

My goal is eventually to use this potential for Fe–Cr solid solutions and Cr23C6 precipitates/interfaces in LAMMPS.

The publicly available NIST file contains only the Fe–C subset, so I reconstructed the remaining Cr-containing terms from the published parameters.

I first used the straightforward mapping of the published ABOP parameters to tersoff/zbl. The Fe–C, Cr–Cr, and Fe–Cr parts behaved well, but the Cr–C compounds did not reproduce the published benchmarks.

After contacting researchers familiar with the original implementation, I learned that there are subtle differences in the three-body index conventions between the original PARCAS implementation and LAMMPS. I therefore revised the mapping of the Cr–C three-body terms.

With the revised mapping, the equilibrium lattice constants improved substantially:

CrC-B2:

LAMMPS = 2.40504 Å

paper = 2.405 Å

CrC-B1:

LAMMPS = 3.95392 Å

paper = 3.904 Å

Cr23C6:

LAMMPS = 10.41719 Å

paper = 10.431 Å

So the equilibrium structures, especially CrC-B2 and Cr23C6, are reproduced quite well.

However, the bulk moduli still show noticeable discrepancies:

CrC-B2:

LAMMPS = 595.5 GPa

paper = 537 GPa

CrC-B1:

LAMMPS = 401.4 GPa

paper = 424 GPa

Cr23C6:

LAMMPS ≈ 305 GPa

paper = 358 GPa

For Cr23C6 I checked both:

  1. E(V) scans with internal atomic coordinates relaxed at every volume;

  2. E(V) scans with fixed internal fractional coordinates.

Both gave approximately 305 GPa, so the discrepancy does not appear to come from internal-coordinate relaxation.

For CrC-B2, which has essentially no internal positional degree of freedom, the lattice constant is reproduced almost exactly while the bulk modulus is still about 11% higher than the published value.

This makes me suspect that the original PARCAS/TULIP implementation may not be completely equivalent to the standard LAMMPS tersoff/zbl implementation, even when the same nominal parameters are used.

My questions are:

  1. Has anyone implemented the complete Henriksson Fe–Cr–C potential in LAMMPS and successfully reproduced the CrC / Cr23C6 benchmarks?

  2. Is there a complete Fe–Cr–C LAMMPS parameter file available somewhere beyond the published Fe–C NIST subset?

  3. Are there known implementation differences between the original Albe/ABOP form used in PARCAS/TULIP and LAMMPS tersoff/zbl that could affect second derivatives such as elastic constants or bulk moduli while leaving equilibrium lattice constants nearly unchanged?

  4. In particular, could differences in the ZBL/repulsive-potential blending or three-body indexing produce this behavior?

I would be very grateful for any advice, reference input files, or previous implementations that I could compare against.

Thank you very much.