API reference

DFTK –- The density-functional toolkit. Provides functionality for experimenting with plane-wave density-functional theory algorithms.

The default search location for Pseudopotential data files

TimerOutput object used to store DFTK timings.

Full χ0 application, optionally dropping terms or disabling Sternheimer. All keyword arguments passed to apply_χ0.

Atomic local potential defined by model.atoms.

Nonlocal term coming from norm-conserving pseudopotentials in Kleinmann-Bylander form. $\text{Energy} = \sum_a \sum_{ij} \sum_{n} f_n <ψ_n|p_{ai}> D_{ij} <p_{aj}|ψ_n>.$

We use a simplification of the Resta model DOI 10.1103/physrevb.16.2717 and set $χ_0(q) = \frac{C_0 G^2}{4π (1 - C_0 G^2 / k_{TF}^2)}$ where $C_0 = 1 - ε_r$ with $ε_r$ being the macroscopic relative permittivity. We neglect $K_\text{xc}$, such that $J^{-1} ≈ \frac{k_{TF}^2 - C_0 G^2}{ε_r k_{TF}^2 - C_0 G^2}$

By default it assumes a relative permittivity of 10 (similar to Silicon). εr == 1 is equal to SimpleMixing and εr == Inf to KerkerMixing. The mixing is applied to $ρ$ and $ρ_\text{spin}$ in the same way.

A localised dielectric model for $χ_0$:

\[\sqrt{L(x)} \text{IFFT} \frac{C_0 G^2}{4π (1 - C_0 G^2 / k_{TF}^2)} \text{FFT} \sqrt{L(x)}\]

where $C_0 = 1 - ε_r$, L(r) is a real-space localization function and otherwise the same conventions are used as in DielectricMixing.

Nonlocal "divAgrad" operator $-½ ∇ ⋅ (A ∇)$ where $A$ is a scalar field on the real-space grid. The $-½$ is included, such that this operator is a generalisation of the kinetic energy operator (which is obtained for $A=1$).

Element where the interaction with electrons is modelled as in CohenBergstresser1966. Only the homonuclear lattices of the diamond structure are implemented (i.e. Si, Ge, Sn).

key may be an element symbol (like :Si), an atomic number (e.g. 14) or an element name (e.g. "silicon")

Element interacting with electrons via a bare Coulomb potential (for all-electron calculations) key may be an element symbol (like :Si), an atomic number (e.g. 14) or an element name (e.g. "silicon")

Element interacting with electrons via a pseudopotential model. key may be an element symbol (like :Si), an atomic number (e.g. 14) or an element name (e.g. "silicon")

A simple struct to contain a vector of energies, and utilities to print them in a nice format.

Entropy term -TS, where S is the electronic entropy. Turns the energy E into the free energy F=E-TS. This is in particular useful because the free energy, not the energy, is minimized at self-consistency.

Ewald term: electrostatic energy per unit cell of the array of point charges defined by model.atoms in a uniform background of compensating charge yielding net neutrality.

External potential from the (unnormalized) Fourier coefficients V(G) G is passed in cartesian coordinates

External potential from an analytic function V (in cartesian coordinates). No low-pass filtering is performed.

Fourier space multiplication, like a kinetic energy term: (Hψ)(G) = multiplier(G) ψ(G)

Hartree term: for a decaying potential V the energy would be

1/2 ∫ρ(x)ρ(y)V(x-y) dxdy

with the integral on x in the unit cell and of y in the whole space. For the Coulomb potential with periodic boundary conditions, this is rather

1/2 ∫ρ(x)ρ(y) G(x-y) dx dy

where G is the Green's function of the periodic Laplacian with zero mean (-Δ G = sum{R} 4π δR, integral of G zero on a unit cell).

The same as KerkerMixing, but the Thomas-Fermi wavevector is computed from the current density of states at the Fermi level.

Kerker mixing: $J^{-1} ≈ \frac{|G|^2}{k_{TF}^2 + |G|^2}$ where $k_{TF}$ is the Thomas-Fermi wave vector. For spin-polarized calculations by default the spin density is not preconditioned. Unless a non-default value for $ΔDOS_Ω$ is specified. This value should roughly be the expected difference in density of states (per unit volume) between spin-up and spin-down.

Notes:

• Abinit calls $1/k_{TF}$ the dielectric screening length (parameter dielng)

Kinetic energy: 1/2 sumn fn ∫ |∇ψn|^2.

Discretization information for $k$-point-dependent quantities such as orbitals. More generally, a $k$-point is a block of the Hamiltonian; eg collinear spin is treated by doubling the number of kpoints.

Represents the LDOS-based $χ_0$ model

