Persistent SCF non-convergence in a collinear DFT+U calculation (CeIn₃ + µ⁺)

Date: 2026-08-06
System: CeIn₃ with an interstitial positive muon, 2×2×2 supercell
Code: Quantum ESPRESSO 7.5 (pw.x), via AiiDA 2.7.3
Problem: Phase 5, AFM SCF (collinear): Tags Ce FM planes, AFM stacking, runs collinear SCF, angle1/angle2 set to encode the [1,1,1]/√3 moment direction, does not converge, and get traped between 2 values oscillating between them.

If you would like, you can comment directly on this file, then save the HTML with the changes.

I want first to set a definition for the terms collinear and non-collinear, it's worth being precise about what that means here. CeIn₃'s real magnetic structure genuinely is collinear in the physics sense: every Ce moment is either exactly parallel or exactly antiparallel to a single axis, [111]; no canting between sites. On the other hand, and please correct me if I'm wrong, QE's standard collinear mode (nspin=2) can only put that axis along z of the simulation cell. Since the cubic cell here isn't rotated to align [111] with z, the only way to point the moments along the true [111] direction without distorting the cell is to use the non-collinear line (noncolin=True), which allows an arbitrary quantization axis per atom via angle1/angle2; even though there's no actual non-collinearity (spin canting) between the Ce sites themselves. It's being used as a "collinear-but-tilted" trick, not because the magnetism is truly non-collinear.

I will use the term collinear to describe our system, but you will see non-collinear in the codes.

What we are trying to compute

We are simulating a µSR experiment on CeIn₃. The pipeline:

1. In phase 3, we relax a positive muon (modeled as a charged interstitial H atom, tot_charge=1.0) into a candidate stopping site inside a 2×2×2 (and separately 3×3×3, 4×4×4, etc.) supercell of CeIn₃; this step (structural relaxation, no Hubbard U) converges without issue.
2. The step we cannot get past: in phase 5, we take that relaxed geometry and run a fully self-consistent, collinear, DFT+U calculation to obtain the type-II AFM (111) spin density of the CeIn₃ lattice, from which we later extract the Fermi contact field and dipolar field the muon would experience.

Step 2 is the one described in this letter.

The problem, precisely

Target: conv_thr = 1.0e-5 Ry (energy accuracy).

Observed: after electron_maxstep = 300 iterations, the SCF accuracy reaches a flat, oscillating plateau and goes no further. The plateau level depends somewhat on the numerical settings, and our best result to date (combining our three strongest levers, Round 8) reaches ~0.0026 Ry at its lowest point, still roughly 260× above the target, essentially unchanged over the last ~200 of 300 iterations in every case, across 18 tests so far. This is not slow convergence that more iterations would fix; it is stcuk.

Recurring, throughout every run: c_bands: N eigenvalues not converged warnings, QE's iterative diagonalization repeatedly failing to fully resolve some eigenvalues.

The exact calculation being run

This is pw.x input file for our Round 6, pulled directly from the AiiDA-submitted job (verdi calcjob inputcat, not reconstructed):

&CONTROL
  calculation = 'scf'
  disk_io = 'high'
  outdir = './out/'
  prefix = 'aiida'
  pseudo_dir = './pseudo/'
  verbosity = 'high'
/
&SYSTEM
  angle1(1) =   5.4735600000d+01
  angle1(2) =   1.2526440000d+02
  angle2(1) =   4.5000000000d+01
  angle2(2) =   2.2500000000d+02
  degauss =   2.0000000000d-02
  ecutrho =   4.0000000000d+02
  ecutwfc =   5.0000000000d+01
  ibrav = 0
  lspinorb = .false.
  nat = 33
  noncolin = .true.
  nosym = .true.
  ntyp = 4
  occupations = 'smearing'
  smearing = 'marzari-vanderbilt'
  starting_charge(3) =   1.0000000000d+00
  starting_magnetization(1) =   7.0000000000d-01
  starting_magnetization(2) =   7.0000000000d-01
  tot_charge =   1.0000000000d+00
/
&ELECTRONS
  conv_thr =   1.0000000000d-05
  diagonalization = 'cg'
  electron_maxstep = 300
  mixing_beta =   5.0000000000d-02
  mixing_mode = 'plain'
  mixing_ndim = 12
