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