\[χ_0(r, r') = (-D_\text{loc}(r) δ(r, r') + D_\text{loc}(r) D_\text{loc}(r') / D)\]

where $D_\text{loc}$ is the local density of states and $D$ the density of states. For details see Herbst, Levitt 2020 arXiv:2009.01665

Compute density in real space and its derivatives starting from ρ

Local nonlinearity, with energy ∫f(ρ) where ρ is the density

Magnetic term $A⋅(-i∇)$. It is assumed (but not checked) that $∇⋅A = 0$.

Magnetic field operator A⋅(-i∇).

Model(lattice; n_electrons, atoms, magnetic_moments, terms, temperature,
               smearing, spin_polarization, symmetries)

Creates the physical specification of a model (without any discretization information).

n_electrons is taken from atoms if not specified

spin_polarization is :none by default (paired electrons) unless any of the elements has a non-zero initial magnetic moment. In this case the spin_polarization will be :collinear.

magnetic_moments is only used to determine the symmetry and the spin_polarization; it is not stored inside the datastructure.

smearing is Fermi-Dirac if temperature is non-zero, none otherwise

The symmetries kwarg allows (a) to pass true / false to enable / disable the automatic determination of lattice symmetries or (b) to pass an explicit list of symmetry operations to use for lowering the computational effort. The default behaviour is equal to true, namely that the code checks the specified model in form of the Hamiltonian terms, lattice, atoms and magnetic_moments parameters and from these automatically determines a set of symmetries it can safely use. If you want to pass custom symmetry operations (e.g. a reduced or extended set) use the symmetry_operations function. Notice that this may lead to wrong results if e.g. the external potential breaks some of the passed symmetries. Use false to turn off symmetries completely.

Nonlocal operator in Fourier space in Kleinman-Bylander format, defined by its projectors P matrix and coupling terms D: Hψ = PDP' ψ

Noop operation: don't do anything. Useful for energy terms that don't depend on the orbitals at all (eg nuclei-nuclei interaction).

A plane-wave discretized Model. Normalization conventions:

• Things that are expressed in the G basis are normalized so that if $x$ is the vector, then the actual function is $\sum_G x_G e_G$ with $e_G(x) = e^{iG x} / \sqrt(\Omega)$, where $\Omega$ is the unit cell volume. This is so that, eg $norm(ψ) = 1$ gives the correct normalization. This also holds for the density and the potentials.
• Quantities expressed on the real-space grid are in actual values.

G_to_r and r_to_G convert between these representations.

Creates a PlaneWaveBasis using the kinetic energy cutoff Ecut and a Monkhorst-Pack $k$-point grid. The MP grid can either be specified directly with kgrid providing the number of points in each dimension and kshift the shift (0 or 1/2 in each direction). If not specified a grid is generated using kgrid_from_minimal_spacing with a minimal spacing of 2π * 0.022 per Bohr.

Creates a new basis identical to basis, but with a custom set of kpoints

No preconditioning

(simplified version of) Tetter-Payne-Allan preconditioning ↑ M.P. Teter, M.C. Payne and D.C. Allan, Phys. Rev. B 40, 12255 (1989).

Pseudopotential correction energy. TODO discuss the need for this.

PspHgh(Zion::Number, rloc::Number, cloc::Vector, rp::Vector, h::Vector;
       identifier="", description="")

Construct a Hartwigsen, Goedecker, Teter, Hutter separable dual-space Gaussian pseudopotential (1998). The required parameters are the ionic charge Zion (total charge - valence electrons), the range for the local Gaussian charge distribution rloc, the coefficients for the local part cloc, the projector radius rp (one per AM channel) and the non-local coupling coefficients between the projectors h (one matrix per AM channel).

Linear operators that act on tuples (real, fourier) The main entry point is apply!(out, op, in) which performs the operation out += op*in where out and in are named tuples (real, fourier) They also implement mul! and Matrix(op) for exploratory use.

Real space multiplication by a potential: (Hψ)(r) V(r) ψ(r)

Simple mixing: $J^{-1} ≈ 1$

A term with a constant zero energy.

Exchange-correlation term, defined by a list of functionals and usually evaluated through libxc.

Generic mixing function using a model for the susceptibility composed of the sum of the χ0terms. For valid χ0terms See the subtypes of χ0Model. The dielectric model is solved in real space using a GMRES. Either the full kernel (RPA=false) or only the Hartree kernel (RPA=true) are employed. verbose=true lets the GMRES run in verbose mode (useful for debugging).

CROP-accelerated root-finding iteration for f, starting from x0 and keeping a history of m steps. Optionally warming specifies the number of non-accelerated steps to perform for warming up the history.

In-place version of G_to_r.

G_to_r(basis::PlaneWaveBasis, [kpt::Kpoint, ] f_fourier)

Perform an iFFT to obtain the quantity defined by f_fourier defined on the k-dependent spherical basis set (if kpt is given) or the k-independent cubic (if it is not) on the real-space grid.

G_vectors(basis::PlaneWaveBasis)
G_vectors(basis::PlaneWaveBasis, kpt::Kpoint)

The list of wave vectors $G$ in reduced (integer) coordinates of a basis or a $k$-point kpt.

G_vectors(fft_size::Tuple)

The wave vectors G in reduced (integer) coordinates for a cubic basis set of given sizes.

G_vectors_cart(basis::PlaneWaveBasis)
G_vectors_cart(basis::PlaneWaveBasis, kpt::Kpoint)

The list of $G$ vectors of a given basis or kpt, in cartesian coordinates.

Gplusk_vectors(basis::PlaneWaveBasis, kpt::Kpoint)

The list of $G + k$ vectors, in reduced coordinates.

Gplusk_vectors_cart(basis::PlaneWaveBasis, kpt::Kpoint)

The list of $G + k$ vectors, in cartesian coordinates.

The model for the susceptibility is

\[\begin{aligned} χ_0(r, r') &= (-D_\text{loc}(r) δ(r, r') + D_\text{loc}(r) D_\text{loc}(r') / D) \\ &+ \sqrt{L(x)} \text{IFFT} \frac{C_0 G^2}{4π (1 - C_0 G^2 / k_{TF}^2)} \text{FFT} \sqrt{L(x)} \end{aligned}\]

where $C_0 = 1 - ε_r$, $D_\text{loc}$ is the local density of states, $D$ is the density of states and the same convention for parameters are used as in DielectricMixing. Additionally there is the real-space localization function L(r). For details see Herbst, Levitt 2020 arXiv:2009.01665

Important kwargs passed on to χ0Mixing

• RPA: Is the random-phase approximation used for the kernel (i.e. only Hartree kernel is used and not XC kernel)
• verbose: Run the GMRES in verbose mode.
• reltol: Relative tolerance for GMRES

Increase the temperature used for computing the SCF preconditioners. Initially the temperature is increased by a factor, which is then smoothly lowered towards the temperature used within the model as the SCF converges. Once the density change is below above_ρdiff the mixing temperature is equal to the model temperature.

The model for the susceptibility is

\[\begin{aligned} χ_0(r, r') &= (-D_\text{loc}(r) δ(r, r') + D_\text{loc}(r) D_\text{loc}(r') / D) \end{aligned}\]

where $D_\text{loc}$ is the local density of states, $D$ is the density of states. For details see Herbst, Levitt 2020 arXiv:2009.01665.

Important kwargs passed on to χ0Mixing

• RPA: Is the random-phase approximation used for the kernel (i.e. only Hartree kernel is used and not XC kernel)
• verbose: Run the GMRES in verbose mode.
• reltol: Relative tolerance for GMRES

Accept a step if the energy is at most increasing by max_energy and the residual is at most max_relative_residual times the residual in the previous step.

Flag convergence by using the L2Norm of the change between input density and unpreconditioned output density (ρout)

Flag convergence as soon as total energy change drops below tolerance

Default callback function for self_consistent_field and newton, which prints a convergence table.

Determine the tolerance used for the next diagonalization. This function takes $|ρnext - ρin|$ and multiplies it with ratio_ρdiff to get the next diagtol, ensuring additionally that the returned value is between diagtol_min and diagtol_max and never increases.

Plot the trace of an SCF, i.e. the absolute error of the total energy at each iteration versus the converged energy in a semilog plot. By default a new plot canvas is generated, but an existing one can be passed and reused along with kwargs for the call to plot!. Requires Plots to be loaded.

Adds simplistic checkpointing to a DFTK self-consistent field calculation. Requires JLD2 to be loaded.

apply_K(basis::PlaneWaveBasis, δψ, ψ, ρ, occupation)

Compute the application of K defined at ψ to δψ. ρ is the density issued from ψ. δψ also generates a δρ, computed with compute_δρ.

apply_kernel(basis::PlaneWaveBasis, δρ; kwargs...)

Computes the potential response to a perturbation δρ in real space, as a 4D (i,j,k,σ) array.

Apply a symmetry operation to eigenvectors ψk at a given kpoint to obtain an equivalent point in [-0.5, 0.5)^3 and associated eigenvectors (expressed in the basis of the new $k$-point).

Apply a k-point symmetry operation (the tuple (S, τ)) to a partial density.

apply_Ω(δψ, ψ, H::Hamiltonian, Λ)

Compute the application of Ω defined at ψ to δψ. H is the Hamiltonian computed from ψ and Λ is the set of Rayleigh coefficients ψk' * Hk * ψk at each k-point.

Get the density variation δρ corresponding to a total potential variation δV.

Note: This function assumes that all bands contained in ψ and eigenvalues are sufficiently converged. By default the self_consistent_field routine of DFTK returns 3 extra bands, which are not converged by the eigensolver (see n_ep_extra parameter). These should be discarded from ψ before using this function. However, one can still use additional information from nonconverged bands stored in ψ_extra (typically, the 3 extra bands returned by default in SCF routine).

We use here the full scfres, with a distinction between converged bands and nonconverged extra bands used in the SCF.

Map translating between the atoms datastructure used in DFTK and the indices of the respective atoms in the ASE.Atoms object in pyobj.

Get the lengthscale of the valence density for an atom with n_elec_core core and n_elec_valence valence electrons.

Chemical symbol corresponding to an element

Plan a FFT of type T and size fft_size, spending some time on finding an optimal algorithm. (Inplace, out-of-place) x (forward, backward) FFT plans are returned.

Build form factors (Fourier transforms of projectors) for an atom centered at 0.

Build projection vectors for a atoms array generated by term_nonlocal

Hat = sumij Cij |pi> <pj| Hper = sumR sumij Cij |pi(x-R)> <pj(x-R)| = sumR sum_ij Cij |pi(x-R)> <pj(x-R)|

<ekG'|Hper|ekG> = ... = 1/Ω sumij Cij pihat(k+G') pjhat(k+G)^*

where pihat(q) = ∫_R^3 pi(r) e^{-iqr} dr

We store 1/√Ω pihat(k+G) in proj_vectors.

 bzmesh_ir_wedge(kgrid_size, symmetries; kshift=[0, 0, 0])

Construct the irreducible wedge of a uniform Brillouin zone mesh for sampling $k$-points. The function returns a tuple (kcoords, ksymops), where kcoords are the list of irreducible $k$-points and ksymops are a list of symmetry operations for regenerating the full mesh. symmetries is the tuple returned from symmetry_operations(lattice, atoms, magnetic_moments). tol_symmetry is the tolerance used for searching for symmetry operations.

bzmesh_uniform(kgrid_size; kshift=[0, 0, 0])

Construct a (shifted) uniform Brillouin zone mesh for sampling the $k$-points. The function returns a tuple (kcoords, ksymops), where kcoords are the list of $k$-points and ksymops are a list of symmetry operations (for interface compatibility with PlaneWaveBasis and bzmesh_irreducible. No symmetry reduction is attempted, such that there will be prod(kgrid_size) $k$-points returned and all symmetry operations are the identity.

Return the total ionic charge of an atom type (nuclear charge - core electrons)

Return the total nuclear charge of an atom type

Returns the sum formula of the atoms list as a string.

Compute the matrix $[A_k]_{m,n} = \langle ψ_m^k | g^{\text{per}}_n \rangle$

$g^{per}_n$ are periodized gaussians whose respective centers are given as an (num_bands,1) array [ [center 1], ... ].

Centers are to be given in lattice coordinates and G_vectors in reduced coordinates. The dot product is computed in the Fourier space.

Given an orbital $g_n$, the periodized orbital is defined by : $g^{per}_n = \sum\limits_{R \in {\rm lattice}} g_n( \cdot - R)$. The Fourier coefficient of $g^{per}_n$ at any G is given by the value of the Fourier transform of $g_n$ in G.

Computes the probability (not charge) current, ∑ fn Im(ψn* ∇ψn)

compute_density(basis::PlaneWaveBasis, ψ::AbstractVector, occupation::AbstractVector)

Compute the density for a wave function ψ discretized on the plane-wave grid basis, where the individual k-points are occupied according to occupation. ψ should be one coefficient matrix per $k$-point.

Total density of states at energy ε

Determine the minimal grid size for the cubic basis set to be able to represent product of orbitals (with the default supersampling=2).

Optionally optimize the grid afterwards for the FFT procedure by ensuring factorization into small primes.

The function will determine the smallest parallelepiped containing the wave vectors $|G|^2/2 \leq E_\text{cut} ⋅ \text{supersampling}^2$. For an exact representation of the density resulting from wave functions represented in the spherical basis sets, supersampling should be at least 2.

Compute the forces of an obtained SCF solution. Returns the forces wrt. the fractional lattice vectors. To get cartesian forces use compute_forces_cart. Returns a list of lists of forces [[force for atom in positions] for (element, positions) in atoms] which has the same structure as the atoms object passed to the underlying Model.

Compute the cartesian forces of an obtained SCF solution in Hartree / Bohr. Returns a list of lists of forces [[force for atom in positions] for (element, positions) in atoms] which has the same structure as the atoms object passed to the underlying Model.

Compute the inverse of the lattice. Takes special care of 1D or 2D cases.

compute_kernel(basis::PlaneWaveBasis; kwargs...)

Computes a matrix representation of the full response kernel (derivative of potential with respect to density) in real space. For non-spin-polarized calculations the matrix dimension is prod(basis.fft_size) × prod(basis.fft_size) and for collinear spin-polarized cases it is 2prod(basis.fft_size) × 2prod(basis.fft_size). In this case the matrix has effectively 4 blocks

\[\left(\begin{array}{cc} K_{αα} & K_{αβ}\\ K_{βα} & K_{ββ} \end{array}\right)\]

Local density of states, in real space

compute_nos(ε, basis, eigenvalues; smearing=basis.model.smearing,
            temperature=basis.model.temperature)

The number of Kohn-Sham states in a temperature window of width temperature around the energy ε contributing to the DOS at temperature T.

This quantity is not a physical quantity, but rather a dimensionless approximate measure for how well properties near the Fermi surface are sampled with the passed smearing and temperature T. It increases with both T and better sampling of the BZ with $k$-points. A value $\gg 1$ indicates a good sampling of properties near the Fermi surface.

Find the occupation and Fermi level.

Find Fermi level and occupation for the given parameters, assuming a band gap and zero temperature. This function is for DEBUG purposes only, and the finite-temperature version with 0 temperature should be preferred.

Compute the reciprocal lattice. We use the convention that the reciprocal lattice is the set of G vectors such that G ⋅ R ∈ 2π ℤ for all R in the lattice.

Compute the stresses (= 1/Vol dE/d(M*lattice), taken at M=I) of an obtained SCF solution.

Return a sparse matrix that maps quantities given on basis_in and kpt_in to quantities on basis_out and kpt_out.

Return a list of sparse matrices (one per $k$-point) that map quantities given in the basis_in basis to quantities given in the basis_out basis.

Compute unit cell volume volume. In case of 1D or 2D case, the volume is the length/surface.

Compute the independent-particle susceptibility. Will blow up for large systems. For non-spin-polarized calculations the matrix dimension is prod(basis.fft_size) × prod(basis.fft_size) and for collinear spin-polarized cases it is 2prod(basis.fft_size) × 2prod(basis.fft_size). In this case the matrix has effectively 4 blocks, which are:

\[\left(\begin{array}{cc} (χ_0)_{αα} & (χ_0)_{αβ} \\ (χ_0)_{βα} & (χ_0)_{ββ} \end{array}\right)\]

Return the data directory with pseudopotential files

:none if no element has a magnetic moment, else :collinear or :full

Default logic to determine the symmetry operations to be used in the model.

Default random Gaussian guess for maximally-localised wannier functions generated in reduced coordinates.

Function for diagonalising each $k$-Point blow of ham one step at a time. Some logic for interpolating between $k$-points is used if interpolate_kpoints is true and if no guesses are given. eigensolver is the iterative eigensolver that really does the work, operating on a single $k$-Block. eigensolver should support the API eigensolver(A, X0; prec, tol, maxiter) prec_type should be a function that returns a preconditioner when called as prec(ham, kpt)

Compute the diameter of the unit cell

Computes the ground state by direct minimization. kwargs... are passed to Optim.Options(). Note that the resulting ψ are not necessarily eigenvectors of the Hamiltonian.

Convenience function to disable all threading in DFTK and assert that Julia threading is off as well.

Compute divergence of an operand function, which returns the cartesian x,y,z components in real space when called with the arguments 1 to 3. The divergence is also returned as a real-space array.

Compute the electrostatic interaction energy per unit cell between point charges in a uniform background of compensating charge to yield net neutrality. the lattice and recip_lattice should contain the lattice and reciprocal lattice vectors as columns. charges and positions are the point charges and their positions (as an array of arrays) in fractional coordinates. If forces is not nothing, minus the derivatives of the energy with respect to positions is computed.

Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)

del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)

Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)

Sarmiento-Perez, Silvana, Marques 2015 (DOI 10.1021/acs.jctc.5b00529)

Perdew, Burke, Ernzerhof (DOI: 10.1103/PhysRevLett.77.3865)

Xu, Goddard 2004 (DOI 10.1063/1.1771632)

Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)

del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)

Zhang, Yang 1998 (DOI 10.1103/physrevlett.80.890)

Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)