/
ATOMIC_SPECIES
Ce1    140.116 Ce.paw.z_12.atompaw.wentzcovitch.v1.2.upf
Ce2    140.116 Ce.paw.z_12.atompaw.wentzcovitch.v1.2.upf
H      1.00794 H.pbe-rrkjus_psl.1.0.0.UPF
In     114.818 In.pbe-dn-rrkjus_psl.0.2.2.UPF
ATOMIC_POSITIONS angstrom
Ce1         -0.0000014650       0.0000042047       0.0113245327
In           2.3404671397       2.3408373457      -0.0111667241
In           2.3538760775       0.0000581675       2.3445364702
In          -0.0000036259       2.3541061976       2.3443660555
Ce2         -0.0000098359      -0.0000004291       4.6773105022
In           2.3405845294       2.3405313325       4.6998251742
In           2.3155787117       0.0000443012       7.0334727083
In          -0.0000003859       2.3158148213       7.0336226752
Ce2         -0.0000201594       4.6890071500      -0.0082627498
In           2.3404719781       7.0372023737      -0.0110411675
In           2.3417474416       4.6889745047       2.3444545067
In           0.0000287709       7.0239011425       2.3444070696
Ce1         -0.0000185428       4.6890031431       4.6971554190
In           2.3405833673       7.0374758630       4.6998334062
In           2.2711762388       4.6889918828       7.0334507716
In           0.0000315026       7.0621911002       7.0336630146
Ce2          4.6889825995      -0.0000003553      -0.0084447021
In           7.0375046123       2.3407860165      -0.0111979350
In           7.0241422308      -0.0000573656       2.3444953475
In           4.6890693235       2.3420197680       2.3447781119
Ce1          4.6889907413       0.0000064449       4.6974381857
In           7.0373938781       2.3405208812       4.6998308457
In           7.0623700179      -0.0000591905       7.0334459281
In           4.6888939299       2.2716128524       7.0339250884
Ce1          4.6889849904       4.6890111328      -0.0613078284
In           7.0375052009       7.0371992744      -0.0110929978
In           7.0362898101       4.6890012970       2.3444910565
In           4.6890277687       7.0359910580       2.3448146420
Ce2          4.6889879501       4.6889983518       4.7515777577
In           7.0373868083       7.0374775740       4.6998291926
In           7.1068099099       4.6890289590       7.0334725756
In           4.6888505429       7.1063835364       7.0339411769
H            4.6893179426       4.6889366638       7.0325518903
K_POINTS automatic
1 1 1 0 0 0
CELL_PARAMETERS angstrom
      9.3780000000       0.0000000000       0.0000000000
      0.0000000000       9.3780000000       0.0000000000
      0.0000000000       0.0000000000       9.3780000000
HUBBARD	ortho-atomic
 U	Ce1-4f	5.0
 U	Ce2-4f	5.0

Everything we have tried, and the actual results

All tests below use the same input structure (one specific relaxed geometry, from our edge_exact branch) so results are directly comparable. Each row changes only what is stated relative to the row above it (Round 3 is our informal "baseline"), except where noted as "combined."

Even our best result (Round 8, ~0.0026 Ry at its lowest point) is an improvement on the default, but is still 260× above the 1e-5 Ry target.

Representative raw output (Round 6, the input shown in the table above)

Early iterations:

 estimated scf accuracy    <      39.74242776 Ry   (iter 1)
 estimated scf accuracy    <      25.85845160 Ry
 estimated scf accuracy    <       7.97410865 Ry
 estimated scf accuracy    <       4.68543398 Ry
 ...
 estimated scf accuracy    <       0.42090813 Ry   (iter ~26)
 estimated scf accuracy    <       0.40048649 Ry
 estimated scf accuracy    <       0.39296991 Ry

