How to run a DFT+U calculation

Plain density functional theory (DFT) calculations based on the local density approximation (LDA) or the generalised gradient approximation (GGA) usually fail to describe strongly correlated system properly. LDA and GGA often incorrectly predict a metallic ground state for systems like FeO, CoO or UO2. A possible remedy for this deficit is the addition of a Hubbard U correction which is known as DFT+U (or also LDA+U or GGA+U) approximation. Here we show how DFT+U calculation can be performed with CP2K.

EOS calculation for cubic FeO using DFT+U

The following CP2K input file shows how an equation of state (EOS) calculation for cubic FeO (space group #225, wüstite) using 2x2x2 unit cells to account for the antiferromagnetic ordering along [111], can be performed with CP2K using DFT+U:

# FeO (#225, Fm-3m, wustite, antiferromagnetic along [111])
 
# f094 = 0.94^(1/3) .. f120 = 1.20^(1/3)
@SET f094      0.9795861087155615099
@SET f096      0.9864848297321880954
@SET f098      0.9932883883792685831
@SET f100      1.0000000000000000000
@SET f102      1.0066227095601130159
@SET f104      1.0131594038201774399
@SET f106      1.0196128224222165137
@SET f108      1.0259855680060181936
@SET f110      1.0322801154563671592
@SET f112      1.0384988203702208052
@SET f120      1.0626585691826110661
 
@SET run_type  energy_force
@SET system    FeO
@SET project   ${system}-2x2x2
 
@SET a         8.56*${f100}
@SET b         ${a}
@SET c         ${a}
 
@SET a_ref     ${a}*${f120}
@SET b_ref     ${a_ref}
@SET c_ref     ${a_ref}
 
@SET PLUS_U    on
@SET U_Fe      1.9
 
&GLOBAL
  PRINT_LEVEL low
  PROJECT ${project}
  RUN_TYPE ${run_type}
  WALLTIME 1800
&END GLOBAL
 
&FORCE_EVAL
  METHOD Quickstep
  &DFT
    LSD
    PLUS_U_METHOD Mulliken
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    &MGRID
      CUTOFF 600.0
      REL_CUTOFF 60.0
    &END MGRID
    &QS
      EPS_DEFAULT 1.0E-12
    &END QS
    &SCF
      EPS_SCF 3.0E-7
      MAX_SCF 21
      SCF_GUESS restart
      &OT on
        MINIMIZER CG
        PRECONDITIONER FULL_SINGLE_INVERSE
        STEPSIZE 0.1
      &END OT
      &OUTER_SCF on
        EPS_SCF 3.0E-7
        MAX_SCF 20
      &END OUTER_SCF
      &PRINT
        &RESTART
          BACKUP_COPIES 0
          &EACH
            QS_SCF 10
          &END EACH
        &END RESTART
      &END PRINT
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC ${a} ${b} ${c}
      &CELL_REF
        ABC ${a_ref} ${b_ref} ${c_ref}
      &END CELL_REF
    &END CELL
    &COORD
      SCALED
      Fe_a                0                   0                   0
      Fe_b                0                 1/4                 1/4
      Fe_b              1/4                   0                 1/4
      Fe_b              1/4                 1/4                   0
      Fe_b                0                   0                 1/2
      Fe_b              1/2                   0                   0
      Fe_b                0                 1/2                   0
      Fe_a                0                 1/4                 3/4
      Fe_a              1/2                 1/4                 1/4
      Fe_a                0                 3/4                 1/4
      Fe_a              1/4                   0                 3/4
      Fe_a              3/4                   0                 1/4
      Fe_a              1/4                 1/2                 1/4
      Fe_a              1/4                 1/4                 1/2
      Fe_a              3/4                 1/4                   0
      Fe_a              1/4                 3/4                   0
      Fe_a                0                 1/2                 1/2
      Fe_a              1/2                   0                 1/2
      Fe_a              1/2                 1/2                   0
      Fe_b                0                 3/4                 3/4
      Fe_b              1/2                 1/4                 3/4
      Fe_b              1/2                 3/4                 1/4
      Fe_b              1/4                 1/2                 3/4
      Fe_b              3/4                   0                 3/4
      Fe_b              3/4                 1/2                 1/4
      Fe_b              1/4                 3/4                 1/2
      Fe_b              3/4                 1/4                 1/2
      Fe_b              3/4                 3/4                   0
      Fe_b              1/2                 1/2                 1/2
      Fe_a              1/2                 3/4                 3/4
      Fe_a              3/4                 1/2                 3/4
      Fe_a              3/4                 3/4                 1/2
      O                 1/4                 1/4                 1/4
      O                 3/4                 3/4                 3/4
      O                 3/4                 3/4                 1/4
      O                 1/4                 1/4                 3/4
      O                 3/4                 1/4                 3/4
      O                 1/4                 3/4                 1/4
      O                 1/4                 3/4                 3/4
      O                 3/4                 1/4                 1/4
      O                 1/4                 1/2                 1/2
      O                 3/4                   0                   0
      O                 3/4                   0                 1/2
      O                 1/4                 1/2                   0
      O                 3/4                 1/2                   0
      O                 1/4                   0                 1/2
      O                 1/4                   0                   0
      O                 3/4                 1/2                 1/2
      O                 1/2                 1/4                 1/2
      O                   0                 3/4                   0
      O                   0                 3/4                 1/2
      O                 1/2                 1/4                   0
      O                   0                 1/4                   0
      O                 1/2                 3/4                 1/2
      O                 1/2                 3/4                   0
      O                   0                 1/4                 1/2
      O                 1/2                 1/2                 1/4
      O                   0                   0                 3/4
      O                   0                   0                 1/4
      O                 1/2                 1/2                 3/4
      O                   0                 1/2                 3/4
      O                 1/2                   0                 1/4
      O                 1/2                   0                 3/4
      O                   0                 1/2                 1/4
    &END COORD
    &KIND O
      BASIS_SET DZVP-MOLOPT-SR-GTH-q6
      POTENTIAL GTH-PBE-q6
      # O(2-): 2s2 2p6
      &BS on
        &ALPHA
          N    2
          L    1
          NEL  2
        &END ALPHA
        &BETA
          N    2
          L    1
          NEL  2
        &END BETA
      &END BS
    &END KIND
    &KIND Fe_a
      BASIS_SET DZVP-MOLOPT-SR-GTH-q16
      POTENTIAL GTH-PBE-q16
      # Fe(2+): 3d6 (alpha, spin up)
      &BS on
        &ALPHA
          N    4   3
          L    0   2
          NEL -2   4
        &END ALPHA
        &BETA
          N    4   3
          L    0   2
          NEL -2  -4
        &END BETA
      &END BS
      &DFT_PLUS_U on
        L 2
        U_MINUS_J [eV] ${U_Fe}
      &END DFT_PLUS_U
    &END KIND
    &KIND Fe_b
      BASIS_SET DZVP-MOLOPT-SR-GTH-q16
      POTENTIAL GTH-PBE-q16
      # Fe(2+); 3d6 (beta, spin down)
      &BS on
        &ALPHA
          N    4   3
          L    0   2
          NEL -2  -4
        &END ALPHA
        &BETA
          N    4   3
          L    0   2
          NEL -2   4
        &END BETA
      &END BS
      &DFT_PLUS_U on
        L 2
        U_MINUS_J [eV] ${U_Fe}
      &END DFT_PLUS_U
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

The suggested Hubbard Ueff value of 1.9 eV for iron using CP2K has been retrieved from the work of Kéri et al., ES&T 51, 10585-10594 (2017).

The calculations for different cell sizes, e.g. for 98% of the experimental cell volume, can be accomplished by changing the scaling factor f100 in this line

@SET a         8.56*${f100}

to

@SET a         8.56*${f098}

and likewise for other cell sizes. It recommended to start with a cell size close to the assumed minimum and to employ the wavefunction restart files consecutively for the larger and smaller cell sizes to avoid a convergence to different states for compressed and enlarged cells.

Running the example above while applying the scaling factors f094 to f112 using the CP2K version 2023.2 resulted in the following total energies:

# a [Angstrom]  Total energy [Hartree]
  8.385257      -4466.80963727
  8.444310      -4466.85496007
  8.502549      -4466.88640521
  8.560000      -4466.90536056
  8.616690      -4466.91306948
  8.672644      -4466.91064762
  8.727886      -4466.89906783
  8.782436      -4466.87923722
  8.836318      -4466.85199239
  8.889550      -4466.81806793