Sarmiento-Perez, Silvana, Marques 2015 (DOI 10.1021/acs.jctc.5b00529)

Perdew, Burke, Ernzerhof (DOI: 10.1103/PhysRevLett.77.3865)

Xu, Goddard 2004 (DOI 10.1063/1.1771632)

Perdew, Wang correlation from 1992 (10.1103/PhysRevB.45.13244)

VWN5 LDA correlation according to Vosko, Wilk, and Nusair, (DOI 10.1139/p80-159).

LDA Slater exchange (DOI: 10.1017/S0305004100016108 and 10.1007/BF01340281)

Compute the correction term for properly modelling the interaction of the pseudopotential core with the compensating background charge induced by the Ewald term.

eval_psp_energy_correction([T=Float64,] psp, n_electrons)

Evaluate the energy correction to the Ewald electrostatic interaction energy of one unit cell, which is required compared the Ewald expression for point-like nuclei. n_electrons is the number of electrons per unit cell. This defines the uniform compensating background charge, which is assumed here.

Notice: The returned result is the energy per unit cell and not the energy per volume. To obtain the latter, the caller needs to divide by the unit cell volume.

eval_psp_local_fourier(psp, q)

Evaluate the local part of the pseudopotential in reciprocal space: V(q) = ∫R^3 Vloc(r) e^{-iqr} dr = 4π ∫{R+} sin(qr)/q r e^{-iqr} dr