Final iterations (295–300 of 300), flat, oscillating, no net progress:

 iteration #295   total energy = -7335.54430615 Ry   scf accuracy < 0.00799950 Ry
   total magnetization = 0.35  0.76  -1.62 Bohr mag/cell   absolute magnetization = 21.95
 iteration #296   total energy = -7335.54474091 Ry   scf accuracy < 0.00817787 Ry
   total magnetization = 0.40  0.74  -1.64 Bohr mag/cell   absolute magnetization = 21.93
 iteration #297   total energy = -7335.53690104 Ry   scf accuracy < 0.00827333 Ry
   total magnetization = 0.35  0.76  -1.62 Bohr mag/cell   absolute magnetization = 22.01
 iteration #298   total energy = -7335.54044460 Ry   scf accuracy < 0.00813764 Ry
   total magnetization = 0.38  0.76  -1.64 Bohr mag/cell   absolute magnetization = 21.91
 iteration #299   total energy = -7335.54304607 Ry   scf accuracy < 0.00813226 Ry
   total magnetization = 0.40  0.75  -1.65 Bohr mag/cell   absolute magnetization = 21.87
 iteration #300   total energy = -7335.52879847 Ry   scf accuracy < 0.00839354 Ry
   total magnetization = 0.36  0.80  -1.62 Bohr mag/cell

 convergence NOT achieved after 300 iterations: stopping

Representative raw output (Round 8, our current best result)

Same flat-plateau at a lower level, and visibly noisier iteration-to-iteration (consistent with the larger mixing_beta). Final 4 of 300 iterations, pulled directly from the retrieved aiida.out:

 iteration #297   total energy = -7335.96478982 Ry   scf accuracy < 0.00344255 Ry
   total magnetization = -1.96  0.97  3.08 Bohr mag/cell   absolute magnetization = 28.71
 iteration #298   total energy = -7335.96591631 Ry   scf accuracy < 0.00346373 Ry
   total magnetization = -1.97  0.96  3.09 Bohr mag/cell   absolute magnetization = 28.72
 iteration #299   total energy = -7335.97108559 Ry   scf accuracy < 0.00351508 Ry
   total magnetization = -1.97  0.97  3.09 Bohr mag/cell   absolute magnetization = 28.72
 iteration #300   total energy = -7335.97839010 Ry   scf accuracy < 0.00307672 Ry
   total magnetization = -1.97  0.97  3.09 Bohr mag/cell   absolute magnetization = 28.71

 convergence NOT achieved after 300 iterations: stopping

The full output is at: https://muathhamidi.github.io/muath.hamidi/documents/muSR_Aiida_QE/CeIn3_Aiida_QE_5_Output/CeIn3_Aiida_QE_5_Output.html

Where are we now

Now with all the canaries: three numerical levers improved the convergence, higher starting_magnetization (7B, the largest single effect), a more aggressive mixing_beta combined with plain mixing (7D), and a higher Hubbard Ueff (7K). Combined (Round 8), they gave only a modest further gain over 7B alone.

I'm thinking of these 2 candidates too, would like to hear your opintion about:

1. Something about the imposed non-collinear magnetic structure itself (angle1/angle2, the (111)-plane AFM) is responsible, independent of U, e.g. a near-degeneracy between two very close but distinct non-collinear spin configurations that both satisfy our imposed angle1/angle2 constraints loosely, which the SCF alternates between. This might explain why U=0 doesn't help. We have not tested any change to the magnetic-configuration itself. I'm not sure what to do here, redesigning the imposed magnetic configuration itself changes the way we look at the physcis of the system, not only a numerical issue.
2. The 1e-5 Ry threshold itself may be unreasonably tight for this level of theory (collinear + Hubbard U). I've seen people using 1e-8 Ry, which is even tighter (we currently use the looser 1e-5). I do not have independent guidance on what a realistic, defensible threshold is for a calculation of this type.

This is the main thing we would value your judgment on.

My main questions would be, is our issue numerical or physical? And how to move forward? Testing more numerical setups or go ahead and do something about the spin structure? And do you have an intuition how? Is there a defualt way to deal with this type of systems?

Codes

Current production submission script (scripts/phase5_afm_scf.py)

This is the file the pipeline calls in normal operation. Note it currently reflects our Round 4 settings (mixing_mode='local-TF', mixing_ndim=16, degauss=0.03, diago_full_acc=True). We deliberately do not overwrite this file with each new experiment.

