calculateFiniteDifferenceLinearResponseHubbardU

calculateFiniteDifferenceLinearResponseHubbardU(configuration, shell_indices, hubbard_term, supercell_repetitions, potential_displacement=PhysicalQuantity(0.02, eV))

Calculate Hubbard Us for a selection of shells using the finite difference linear response method.

Parameters:

• configuration (BulkConfiguration) – The configuration to calculate the Us for.

• shell_indices (dict) – A map between elements and the list of shell indices to calculate the Us for.

• hubbard_term (Onsite | Dual) – The projection method to use for calculating shell occupations.

• supercell_repetitions (tuple of three integers) – The supercell size.

• potential_displacement (PhysicalQuantity of type energy) – The finite difference potential displacement to use.

Returns:

A list of calculated Us for each atom in the configuration.

Return type:

list

Notes

The Hubbard energy contribution is given by

\[E_\mathrm{Hubbard-U} = \sum_{a,i} \frac{U_{a,i}}{2} \sum_{m,\sigma} \left\{ [n_{a,i}]_{mm}^{\sigma\sigma} - \sum_{m',\sigma'} [n_{a,i}]_{mm'}^{\sigma\sigma'} [n_{a,i}]_{m'm}^{\sigma\sigma'} \right\},\]

where \([n_{a,i}]_{mm'}^{\sigma\sigma'} = \sum_\mu f_\mu \langle \psi_\mu | [\hat{P}_{a,i}]_{mm'}^{\sigma\sigma'} |\psi_\mu \rangle\) with \([\hat{P}_{a,i}]_{mm'}^{\sigma\sigma'}\) the specific projection operator for a shell \(i\) on an atom \(a\).

This function calculates the Hubbard U’s, \(U_{a,i}\), for selected atoms and shells in a system using linear response theory based on constrained DFT[1]. In practice this is done by introducing a term in the Hamiltonian given by

\[\hat{V}_\mathrm{lr-U} = \sum_{a,i} \alpha_{a,i} \hat{P}_{a,i},\]

where \(\hat{P}_{a,i} = \sum_{mm'} \sum_{\sigma\sigma'} [\hat{P}_{a,i}]_{mm'}^{\sigma\sigma'}\) is the projection operator for a given atom + shell (Dual or Onsite) and \(\alpha_{a,i}\) is an arbitrary potential displacement. One then calculates the shell occupation response function corresponding to this potential displacement

\[\chi_{ai;bi} = \frac{\partial n_{a,i}}{\partial \alpha_{b,i}} \approx \frac{n_{a,i}(\Delta\alpha_{b,i}) - n_{a,i}(-\Delta\alpha_{b,i})}{2\Delta\alpha_{b,i}}\]

where \(n_{a,i}(\Delta\alpha_{b,i}) = \sum_m n^{a,i}_{mm}(\Delta\alpha_{b,i})\) is the occupation of shell \(i\) on atom \(a\) calculated for a self-consistent DFT calculation with a small projection potential displacement on shell \(i\) on atom \(b\) of size \(\Delta\alpha_{b,i}\). The Hubbard U can then be shown to be given by

\[U_{a,i} = (\chi_0^{-1} - \chi^{-1})_{ai,ai}\]

where \(\chi_0\) is the response matrix as given above, but where the Hamiltonian is not updated, in practice found by calculating the shell occupations at the first SCF iteration.

For systems with periodic boundary conditions the pertubation potential will get artificically repeated. In theory one should apply the perturbative potential to a single atom/shell atom a time for the entire infinite system. This cannot of course be done in practice with this method, so instead one in practive performs the perturbative calculations in a supercell of a size big enough that the spatial effect of the perturbation has decayed so as to not overlap with that of the neighboring supercells. One will often have to perform calculations of increasing supercell size to check for convergence.

Note

For semiconductors and insulators one should use an occupation function smearing of practically zero. This is because the desired shell may be part of the valence band and when introducing the finite perturbation, this may lead to isolated states above/below the valence band for the two finite perturbations. These states will get a perturbation dependent occupation due to the non-zero smearing, and this artificial occupation difference can become a significant contribution to the shell occupation response function.

Note

It can also be required to decrease the SCF convergence threshold as the small finite perturbation may only lead to small changes in the density and Hamiltonian. As the function calculates the finite difference derivatives of the shell occupations it is important that the shell occupations are converged within a tolerance that is mcuh smaller than the finite difference.

Usage Examples

The following script calculates the U values for anti-ferromagnetic NiO using the Dual projection and a supercell of size \(3 \times 3 \times 3\), prints out the values and saves them to a “pickle” file for later use.

