BUG: `MAGMOM` not preserved in r2SCAN workflow

[ShawnTxY]

Describe the bug I encountered an issue with MAGMOM when defining my own user_incar_settings in atomate2 0.0.14 and pymatgen 2024.3.1. The MAGMOM values I defined in user_incar_settings were not correctly applied throughout all the steps in the workflow. Specifically, there might be a handoff issue where the MAGMOM setting for the r2SCAN INCAR file differ from the initial PBEsol MAGMOM values when transitioning from the PBEsol to the r2SCAN step using MPMetaGGADoubleRelaxMaker.

My initial MAGMOM initialization was to set the individual spin states for the two transition metals (both iron) with its individual spin states of Fe2+ = 3.0 and Fe3+ = 1.0.

Tagging @rkingsbury here.

Thank you all!

POSCAR:

Na8 Fe8 C24 N24
1.0
  10.500000000000000    0.0000000000000000    0.0000000000000000
   0.0000000000000000   10.500000000000000    0.0000000000000000
   0.0000000000000000    0.0000000000000000   10.500000000000000
Na Fe2+ Fe3+ C N
8 4 4 24 24

r2SCAN workflow:

from atomate2.vasp.flows.mp import MPMetaGGADoubleRelaxMaker
from pymatgen.core import Structure
from atomate2.vasp.powerups import update_user_incar_settings
from jobflow import run_locally, JobStore

structure = Structure.from_file('POSCAR')

user_incar_settings = {
        "LVTOT": "False",
        "MAGMOM": {"Fe2+":3.0, "Fe3+":1.0}
        }

relax_job = MPMetaGGADoubleRelaxMaker().make(structure)

relax_job_update = update_user_incar_settings(relax_job,user_incar_settings)

run_locally(relax_job_update, create_folders=True, log=True)

Observed behavior:

• For the PBEsol workflow, the MAGMOM in the output INCAR was set to MAGMOM = 8*0.6 4*3.0 4*1.0 48*0.6

• For the r2SCAN workflow, the MAGMOM in the output INCAR was set to MAGMOM = 64*0.0.

Expected behavior The expected MAGMOM for both PBEsol and r2SCAN should be set at MAGMOM = 8*0.6 4*3.0 4*1.0 48*0.6.

[Andrew-S-Rosen]

Tagging @esoteric-ephemera in case it impacts his ongoing atomate2 work.

[esoteric-ephemera]

I was able to replicate this using a similar flow. This specific behavior occurs because the oxidation state labels are not preserved after the first VASP calculation finishes.

I can see in the fireworks spec of all jobs in the flow that the user_incar_settings are appropriately updated, but subsequent jobs have oxidation state labels purged. VASP doesn't preserve site labels in CONTCAR / vasprun.xml

Another issue to consider here is that the MP flows carry the WAVECAR and CHGCAR through to subsequent calculations. Every subsequent job will take MAGMOMs from a previous wavefunction or charge density

[rkingsbury]

Thanks @esoteric-ephemera ! If I understand correctly then, there is no bug in atomate2 here; this is down to how VASP interacts with oxidation state labels. Two questions:

subsequent jobs have oxidation state labels purged. VASP doesn't preserve site labels in CONTCAR / vasprun.xml

So I guess this part of the bug is really a limitation of VASP, am I understanding that correctly?

Another issue to consider here is that the MP flows carry the WAVECAR and CHGCAR through to subsequent calculations. Every subsequent job will take MAGMOMs from a previous wavefunction or charge density

Does this mean that if CHGCAR or WAVECAR are present, the MAGMOM tag is effectively ignored? If the MAGMOM is properly initialized in the first calculation, then it (or however it evolves during the calculation) should be propagated to subsequent calculations? If that's correct, I guess this is a good thing, and the apparent 0 MAGMOM is a red herring?

[Andrew-S-Rosen]

Does this mean that if CHGCAR or WAVECAR are present, the MAGMOM tag is effectively ignored? If the MAGMOM is properly initialized in the first calculation, then it (or however it evolves during the calculation) should be propagated to subsequent calculations? If that's correct, I guess this is a good thing, and the apparent 0 MAGMOM is a red herring?