import sys
import numpy as np
from aiida.engine import submit
from aiida.orm import load_node, load_code, load_group, Dict, StructureData
from aiida.plugins import DataFactory
from aiida_quantumespresso.common.hubbard import Hubbard, HubbardParameters
from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData
from config import LATTICE_CONST, PW_CODE, PSEUDO_GROUP_SSSP, hpc_resources_afm_scf


if __name__ == '__main__':
    PHASE_META = {
        'title': 'Phase 5: AFM SCF (Non-Collinear)',
        'goal': 'Compute the (111) antiferromagnetic spin density of the relaxed CeIn3+μ supercell using non-collinear DFT.',
        'overview': (
            'Loads the relaxed supercell from Phase 3, tags alternating (111) planes of Ce atoms '
            'as Ce1 (spin ↑) and Ce2 (spin ↓) to model the type-II AFM order observed in CeIn3. '
            'Full SOC is omitted because: (1) spin orientation along [111] is imposed explicitly via '
            'angle1/angle2, not derived from SOC anisotropy; (2) Phase 9 overrides the DFT moment with '
            'the experimental 0.48 µ_B anyway, so SOC precision on the moment is wasted compute; '
            '(3) Note: fully relativistic pseudopotentials drastically worsen Ce 4f SCF convergence. '
            'Submits a non-collinear pw.x SCF without spin-orbit coupling (lspinorb=False). '
            'Scales HPC resources with supercell size N using hpc_resources_afm_scf(N). '
            'The wavefunction output is retrieved by Phase 6 to run pp.x for spin density extraction.'
        ),
        'hpc': True,
        'hpc_resources': (
            f'Scaled by supercell size N via hpc_resources_afm_scf(N): '
            f'{hpc_resources_afm_scf(2)[0]} nodes/{hpc_resources_afm_scf(2)[1]//3600}h at N=2 '
            f'up to {hpc_resources_afm_scf(5)[0]} nodes/{hpc_resources_afm_scf(5)[1]//3600}h at N≥5, '
            f'× 32 MPI cores/node.'
        ),
    }

    # 1. Load the perfectly relaxed supercell from the output trajectory
    pk = int(sys.argv[1])
    relax_calc = load_node(pk)
    ase_struct = relax_calc.inputs.structure.get_ase()
    ase_struct.set_positions(relax_calc.outputs.output_trajectory.get_positions()[-1])

    # Derive supercell size N from cell dimensions.
    a_primitive = LATTICE_CONST
    cell_lengths = ase_struct.get_cell().lengths()
    N = max(1, round(cell_lengths[0] / a_primitive))   # 2 for 2×2×2, 3 for 3×3×3
    print(f"Detected supercell: {N}×{N}×{N}")

    scaled = ase_struct.get_scaled_positions()
    for i, atom in enumerate(ase_struct):
        if atom.symbol == 'Ce':
            x, y, z = scaled[i]
            plane = int(np.round(N * (x + y + z)))
            atom.tag = 1 if plane % 2 == 0 else 2   # Ce1 vs Ce2

    afm_structure = StructureData(ase=ase_struct)

    hubbard_params = []
    for i, site in enumerate(afm_structure.sites):
        if site.kind_name.startswith('Ce'):
            hubbard_params.append(
                HubbardParameters(
                    atom_index=i,
                    atom_manifold='4f',
                    neighbour_index=i,
                    neighbour_manifold='4f',
                    translation=(0, 0, 0),
                    # [ASSUMPTION] Hubbard Ueff = 5.0 eV for Ce 4f states, standard for CeIn3.
                    value=5.0,
                    hubbard_type='Ueff'
                )
            )
    hubbard = Hubbard(parameters=hubbard_params)
    afm_hubbard_structure = HubbardStructureData.from_structure(afm_structure, hubbard)

    builder = load_code(PW_CODE).get_builder()
    builder.structure = afm_hubbard_structure
    builder.pseudos = load_group(PSEUDO_GROUP_SSSP).get_pseudos(structure=afm_hubbard_structure)

    builder.kpoints = DataFactory('core.array.kpoints')()
    # [ASSUMPTION] Gamma-point only SCF.
    builder.kpoints.set_kpoints_mesh([1, 1, 1])

    # 3. Non-collinear (111) AFM SCF parameters
    builder.parameters = Dict(dict={
        'CONTROL': {
            'calculation': 'scf',
            'disk_io': 'high'
        },
        'SYSTEM': {
            'ecutwfc': 50.0, 'ecutrho': 400.0,
            'occupations': 'smearing', 'smearing': 'marzari-vanderbilt', 'degauss': 0.03,
            'lspinorb': False, 'noncolin': True, 'nosym': True,
            'starting_magnetization': {'Ce1': 0.5, 'Ce2': 0.5},
            'angle1': {'Ce1': 54.7356, 'Ce2': 125.2644},
            'angle2': {'Ce1': 45.0, 'Ce2': 225.0},
            'tot_charge': 1.0,
            'starting_charge': {'H': 1.0}
        },
        'ELECTRONS': {
            'conv_thr': 1.0e-5,
            'mixing_beta': 0.05,
            'electron_maxstep': 300,
            'diagonalization': 'cg',
            'diago_full_acc': True,
            'mixing_mode': 'local-TF',
            'mixing_ndim': 16
        }
    })

    # 4. Scale HPC resources with supercell size using config table.
    num_machines, max_wallclock_seconds = hpc_resources_afm_scf(N)

    builder.settings = Dict(dict={'cmdline': ['-nk', '1']})
    builder.metadata.options = {
        'resources': {'num_machines': num_machines, 'num_mpiprocs_per_machine': 32},
        'max_wallclock_seconds': max_wallclock_seconds,
        'custom_scheduler_commands': '#SBATCH --mem=64G\n',
        'withmpi': True,
    }

    node = submit(builder)
    print(f"STEP 5 (AFM): {N}×{N}×{N} AFM SCF Submitted! PK: {node.pk}")