eval_psp_local_real(psp, r)

Evaluate the local part of the pseudopotential in real space.

eval_psp_projector_fourier(psp, i, l, q)

Evaluate the radial part of the i-th projector for angular momentum l at the reciprocal vector with modulus q: p(q) = ∫R^3 p{il}(r) e^{-iqr} dr = 4π ∫{R+} r^2 p{il}(r) j_l(q r) dr

eval_psp_projector_real(psp, i, l, r)

Evaluate the radial part of the i-th projector for angular momentum l in real-space at the vector with modulus r.

Find the Fermi level.

Maximal occupation of a state (2 for non-spin-polarized electrons, 1 otherwise).

Implements a primitive search to find an irreducible subset of kpoints amongst the provided kpoints.

Gather the distributed data of a quantity depending on k-Points on the master process and return it. On the other (non-master) processes nothing is returned.

Gather the distributed $k$-point data on the master process and return it as a PlaneWaveBasis. On the other (non-master) processes nothing is returned. The returned object should not be used for computations and only to extract data for post-processing and serialisation to disk.

Build a superposition of Gaussians as a guess for the density and magnetisation. Expects a list of tuples (coefficient, length, position) for each of the Gaussian, which follow the functional form

\[\hat{ρ}(G) = \text{coefficient} \exp\left(-(2π \text{length} |G|)^2\right)\]