# Set up anti-ferromagnetic NiO
configuration = BulkConfiguration(
    bravais_lattice=Rhombohedral(5.10522 * Angstrom, 34.12513143396482 * Degrees),
    elements=[Nickel, Nickel, Oxygen, Oxygen],
    fractional_coordinates=[
        [ 0.125,  0.125,  0.125],
        [ 0.625,  0.625,  0.625],
        [ 0.375,  0.375,  0.375],
        [ 0.875,  0.875,  0.875]
    ]
)

calculator = LCAOCalculator(
    exchange_correlation=SGGA.PBE,
    basis_set=BasisGGAPseudoDojo.Medium,
    numerical_accuracy_parameters=NumericalAccuracyParameters(
        k_point_sampling=KpointDensity(2.0 * Angstrom),
        occupation_method=FermiDirac(1e-6 * eV) # We need a low smearing to avoid artificial occupations.
    ),
    iteration_control_parameters=IterationControlParameters(
        tolerance=1e-5, # We need a low tolerance to ensure convergence with a small perturbation
        non_convergence_behavior=StopCalculation(),
    ),
    checkpoint_handler=NoCheckpointHandler
)

# Set calculator and antiferromagnetic spin ordering.
configuration.setCalculator(
    calculator,
    initial_spin=InitialSpin(scaled_spins=[1.0, -1.0, 0.0, 0.0]))

# Calculate the standard DFT ground state and save for later use.
configuration.update()
nlsave('NiO_PBE.hdf5', configuration)

# Calculate the Hubbard U's for the 3d shell.
fdlr_us_per_atom = calculateFiniteDifferenceLinearResponseHubbardU(
    configuration=configuration,
    shell_indices={Nickel: [3]}, # Calculate U for the 3d shell of Nickel = 3rd basis orbital
    hubbard_term=Dual, # Dual projection is most local
    supercell_repetitions=(1, 1, 1),
    potential_displacement=0.1 * eV
)

for atom_index, u_values_per_shell in fdlr_us_per_atom:
    nlprint('{:d}: U for 3d shell: {:.3f} eV'.format(atom_index, u_values_per_shell[0].inUnitsOf(eV)))

# Save the calculated U's in a "pickle" file for later use.
import pickle
with open('NiO_fdlr_us_per_atom.pckl', 'wb') as fd:
    pickle.dump(fdlr_us_per_atom, fd)

The following script loads the calculated U values from the pickle file and uses the function createLocalDFTUBasisSet() to perfom a density of states calculation based on DFT+U with the Dual projection and the calculates U values.

# Get the configuration with PBE ground state
pbe_configuration = nlread('NiO_PBE.hdf5', BulkConfiguration)[0]

# Perform reference DOS for PBE
pbe_dos = DensityOfStates(pbe_configuration)
nlsave('NiO_dos.hdf5', pbe_dos, object_id='NiO DOS PBE')

# Get the previously calculated U values.
import pickle
with open('NiO_fdlr_us_per_atom.pckl', 'rb') as fd:
    fdlr_u_per_atom = pickle.load(fd)

pbe_calculator = pbe_configuration.calculator()

# Create a new configuration and basis set with the calculated U values.
configuration_with_tags, basis_set_with_u = createLocalDFTUBasisSet(
    configuration=pbe_configuration,
    basis_set=pbe_calculator.basisSet(),
    shell_indices={Nickel: [3]},
    atom_u_values=fdlr_u_per_atom
)

# Set up spin-polarized PBE with Hubbard U using Dual projection.
xc_with_dftu = ExchangeCorrelation(
    exchange=PerdewBurkeErnzerhofExchange,
    correlation=PerdewBurkeErnzerhofCorrelation,
    hubbard_term=Dual,
    number_of_spins=2
)

# Create a copy of the original calculator but change xc to DFT+U and use the basis set
# with U values.
calculator_with_u = pbe_calculator(
    exchange_correlation=xc_with_dftu,
    basis_set=basis_set_with_u)

configuration_with_tags.setCalculator(
    calculator_with_u,
    initial_spin=InitialSpin(scaled_spins=[1.0, -1.0, 0.0, 0.0]))

# Perform the DFT+U calculation
configuration_with_tags.update()

# Calculate density of states
density_of_states_with_u = DensityOfStates(configuration=configuration_with_tags)
nlsave('NiO_dos.hdf5', density_of_states_with_u, object_id='NiO DOS PBE+U')

[1]

M. Cococcioni and S. de Gironcoli. Linear response approach to the calculation of the effective interaction parameters in the LDA+U method. Phys. Rev. B, 71:035105, Jan 2005. doi:10.1103/PhysRevB.71.035105.