Round 6 script (scripts/phase5_afm_scf_round6.py)

Differs only in the SYSTEM/ELECTRONS blocks (mixing_mode='plain', starting_magnetization=0.7, degauss=0.02, no diago_full_acc, mixing_ndim=12).

import sys
import numpy as np
from aiida.engine import submit
from aiida.orm import load_node, load_code, load_group, Dict, StructureData
from aiida.plugins import DataFactory
from aiida_quantumespresso.common.hubbard import Hubbard, HubbardParameters
from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData
from config import LATTICE_CONST, PW_CODE, PSEUDO_GROUP_SSSP, hpc_resources_afm_scf

if __name__ == '__main__':
    pk = int(sys.argv[1])
    relax_calc = load_node(pk)
    ase_struct = relax_calc.inputs.structure.get_ase()
    ase_struct.set_positions(relax_calc.outputs.output_trajectory.get_positions()[-1])

    a_primitive = LATTICE_CONST
    cell_lengths = ase_struct.get_cell().lengths()
    N = max(1, round(cell_lengths[0] / a_primitive))
    print(f"Detected supercell: {N}x{N}x{N}")

    scaled = ase_struct.get_scaled_positions()
    for i, atom in enumerate(ase_struct):
        if atom.symbol == 'Ce':
            x, y, z = scaled[i]
            plane = int(np.round(N * (x + y + z)))
            atom.tag = 1 if plane % 2 == 0 else 2

    afm_structure = StructureData(ase=ase_struct)

    hubbard_params = []
    for i, site in enumerate(afm_structure.sites):
        if site.kind_name.startswith('Ce'):
            hubbard_params.append(
                HubbardParameters(
                    atom_index=i, atom_manifold='4f', neighbour_index=i,
                    neighbour_manifold='4f', translation=(0, 0, 0),
                    value=5.0, hubbard_type='Ueff'
                )
            )
    hubbard = Hubbard(parameters=hubbard_params)
    afm_hubbard_structure = HubbardStructureData.from_structure(afm_structure, hubbard)

    builder = load_code(PW_CODE).get_builder()
    builder.structure = afm_hubbard_structure
    builder.pseudos = load_group(PSEUDO_GROUP_SSSP).get_pseudos(structure=afm_hubbard_structure)
    builder.kpoints = DataFactory('core.array.kpoints')()
    builder.kpoints.set_kpoints_mesh([1, 1, 1])

    builder.parameters = Dict(dict={
        'CONTROL': {'calculation': 'scf', 'disk_io': 'high'},
        'SYSTEM': {
            'ecutwfc': 50.0, 'ecutrho': 400.0,
            'occupations': 'smearing', 'smearing': 'marzari-vanderbilt', 'degauss': 0.02,
            'lspinorb': False, 'noncolin': True, 'nosym': True,
            'starting_magnetization': {'Ce1': 0.7, 'Ce2': 0.7},
            'angle1': {'Ce1': 54.7356, 'Ce2': 125.2644},
            'angle2': {'Ce1': 45.0, 'Ce2': 225.0},
            'tot_charge': 1.0,
            'starting_charge': {'H': 1.0}
        },
        'ELECTRONS': {
            'conv_thr': 1.0e-5,
            'mixing_beta': 0.05,
            'electron_maxstep': 300,
            'diagonalization': 'cg',
            'mixing_mode': 'plain',
            'mixing_ndim': 12,
        }
    })

    num_machines, max_wallclock_seconds = hpc_resources_afm_scf(N)
    builder.settings = Dict(dict={'cmdline': ['-nk', '1']})
    builder.metadata.options = {
        'resources': {'num_machines': num_machines, 'num_mpiprocs_per_machine': 32},
        'max_wallclock_seconds': max_wallclock_seconds,
        'custom_scheduler_commands': '#SBATCH --mem=64G\n',
        'withmpi': True,
    }

    node = submit(builder)
    print(f"ROUND 6 (plain mix + higher magmom): {N}x{N}x{N} AFM SCF Submitted! PK: {node.pk}")