and are placed at position (in fractional coordinates).

guess_density(basis, magnetic_moments)

Build a superposition of atomic densities (SAD) guess density.

We take for the guess density a gaussian centered around the atom, of length specified by atom_decay_length, normalized to get the right number of electrons

\[\hat{ρ}(G) = Z \exp\left(-(2π \text{length} |G|)^2\right)\]

When magnetic moments are provided, construct a symmetry-broken density guess. The magnetic moments should be specified in units of $μ_B$.

Returns a new Hamiltonian with local potential replaced by the given one

Extract the high-symmetry $k$-point path corresponding to the passed model using Brillouin.jl. Uses the conventions described in the reference work by Cracknell, Davies, Miller, and Love (CDML). Of note, this has minor differences to the $k$-path reference (Y. Himuma et. al. Comput. Mater. Sci. 128, 140 (2017)) underlying the path-choices of Brillouin.jl, specifically for oA and mC Bravais types. The kline_density is given in number of $k$-points per inverse bohrs (i.e. overall in units of length).

Issues a warning in case the passed lattice does not match the expected primitive.

Return the index tuple I such that G_vectors(basis)[I] == G or the index i such that G_vectors(basis, kpoint)[i] == G. Returns nothing if outside the range of valid wave vectors.

Interpolate a function expressed in a basis basis_in to a basis basis_out This interpolation uses a very basic real-space algorithm, and makes a DWIM-y attempt to take into account the fact that basisout can be a supercell of basisin

