Computing U with hp.x

Goal

In E11 the U of 4.6 eV was "a given". This time, compute the U by linear response (DFPT): no empirical parameters, and a value that belongs to this system, this pseudopotential, and this projector. Background in Chapter 14.

pw.x (scf with U ≈ 0, conv_thr 1e-12) → hp.x → FeO.Hubbard_parameters.dat

New cards and variables

Input files

feo_hp_scf.in · feo.hp.in · run.sh

The scf input is the cell of E10 with two changes:

&ELECTRONS
  conv_thr         = 1.0d-12    ! very tight: the response matrix inherits this quality
  mixing_beta      = 0.2
  mixing_mode      = 'local-TF'
  electron_maxstep = 300
/
...
! effectively-zero U: registers the Fe-3d manifolds with hp.x without
! biasing the ground state
HUBBARD (ortho-atomic)
U Fe1-3d 1.0d-8
U Fe2-3d 1.0d-8

The hp.x input in full:

! E12 step 2: hp.x computes the response matrices chi0 and chi by DFPT
! and returns U = chi0^-1 - chi^-1 per Hubbard atom.
&INPUTHP
  prefix       = 'FeO'       ! must match the scf of step 1
  outdir       = './tmp/'
  nq1 = 2, nq2 = 2, nq3 = 2  ! q-point grid for the perturbation; CONVERGE THIS
                             ! (a 1x1x1 U is not trustworthy)
  conv_thr_chi = 1.0d-6      ! chi convergence; 1.0d-8 is unreachable on a metallic
                             ! (GGA) ground state, noise floor ~1e-7 (measured)
  iverbosity   = 2           ! print per-iteration chi values (worth watching)
/

A measured note on conv_thr_chi: we first ran with 1.0d-8, and even after the χ values had stabilized to seven digits, the residual kept bouncing around a noise floor near 10⁻⁷ (checked out to 46 iterations) without ever crossing the threshold. For systems whose GGA ground state is metallic, like FeO, that is where the numerical noise of the response function lives, so we relaxed the threshold to 1.0d-6. The relative scatter of χ is below 0.1%, which moves U by less than 0.01 eV.

Run

#!/bin/bash
# tightly converged scf first, then the hp.x linear-response run
set -e
pw.x -in feo_hp_scf.in > feo_hp_scf.out   # ~3 min on 8 ranks (measured)
hp.x -in feo.hp.in     > feo.hp.out       # ~2 h on 8 ranks, 4 irreducible q (measured)
echo "--- computed Hubbard parameters ---"
cat FeO.Hubbard_parameters.dat

hp.x runs one perturbation series per inequivalent Hubbard atom, so it takes a while (-nk pools apply, and q-points can be split with start_q/last_q).

What to check: measured

FeO.Hubbard_parameters.dat lists the per-atom U along with the full χ₀ and χ matrices. The computed U = 5.22 eV sits 0.6 eV above the conventional 4.6 eV used in E11. Neither number is "the right one" in isolation; the pair illustrates that U only has meaning together with its projector, pseudopotential, and magnetic order. The principle: for your own system, use the U computed for it.

Exercises

1. Change nq from 1×1×1 to 2×2×2 and watch U shift. A U without q-grid convergence cannot be trusted.
2. Put the computed U into the HUBBARD card of E11 and see how the moments and DOS respond.
3. Iterate scf → hp.x once more with the new U (self-consistent U) and check that the value stabilizes.

Related chapters

14 Computing U with hp.x · 13 DFT+U and the HUBBARD card