This is indeed correct. In such scenarios, the MAGMOM flag can be "misleading" to some degree. However, there is an edge case worth mentioning. If in Step 1 of your workflow you are running with the gamma-point only version of VASP (1x1x1 k-points) and in Step 2 of your workflow you are running with the standard version of VASP (>1x1x1 k-points), the standard version of VASP will crash because it can't read the gamma-point only version of VASP's WAVECAR. In doing so, Custodian will automatically remove the WAVECAR and re-launch VASP so it can proceed normally. However, if the MAGMOM flag is "wrong", then you may get unexpected behavior at this step since you are now starting from scratch, and VASP will initialize the magnetic moments from the INCAR.

[esoteric-ephemera]

@rkingsbury:

So I guess this part of the bug is really a limitation of VASP, am I understanding that correctly?

Yeah a limitation on the VASP side, but it should be possible to add some glue code that reassigns oxidation states at each workflow step. (Or maybe better, uses bader / OUTCAR to assign better charge states?) Is that something you'd be interested in having for workflows?

+1 on @Andrew-S-Rosen's explanation of the MAGMOMs / WAVECAR interplay. I think they similarly get overwritten if a CHGCAR is present

[Andrew-S-Rosen]

(Or maybe better, uses bader / OUTCAR to assign better charge states?)

Please don't 😅 Bader is actually awful for predicting oxidations states, especially in complex materials (don't even look at it for a MOF...). I would not rely on any population analysis for this, as the variance is huge. Here is one of many examples discussing the topic. Or the following figure from yours truly.

For this scenario, I'd suggest sticking with the user-defined formal oxidation states as a matter of consistency and following the principle of least surprise.

[rkingsbury]

Yeah a limitation on the VASP side, but it should be possible to add some glue code that reassigns oxidation states at each workflow step. (Or maybe better, uses bader / OUTCAR to assign better charge states?) Is that something you'd be interested in having for workflows?

My main concern is to avoid confusing users (as prompted this question from @ShawnTxY ). It seems clear that, for the r2SCAN workflow at least, the incorrect MAGMOM setting in the INCAR is innocuous but misleading, so having glue code that reassigns oxidation states would prevent that. It might also prevent surprises in edge cases like @Andrew-S-Rosen raised where, e.g., custodian removes the WAVECAR and VASP falls back to the MAGMOM in the INCAR.

Maybe an alternative stopgap would be to have the workflow remove the MAGMOM flag if a WAVECAR/CHGCAR is present?

@ShawnTxY , feel free to chime in if you have other ideas about how the workflow should behave in this situation.

[esoteric-ephemera]

There is a use case for applying MAGMOMs even when restarting from WAVECAR/CHGCAR, which is to enforce a specific magnetic symmetry group. To be clear, this is a quote from the VASP manual:

When restarting a magnetic calculation, MAGMOM is only used to determine the symmetry of the system and not to set the on-site magnetic moment. Therefore, if you remove the MAGMOM tag before restarting from a converged WAVECAR or CHGCAR, the magnetization is likely to be symmetrized away.

So you can still speed up a calculation using WAVECAR/CHGCAR, but use the magnetic space group specified by MAGMOM. The actual values of on-site magmoms will get overwritten in the first self-consistency step

All this to say I'd lean towards writing glue code to carry forward custom magmoms

[mkhorton]

Adding a note on Bader, since I do my best whenever it comes up :)

• Yes, it is not useful for assigning formal valences. Bader basin charge integration != oxidation states. This has been well demonstrated in literature.
  • However, I enjoy its simplicity and mathematical elegance, and the Bader basins do contain physical meaning; they are not useless nor unreliable quantities
• Charge != magnetic moment. I had good results testing Bader basin integration to assign magnetic moment based on integrating the spin density based on the charge density basins, however without a good benchmark to share I would take this only as an anecdote

For the very nice plots @Andrew-S-Rosen presented, I could easily believe there is a transformation between Bader and DDEC6 charges, and imagine both contain much of the same information.

In any case, this is a tangent to the present issue, just wanted to add the note so as not to mislead. I concur with @esoteric-ephemera in terms of issue resolution.

So you can still speed up a calculation using WAVECAR/CHGCAR, but use the magnetic space group specified by MAGMOM.

On this note, I do not know how if VASP handles this differently for collinear calculations from non-collinear calculations. For collinear calculations, the magnetic space group cannot be determined without assuming a preferred axis, so I don't know how exactly they handle this.