Interpolate some data from one $k$-point to another. The interpolation is fast, but not necessarily exact or even normalized. Intended only to construct guesses for iterative solvers

is_metal(band_data, εF, tol)

Determine whether the provided bands indicate the material is a metal, i.e. where bands are cut by the Fermi level.

Selects a kgrid size which ensures that at least a n_kpoints total number of $k$-points are used. The distribution of $k$-points amongst coordinate directions is as uniformly as possible, trying to achieve an identical minimal spacing in all directions.

Selects a kgrid size to ensure a minimal spacing (in inverse Bohrs) between kpoints. A reasonable spacing is 0.13 inverse Bohrs (around $2π * 0.04 \AA^{-1}$).

Construct the coordinates of the $k$-points in a (shifted) Monkorst-Pack grid

Return the index range of $k$-points that have a particular spin component.

list_psp(element; functional, family, core, datadir_psp)

List the pseudopotential files known to DFTK. Allows various ways to restrict the displayed files.

Examples

julia> list_psp(family="hgh")

will list all HGH-type pseudopotentials and

julia> list_psp(family="hgh", functional="lda")

will only list those for LDA (also known as Pade in this context) and

julia> list_psp(:O, core=:semicore)

will list all oxygen semicore pseudopotentials known to DFTK.

Return a DFTK-compatible atoms object loaded from an ASE Atoms, a pymatgen Structure or a compatible file (e.g. xyz file, cif file etc.)