7 parallel-sweep script (scripts/phase5_afm_scf_round7_sweep.py)

Defining all 12 variants tabulated:

import sys
import copy
import numpy as np
from aiida.engine import submit
from aiida.orm import load_node, load_code, load_group, Dict, StructureData
from aiida.plugins import DataFactory
from aiida_quantumespresso.common.hubbard import Hubbard, HubbardParameters
from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData
from config import LATTICE_CONST, PW_CODE, PSEUDO_GROUP_SSSP, hpc_resources_afm_scf

if __name__ == '__main__':
    pk = int(sys.argv[1])  # Phase 3 relaxed structure PK, e.g. 8581 (edge_exact)
    relax_calc = load_node(pk)
    ase_struct = relax_calc.inputs.structure.get_ase()
    ase_struct.set_positions(relax_calc.outputs.output_trajectory.get_positions()[-1])

    a_primitive = LATTICE_CONST
    cell_lengths = ase_struct.get_cell().lengths()
    N = max(1, round(cell_lengths[0] / a_primitive))
    print(f"Detected supercell: {N}x{N}x{N}")

    scaled = ase_struct.get_scaled_positions()
    for i, atom in enumerate(ase_struct):
        if atom.symbol == 'Ce':
            x, y, z = scaled[i]
            plane = int(np.round(N * (x + y + z)))
            atom.tag = 1 if plane % 2 == 0 else 2

    afm_structure = StructureData(ase=ase_struct)

    def build_hubbard_structure(ueff):
        hubbard_params = []
        for i, site in enumerate(afm_structure.sites):
            if site.kind_name.startswith('Ce'):
                hubbard_params.append(
                    HubbardParameters(
                        atom_index=i, atom_manifold='4f', neighbour_index=i,
                        neighbour_manifold='4f', translation=(0, 0, 0),
                        value=ueff, hubbard_type='Ueff'
                    )
                )
        hubbard = Hubbard(parameters=hubbard_params)
        return HubbardStructureData.from_structure(afm_structure, hubbard)

    # Baseline: round-3 base + 5A's validated mixing_mode='plain'.
    BASELINE_SYSTEM = {
        'ecutwfc': 50.0, 'ecutrho': 400.0,
        'occupations': 'smearing', 'smearing': 'marzari-vanderbilt', 'degauss': 0.02,
        'lspinorb': False, 'noncolin': True, 'nosym': True,
        'starting_magnetization': {'Ce1': 0.5, 'Ce2': 0.5},
        'angle1': {'Ce1': 54.7356, 'Ce2': 125.2644},
        'angle2': {'Ce1': 45.0, 'Ce2': 225.0},
        'tot_charge': 1.0,
        'starting_charge': {'H': 1.0}
    }
    BASELINE_ELECTRONS = {
        'conv_thr': 1.0e-5,
        'mixing_beta': 0.05,
        'electron_maxstep': 300,
        'diagonalization': 'cg',
        'mixing_mode': 'plain',
        'mixing_ndim': 12,
    }
    DEFAULT_UEFF = 5.0
    DEFAULT_KPOINTS = [1, 1, 1]
    DEFAULT_KPOINTS_OFFSET = [0, 0, 0]

    # (label, system_overrides, electrons_overrides, ueff_override, kpoints_override, resource_override, hypothesis)
    VARIANTS = [
        ('7A_startingwfc_atomic_random', {}, {'startingwfc': 'atomic+random'}, None, None, None, "..."),
        ('7B_magmom_0.9', {'starting_magnetization': {'Ce1': 0.9, 'Ce2': 0.9}}, {}, None, None, None, "..."),
        ('7C_magmom_0.3', {'starting_magnetization': {'Ce1': 0.3, 'Ce2': 0.3}}, {}, None, None, None, "..."),
        ('7D_beta_0.3_plain', {}, {'mixing_beta': 0.3}, None, None, None, "..."),
        ('7E_beta_0.15_plain', {}, {'mixing_beta': 0.15}, None, None, None, "..."),
        ('7F_mixing_mode_TF', {}, {'mixing_mode': 'TF'}, None, None, None, "..."),
        ('7H_diagonalization_rmm', {}, {'diagonalization': 'rmm-davidson'}, None, None, None, "..."),
        ('7I_ndim20_plain', {}, {'mixing_ndim': 20}, None, None, None, "..."),
        ('7J_ueff_4.0', {}, {}, 4.0, None, None, "..."),
        ('7K_ueff_6.0', {}, {}, 6.0, None, None, "..."),
        ('7L_startingwfc_random', {}, {'startingwfc': 'random'}, None, None, None, "..."),
        ('7M_kpoints_2x1x1', {}, {}, None, ([2, 1, 1], [0, 0, 0]), (2, 10 * 3600), "..."),
    ]

    results = []
    for label, sys_override, elec_override, ueff_override, kpoints_override, resource_override, hypothesis in VARIANTS:
        system_params = copy.deepcopy(BASELINE_SYSTEM)
        system_params.update(sys_override)
        electrons_params = copy.deepcopy(BASELINE_ELECTRONS)
        electrons_params.update(elec_override)

        ueff = ueff_override if ueff_override is not None else DEFAULT_UEFF
        variant_structure = build_hubbard_structure(ueff)

        builder = load_code(PW_CODE).get_builder()
        builder.structure = variant_structure
        builder.pseudos = load_group(PSEUDO_GROUP_SSSP).get_pseudos(structure=variant_structure)
        builder.kpoints = DataFactory('core.array.kpoints')()
        if kpoints_override is not None:
            mesh, offset = kpoints_override
            builder.kpoints.set_kpoints_mesh(mesh, offset=offset)
        else:
            builder.kpoints.set_kpoints_mesh(DEFAULT_KPOINTS, offset=DEFAULT_KPOINTS_OFFSET)

        builder.parameters = Dict(dict={
            'CONTROL': {'calculation': 'scf', 'disk_io': 'high'},
            'SYSTEM': system_params,
            'ELECTRONS': electrons_params,
        })

        if resource_override is not None:
            num_machines, max_wallclock_seconds = resource_override
        else:
            num_machines, max_wallclock_seconds = hpc_resources_afm_scf(N)

        builder.settings = Dict(dict={'cmdline': ['-nk', '1']})
        builder.metadata.options = {
            'resources': {'num_machines': num_machines, 'num_mpiprocs_per_machine': 32},
            'max_wallclock_seconds': max_wallclock_seconds,
            'custom_scheduler_commands': '#SBATCH --mem=64G\n',
            'withmpi': True,
        }

        node = submit(builder)
        print(f"VARIANT [{label}] PK: {node.pk}  -- {hypothesis}")
        results.append((label, node.pk))