Load a DFTK-compatible atoms representation from a supported pymatgen object. All atoms are using a Coulomb model.

Return a DFTK-compatible lattice object loaded from an ASE Atoms, a pymatgen Structure or a compatible file (e.g. xyz file, cif file etc.)

Load a DFTK-compatible lattice object from a supported python object (e.g. pymatgen or ASE)

Load a pseudopotential file from the library of pseudopotentials. The file is searched in the directory datadir_psp() and by the key. If the key is a path to a valid file, the extension is used to determine the type of the pseudopotential file format and a respective class is returned.

load_scfres(filename)

Load back an scfres, which has previously been stored with save_scfres. Note the warning in save_scfres.

Radial local potential, in Fourier space: V(q) = int_{R^3} V(x) e^{-iqx} dx.

Radial local potential, in real space.

Build a DFT model from the specified atoms, with the specified functionals.

Build an LDA model (Teter93 parametrization) from the specified atoms.

Build an PBE-GGA model from the specified atoms. DOI:10.1103/PhysRevLett.77.3865

Build a SCAN meta-GGA model from the specified atoms. DOI:10.1103/PhysRevLett.115.036402

Convenience constructor, which builds a standard atomic (kinetic + atomic potential) model. Use extra_terms to add additional terms.

Number of processors used in MPI. Can be called without ensuring initialization.

Return the number of core electrons

Return the number of valence electrons

newton(basis::PlaneWaveBasis{T}; ψ0=nothing,
       tol=1e-6, tol_cg=1e-10, maxiter=20, verbose=false,
       callback=NewtonDefaultCallback(),
       is_converged=NewtonConvergenceDensity(tol))

Newton algorithm. Be careful that the starting point needs to be not too far from the solution.

Obtain new density ρ by diagonalizing ham.

Bring $k$-point coordinates into the range [-0.5, 0.5)

Computes the matrix $[M^{k,b}]_{m,n} = \langle u_{m,k} | u_{n,k+b} \rangle$ for given k, kpb = k+b.

G_shift is the "shifting" vector, correction due to the periodicity conditions imposed on k -> ψk. It is non zero if kpb is taken in another unit cell of the reciprocal lattice. We use here that : ``u{n(k + Gshift)}(r) = e^{-i*\langle Gshift,r \rangle} u_{nk}``

Compute and plot the band structure. n_bands selects the number of bands to compute. If this value is absent and an scfres is used to start the calculation a default of n_bands_scf + 5sqrt(n_bands_scf) is used. The unit used to plot the bands can be selected using the unit parameter. Like in the rest of DFTK Hartree is used by default. Another standard choices is unit=u"eV" (electron volts). The kline_density is given in number of $k$-points per inverse bohrs (i.e. overall in units of length).

Plot the density of states over a reasonable range

The local potential of a HGH pseudopotentials in reciprocal space can be brought to the form $Q(t) / (t^2 exp(t^2 / 2))$ where $t = r_\text{loc} q$ and Q is a polynomial of at most degree 8. This function returns Q.

The nonlocal projectors of a HGH pseudopotentials in reciprocal space can be brought to the form $Q(t) exp(-t^2 / 2)$ where $t = r_l q$ and Q is a polynomial. This function returns Q.

Estimate an upper bound for the argument q after which abs(eval_psp_local_fourier(psp, q)) is a strictly decreasing function.

Estimate an upper bound for the argument q after which eval_psp_projector_fourier(psp, q) is a strictly decreasing function.

In-place version of r_to_G!. NOTE: If kpt is given, not only f_fourier but also f_real is overwritten.

r_to_G(basis::PlaneWaveBasis, [kpt::Kpoint, ] f_real)

Perform an FFT to obtain the Fourier representation of f_real. If kpt is given, the coefficients are truncated to the k-dependent spherical basis set.

r_vectors(basis::PlaneWaveBasis)

The list of $r$ vectors, in reduced coordinates. By convention, this is in [0,1)^3.

r_vectors_cart(basis::PlaneWaveBasis)

The list of $r$ vectors, in cartesian coordinates.

Generate a physically valid random density integrating to the given number of electrons.

Read the .nnkp file provided by the preprocessing routine of Wannier90 (i.e. "wannier90.x -pp prefix") Returns:

1. the array 'nnkpts' of k points, their respective nearest neighbors and associated shifing vectors (non zero if the neighbor is located in another cell).
2. the number 'nntot' of neighbors per k point.

TODO: add the possibility to exclude bands

Wannerize the obtained bands using wannier90. By default all converged bands from the scfres are employed (change with n_bands kwargs) and n_wannier = n_bands wannier functions are computed using random Gaussians as guesses. All keyword arguments supported by Wannier90 for the disentanglement may be added as keyword arguments. The function returns the fileprefix.

save_scfres(filename, scfres)

Save an scfres obtained from self_consistent_field to a file. The format is determined from the file extension. Currently the following file extensions are recognized and supported:

• jld2: A JLD2 file. Stores the complete state and can be used (with load_scfres) to restart an SCF from a checkpoint or post-process an SCF solution. See Saving SCF results on disk and SCF checkpoints for details.
• vts: A VTK file for visualisation e.g. in paraview. Stores the density, spin density and some metadata (energy, Fermi level, occupation etc.). Supports these keyword arguments:
  • save_ψ: Save the real-space representation of the orbitals as well (may lead to larger files).
  • extra_data: Dict{String,Array} with additional data on the 3D real-space grid to store into the VTK file.

Create a simple anderson-accelerated SCF solver. m specifies the number of steps to keep the history of.

Use the two iteration states info and info_next to find a damping value from a quadratic model for the SCF energy. Returns nothing if the constructed model is not considered trustworthy, else returns the suggested damping.

Create a damped SCF solver updating the density as x = β * x_new + (1 - β) * x

Create a NLSolve-based SCF solver, by default using an Anderson-accelerated fixed-point scheme, keeping m steps for Anderson acceleration. See the NLSolve documentation for details about the other parameters and methods.

Function to select a subset of eigenpairs on each $k$-Point. Works on the Tuple returned by diagonalize_all_kblocks.

Solve the Kohn-Sham equations with a SCF algorithm, starting at ρ.

solve_ΩplusK(basis::PlaneWaveBasis{T}, ψ, res, occupation;
             tol_cg=1e-10, verbose=false) where T

Return δψ where (Ω+K) δψ = rhs

Convert the DFTK atoms datastructure into a tuple of datastructures for use with spglib. positions contains positions per atom, numbers contains the mapping atom to a unique number for each indistinguishable element, spins contains the $z$-component of the initial magnetic moment on each atom, mapping contains the mapping of the numbers to the element objects in DFTK and collinear whether the atoms mark a case of collinear spin or not. Notice that if collinear is false then spins is garbage.

Returns crystallographic conventional cell according to the International Table of Crystallography Vol A (ITA) in case primitive=false. If primitive=true the primitive lattice is returned in the convention of the reference work of Cracknell, Davies, Miller, and Love (CDML). Of note this has minor differences to the primitive setting choice made in the ITA.

Explicit spin components of the KS orbitals and the density

Split an iterable evenly into N chunks, which will be returned.

Apply various standardisations to a lattice and a list of atoms. It uses spglib to detect symmetries (within tol_symmetry), then cleans up the lattice according to the symmetries (unless correct_symmetry is false) and returns the resulting standard lattice and atoms. If primitive is true (default) the primitive unit cell is returned, else the conventional unit cell is returned.

Find the subset of symmetries compatible with the grid induced by the given kcoords and ksymops

Filter out the symmetry operations that respect the symmetries of the discrete BZ grid

Symmetrize the forces in reduced coordinates, forces given as an array forces[iel][α,i]

Symmetrize the stress tensor, given as a Matrix in cartesian coordinates

Symmetrize a density by applying all the basis (by default) symmetries and forming the average.

Return the $k$-point symmetry operations associated to a lattice and atoms.

Get the total local potential of the given Hamiltonian, in real space in the spin components.

Transfer Bloch wave between two basis sets. Limited feature set.

Transfer an array ψk defined on basisin $k$-point kptin to basisout $k$-point kptout.

Compute the index mapping between two bases. Returns two arrays idcs_in and idcs_out such that ψkout[idcs_out] = ψkin[idcs_in] does the transfer from ψkin (defined on basis_in and kpt_in) to ψkout (defined on basis_out and kpt_out).

" Convert a basis into one that doesn't use BZ symmetry. This is mainly useful for debug purposes (e.g. in cases we don't want to bother thinking about symmetries).

Sum an array over kpoints, taking weights into account

Write the eigenvalues in a format readable by Wannier90.

Write a win file at the indicated prefix. Parameters to Wannier90 can be added as kwargs: e.g. num_iter=500.

Returns the (l,m) real spherical harmonic Ylm(r). Consistent with https://en.wikipedia.org/wiki/Tableofsphericalharmonics#Realsphericalharmonics

Shortened version of the @timeit macro from TimerOutputs, which writes to the DFTK timer.

A term in the Hermite delta expansion

Standard Hermite function using physicist's convention.

Entropy. Note that this is a function of the energy x, not of occupation(x). This function satisfies s' = x f' (see https://www.vasp.at/vasp-workshop/k-points.pdf p. 12 and https://arxiv.org/pdf/1805.07144.pdf p. 18.

Occupation at x, where in practice x = (ε - εF) / T.

Derivative of the occupation function, approximation to minus the delta function.

(f(x) - f(y))/(x - y), computed stably in the case where x and y are close