Round 8 script (scripts/phase5_afm_scf_round8.py)

Combines 7B + 7D + 7K's changes on top of the Round 6 baseline.

import sys
import numpy as np
from aiida.engine import submit
from aiida.orm import load_node, load_code, load_group, Dict, StructureData
from aiida.plugins import DataFactory
from aiida_quantumespresso.common.hubbard import Hubbard, HubbardParameters
from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData
from config import LATTICE_CONST, PW_CODE, PSEUDO_GROUP_SSSP, hpc_resources_afm_scf

if __name__ == '__main__':
    pk = int(sys.argv[1])
    relax_calc = load_node(pk)
    ase_struct = relax_calc.inputs.structure.get_ase()
    ase_struct.set_positions(relax_calc.outputs.output_trajectory.get_positions()[-1])

    a_primitive = LATTICE_CONST
    cell_lengths = ase_struct.get_cell().lengths()
    N = max(1, round(cell_lengths[0] / a_primitive))

    scaled = ase_struct.get_scaled_positions()
    for i, atom in enumerate(ase_struct):
        if atom.symbol == 'Ce':
            x, y, z = scaled[i]
            plane = int(np.round(N * (x + y + z)))
            atom.tag = 1 if plane % 2 == 0 else 2

    afm_structure = StructureData(ase=ase_struct)

    # 7K's win: Ueff 5.0 -> 6.0 eV
    hubbard_params = []
    for i, site in enumerate(afm_structure.sites):
        if site.kind_name.startswith('Ce'):
            hubbard_params.append(
                HubbardParameters(
                    atom_index=i, atom_manifold='4f', neighbour_index=i,
                    neighbour_manifold='4f', translation=(0, 0, 0),
                    value=6.0, hubbard_type='Ueff'
                )
            )
    hubbard = Hubbard(parameters=hubbard_params)
    afm_hubbard_structure = HubbardStructureData.from_structure(afm_structure, hubbard)

    builder = load_code(PW_CODE).get_builder()
    builder.structure = afm_hubbard_structure
    builder.pseudos = load_group(PSEUDO_GROUP_SSSP).get_pseudos(structure=afm_hubbard_structure)
    builder.kpoints = DataFactory('core.array.kpoints')()
    builder.kpoints.set_kpoints_mesh([1, 1, 1])

    builder.parameters = Dict(dict={
        'CONTROL': {'calculation': 'scf', 'disk_io': 'high'},
        'SYSTEM': {
            'ecutwfc': 50.0, 'ecutrho': 400.0,
            'occupations': 'smearing', 'smearing': 'marzari-vanderbilt', 'degauss': 0.02,
            'lspinorb': False, 'noncolin': True, 'nosym': True,
            # 7B's win: 0.5 -> 0.9 µ_B
            'starting_magnetization': {'Ce1': 0.9, 'Ce2': 0.9},
            'angle1': {'Ce1': 54.7356, 'Ce2': 125.2644},
            'angle2': {'Ce1': 45.0, 'Ce2': 225.0},
            'tot_charge': 1.0,
            'starting_charge': {'H': 1.0}
        },
        'ELECTRONS': {
            'conv_thr': 1.0e-5,
            # 7D's win: 0.05 -> 0.3
            'mixing_beta': 0.3,
            'electron_maxstep': 300,
            'diagonalization': 'cg',
            'mixing_mode': 'plain',
            'mixing_ndim': 12,
        }
    })

    num_machines, max_wallclock_seconds = hpc_resources_afm_scf(N)
    builder.settings = Dict(dict={'cmdline': ['-nk', '1']})
    builder.metadata.options = {
        'resources': {'num_machines': num_machines, 'num_mpiprocs_per_machine': 32},
        'max_wallclock_seconds': max_wallclock_seconds,
        'custom_scheduler_commands': '#SBATCH --mem=64G\n',
        'withmpi': True,
    }

    node = submit(builder)
    print(f"ROUND 8 (magmom 0.9 + beta 0.3 + Ueff 6.0, combined): {N}x{N}x{N} AFM SCF Submitted! PK: {node.pk}")

---

Thank you very much for taking the time to look at this, and would appreciate any guidance you can offer.

Logs URL: https://muathhamidi.github.io/muath.hamidi/documents/muSR_Aiida_QE/CeIn3_Aiida_QE_5_Output/CeIn3_Aiida_QE_5_Output.html