/
search_index.js
3 lines (3 loc) · 195 KB
/
search_index.js
1
2
3
var documenterSearchIndex = {"docs":
[{"location":"publications/#Publications","page":"Publications","title":"Publications","text":"","category":"section"},{"location":"publications/","page":"Publications","title":"Publications","text":"This page contains a list of publications dealing with DFTK or which make use of DFTK. Feel free to drop us a line if you want your work to be added here. To cite DFTK please use (Image: DOI).","category":"page"},{"location":"publications/","page":"Publications","title":"Publications","text":"M. F. Herbst and A. Levitt. Black-box inhomogeneous preconditioning for self-consistent field iterations in density functional theory (2020). ArXiv:2009.01665 (Supplementary material and computational scripts).\nM. F. Herbst, A. Levitt and E. Cancès. A posteriori error estimation for the non-self-consistent Kohn-Sham equations. Faraday discussions, in press (2020). ArXiv:2004.13549. (Reference implementation).\nE. Cancès, G. Kemlin and A. Levitt. Convergence analysis of direct minimization and self-consistent iterations (2020). ArXiv:2004.09088.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/collinear_magnetism.jl\"","category":"page"},{"location":"examples/collinear_magnetism/#Collinear-spin-and-magnetic-systems","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"","category":"section"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"In this example we consider iron in the BCC phase. To show that this material is ferromagnetic we will model it once allowing collinear spin polarization and once without and compare the resulting SCF energies. In particular the ground state can only be found if collinear spins are allowed.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"First we setup BCC iron without spin polarization using a single iron atom inside the unit cell.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"using DFTK\n\na = 5.42352 # Bohr\nlattice = a / 2 * [[-1 1 1];\n [ 1 -1 1];\n [ 1 1 -1]]\nFe = ElementPsp(:Fe, psp=load_psp(\"hgh/lda/Fe-q8.hgh\"))\natoms = [Fe => [zeros(3)]];\nnothing #hide","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"To get the ground-state energy we use an LDA model and rather moderate discretisation parameters.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"kgrid = [3, 3, 3] # k-point grid (Regular Monkhorst-Pack grid)\nEcut = 15 # kinetic energy cutoff in Hartree\nmodel_nospin = model_LDA(lattice, atoms, temperature=0.01)\nbasis_nospin = PlaneWaveBasis(model_nospin, Ecut; kgrid=kgrid)\n\nscfres_nospin = self_consistent_field(basis_nospin, tol=1e-6, mixing=KerkerMixing());\nnothing #hide","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"scfres_nospin.energies","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"Since we did not specify any initial magnetic moment on the iron atom, DFTK will automatically assume that a calculation with only spin-paired electrons should be performed. As a result the obtained ground state features no spin-polarization.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"Now we repeat the calculation, but give the iron atom an initial magnetic moment. For specifying the magnetic moment pass the desired excess of spin-up over spin-down electrons at each centre to the Model and the guess density functions. In this case we seek the state with as many spin-parallel d-electrons as possible. In our pseudopotential model the 8 valence electrons are 2 pair of s-electrons, 1 pair of d-electrons and 4 unpaired d-electrons giving a desired magnetic moment of 4 at the iron centre. The structure (i.e. pair mapping and order) of the magnetic_moments array needs to agree with the atoms array and 0 magnetic moments need to be specified as well.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"magnetic_moments = [Fe => [4, ]];\nnothing #hide","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"tip: Units of the magnetisation and magnetic moments in DFTK\nUnlike all other quantities magnetisation and magnetic moments in DFTK are given in units of the Bohr magneton μ_B, which in atomic units has the value frac12. Since μ_B is (roughly) the magnetic moment of a single electron the advantage is that one can directly think of these quantities as the excess of spin-up electrons or spin-up electron density.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"We repeat the calculation using the same model as before. DFTK now detects the non-zero moment and switches to a collinear calculation.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"model = model_LDA(lattice, atoms, magnetic_moments=magnetic_moments, temperature=0.01)\nbasis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\nρspin = guess_spin_density(basis, magnetic_moments)\nscfres = self_consistent_field(basis, tol=1e-6, ρspin=ρspin, mixing=KerkerMixing());\nnothing #hide","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"scfres.energies","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"note: Model and magnetic moments\nDFTK does not store the magnetic_moments inside the Model, but only uses them to determine the lattice symmetries. This step was taken to keep Model (which contains the physical model) independent of the details of the numerical details such as the initial guess for the spin density.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"In direct comparison we notice the first, spin-paired calculation to be a little higher in energy","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"println(\"No magnetization: \", scfres_nospin.energies.total)\nprintln(\"Magnetic case: \", scfres.energies.total)\nprintln(\"Difference: \", scfres.energies.total - scfres_nospin.energies.total);\nnothing #hide","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"Notice that with the small cutoffs we use to generate the online documentation the calculation is far from converged. With more realistic parameters a larger energy difference of about 0.1 Hartree is obtained.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"The spin polarization in the magnetic case is visible if we consider the occupation of the spin-up and spin-down Kohn-Sham orbitals. Especially for the d-orbitals these differ rather drastically. For example for the first k-point:","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"iup = 1\nidown = iup + length(scfres.basis.kpoints) ÷ 2\n@show scfres.occupation[iup][1:7]\n@show scfres.occupation[idown][1:7];\nnothing #hide","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"Similarly the eigenvalues differ","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"@show scfres.eigenvalues[iup][1:7]\n@show scfres.eigenvalues[idown][1:7];\nnothing #hide","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"note: k-points in collinear calculations\nFor collinear calculations the kpoints field of the PlaneWaveBasis object contains each k-point coordinate twice, once associated with spin-up and once with down-down. The list first contains all spin-up k-points and then all spin-down k-points, such that iup and idown index the same k-point, but differing spins.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"We can observe the spin-polarization by looking at the density of states (DOS) around the Fermi level, where the spin-up and spin-down DOS differ.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"using Plots\nplot_dos(scfres)","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"Similarly the band structure shows clear differences between both spin components.","category":"page"},{"location":"examples/collinear_magnetism/","page":"Collinear spin and magnetic systems","title":"Collinear spin and magnetic systems","text":"plot_bandstructure(scfres, kline_density=3, unit=:eV)","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/custom_solvers.jl\"","category":"page"},{"location":"examples/custom_solvers/#Custom-solvers","page":"Custom solvers","title":"Custom solvers","text":"","category":"section"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"In this example, we show how to define custom solvers. Our system will again be silicon, because we are not very imaginative","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"using DFTK, LinearAlgebra\n\na = 10.26\nlattice = a / 2 * [[0 1 1.];\n [1 0 1.];\n [1 1 0.]]\nSi = ElementPsp(:Si, psp=load_psp(\"hgh/lda/Si-q4\"))\natoms = [Si => [ones(3)/8, -ones(3)/8]]\n\n# We take very (very) crude parameters\nmodel = model_LDA(lattice, atoms)\nkgrid = [1, 1, 1]\nEcut = 5\nbasis = PlaneWaveBasis(model, Ecut; kgrid=kgrid);\nnothing #hide","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"We define our custom fix-point solver: simply a damped fixed-point","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"function my_fp_solver(f, x0, max_iter; tol)\n mixing_factor = .7\n x = x0\n fx = f(x)\n for n = 1:max_iter\n inc = fx - x\n if norm(inc) < tol\n break\n end\n x = x + mixing_factor * inc\n fx = f(x)\n end\n (fixpoint=x, converged=norm(fx-x) < tol)\nend;\nnothing #hide","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"Our eigenvalue solver just forms the dense matrix and diagonalizes it explicitly (this only works for very small systems)","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"function my_eig_solver(A, X0; maxiter, tol, kwargs...)\n n = size(X0, 2)\n A = Array(A)\n E = eigen(A)\n λ = E.values[1:n]\n X = E.vectors[:, 1:n]\n (λ=λ, X=X, residual_norms=[], iterations=0, converged=true, n_matvec=0)\nend;\nnothing #hide","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"Finally we also define our custom mixing scheme. It will be a mixture of simple mixing (for the first 2 steps) and than default to Kerker mixing. In the mixing interface δF is (ρ_textout - ρ_textin), i.e. the difference in density between two subsequent SCF steps and the mix function returns δρ, which is added to ρ_textin to yield ρ_textnext, the density for the next SCF step.","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"struct MyMixing\n α # Damping parameter\nend\nMyMixing() = MyMixing(0.7)\n\nfunction DFTK.mix(mixing::MyMixing, basis, δF::RealFourierArray, δF_spin=nothing; n_iter, kwargs...)\n if n_iter <= 2\n # Just do simple mixing on total density and spin density (if it exists)\n (mixing.α * δF, isnothing(δF_spin) ? nothing : mixing.α * δF_spin)\n else\n # Use the KerkerMixing from DFTK\n DFTK.mix(KerkerMixing(α=mixing.α), basis, δF, δF_spin; kwargs...)\n end\nend","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"That's it! Now we just run the SCF with these solvers","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"scfres = self_consistent_field(basis;\n tol=1e-8,\n solver=my_fp_solver,\n eigensolver=my_eig_solver,\n mixing=MyMixing());\nnothing #hide","category":"page"},{"location":"examples/custom_solvers/","page":"Custom solvers","title":"Custom solvers","text":"Note that the default convergence criterion is on the difference of energy from one step to the other; when this gets below tol, the \"driver\" self_consistent_field artificially makes the fixpoint solver think it's converged by forcing f(x) = x. You can customize this with the is_converged keyword argument to self_consistent_field.","category":"page"},{"location":"advanced/useful_formulas/#Useful-formulas","page":"Useful formulas","title":"Useful formulas","text":"","category":"section"},{"location":"advanced/useful_formulas/","page":"Useful formulas","title":"Useful formulas","text":"This section holds a collection of formulae, which are helpful when working with DFTK and plane-wave DFT in general. See also Notation and conventions for a description of the conventions used in the equations.","category":"page"},{"location":"advanced/useful_formulas/#Fourier-transforms","page":"Useful formulas","title":"Fourier transforms","text":"","category":"section"},{"location":"advanced/useful_formulas/","page":"Useful formulas","title":"Useful formulas","text":"The Fourier transform is\nwidehatf(q) = int_mathbb R^3 e^-i q cdot x f(x) dx\nFourier transforms of centered functions: If f(x) = R(x) Y_l^m(xx), then\nbeginaligned\n hat f( q)\n = int_mathbb R^3 R(x) Y_l^m(xx) e^-i q cdot x dx \n = sum_l = 0^infty 4 pi i^l\n sum_m = -l^l int_mathbb R^3\n R(x) j_l(q x)Y_l^m(-qq) Y_l^m(xx)\n Y_l^mast(xx)\n dx \n = 4 pi Y_l^m(-qq) i^l\n int_mathbb R^+ r^2 R(r) j_l(q r) dr\n endaligned\nThis also holds true for real spherical harmonics.","category":"page"},{"location":"advanced/useful_formulas/#Spherical-harmonics","page":"Useful formulas","title":"Spherical harmonics","text":"","category":"section"},{"location":"advanced/useful_formulas/","page":"Useful formulas","title":"Useful formulas","text":"Plane wave expansion formula\ne^i q cdot r =\n 4 pi sum_l = 0^infty sum_m = -l^l\n i^l j_l(q r) Y_l^m(qq) Y_l^mast(rr)\nSpherical harmonics orthogonality\nint_mathbbS^2 Y_l^m*(r)Y_l^m(r) dr\n = delta_ll delta_mm\nThis also holds true for real spherical harmonics.","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/metallic_systems.jl\"","category":"page"},{"location":"examples/metallic_systems/#Temperature-and-metallic-systems","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"","category":"section"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"In this example we consider the modeling of a magnesium lattice as a simple example for a metallic system. For our treatment we will use the PBE exchange-correlation functional. First we import required packages and setup the lattice. Again notice that DFTK uses the convention that lattice vectors are specified column by column.","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"using DFTK\nusing Plots\n\na = 3.01794 # bohr\nb = 5.22722 # bohr\nc = 9.77362 # bohr\nlattice = [[-a -a 0]; [-b b 0]; [0 0 -c]]\nMg = ElementPsp(:Mg, psp=load_psp(\"hgh/pbe/Mg-q2\"))\natoms = [Mg => [[2/3, 1/3, 1/4], [1/3, 2/3, 3/4]]];\nnothing #hide","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"Next we build the PBE model and discretize it. Since magnesium is a metal we apply a small smearing temperature to ease convergence using the Fermi-Dirac smearing scheme. Note that both the Ecut is too small as well as the minimal k-point spacing kspacing far too large to give a converged result. These have been selected to obtain a fast execution time. By default PlaneWaveBasis chooses a kspacing of 2π * 0.022 inverse Bohrs, which is much more reasonable.","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"kspacing = 0.5 # Minimal spacing of k-points,\n# in units of wavevectors (inverse Bohrs)\nEcut = 5 # kinetic energy cutoff in Hartree\ntemperature = 0.01 # Smearing temperature in Hartree\n\nmodel = model_DFT(lattice, atoms, [:gga_x_pbe, :gga_c_pbe];\n temperature=temperature,\n smearing=DFTK.Smearing.FermiDirac())\nkgrid = kgrid_size_from_minimal_spacing(lattice, kspacing)\nbasis = PlaneWaveBasis(model, Ecut, kgrid=kgrid);\nnothing #hide","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"Finally we run the SCF. Two magnesium atoms in our pseudopotential model result in four valence electrons being explicitly treated. Nevertheless this SCF will solve for eight bands by default in order to capture partial occupations beyond the Fermi level due to the employed smearing scheme.","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"scfres = self_consistent_field(basis);\nnothing #hide","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"scfres.occupation[1]","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"scfres.energies","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"The fact that magnesium is a metal is confirmed by plotting the density of states around the Fermi level.","category":"page"},{"location":"examples/metallic_systems/","page":"Temperature and metallic systems","title":"Temperature and metallic systems","text":"plot_dos(scfres)","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/scf_callbacks.jl\"","category":"page"},{"location":"examples/scf_callbacks/#Monitoring-self-consistent-field-calculations","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"","category":"section"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"The self_consistent_field function takes as the callback keyword argument one function to be called after each iteration. This function gets passed the complete internal state of the SCF solver and can thus be used both to monitor and debug the iterations as well as to quickly patch it with additional functionality.","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"This example discusses a few aspects of the callback function taking again our favourite silicon example.","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"We setup silicon in an LDA model using the ASE interface to build the silicon lattice, see Creating slabs with ASE for more details.","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"using DFTK\nusing PyCall\n\nsilicon = pyimport(\"ase.build\").bulk(\"Si\")\natoms = load_atoms(silicon)\natoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional=\"lda\")) => position\n for (el, position) in atoms]\nlattice = load_lattice(silicon);\n\nmodel = model_LDA(lattice, atoms)\nkgrid = [3, 3, 3] # k-point grid\nEcut = 5 # kinetic energy cutoff in Hartree\nbasis = PlaneWaveBasis(model, Ecut; kgrid=kgrid);\nnothing #hide","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"DFTK already defines a few callback functions for standard tasks. One example is the usual convergence table, which is defined in the callback ScfDefaultCallback. Another example is ScfPlotTrace, which records the total energy at each iteration and uses it to plot the convergence of the SCF graphically once it is converged. For details and other callbacks see src/scf/scf_callbacks.jl.","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"note: Callbacks are not exported\nCallbacks are not exported from the DFTK namespace as of now, so you will need to use them, e.g., as DFTK.ScfDefaultCallback and DFTK.ScfPlotTrace.","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"In this example we define a custom callback, which plots the change in density at each SCF iteration after the SCF has finished. For this we first define the empty plot canvas and an empty container for all the density differences:","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"using Plots\np = plot(yaxis=:log)\ndensity_differences = Float64[];\nnothing #hide","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"The callback function itself gets passed a named tuple similar to the one returned by self_consistent_field, which contains the input and output density of the SCF step as ρin and ρout. Since the callback gets called both during the SCF iterations as well as after convergence just before self_consistent_field finishes we can both collect the data and initiate the plotting in one function.","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"using LinearAlgebra\n\nfunction plot_callback(info)\n if info.stage == :finalize\n plot!(p, density_differences, label=\"|ρout - ρin|\", markershape=:x)\n else\n push!(density_differences, norm(info.ρout.real - info.ρin.real))\n end\n info\nend\ncallback = DFTK.ScfDefaultCallback() ∘ plot_callback;\nnothing #hide","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"Notice that for constructing the callback function we chained the plot_callback (which does the plotting) with the ScfDefaultCallback, such that when using the plot_callback function with self_consistent_field we still get the usual convergence table printed. We run the SCF with this callback ...","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"scfres = self_consistent_field(basis, tol=1e-8, callback=callback);\nnothing #hide","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"... and show the plot","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"p","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"The info object passed to the callback contains not just the densities but also the complete Bloch wave (in ψ), the occupation, band eigenvalues and so on. See src/scf/self_consistent_field.jl for all currently available keys.","category":"page"},{"location":"examples/scf_callbacks/","page":"Monitoring self-consistent field calculations","title":"Monitoring self-consistent field calculations","text":"tip: Debugging with callbacks\nVery handy for debugging SCF algorithms is to employ callbacks with an @infiltrate from Infiltrator.jl to interactively monitor what is happening each SCF step.","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/pymatgen.jl\"","category":"page"},{"location":"examples/pymatgen/#Creating-supercells-with-pymatgen","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"","category":"section"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"The Pymatgen python library allows to setup solid-state calculations using a flexible set of classes as well as an API to an online data base of structures. Its Structure and Lattice objects are directly supported by the DFTK load_atoms and load_lattice functions, such that DFTK may be readily used to run calculation on systems defined in pymatgen. Using the pymatgen_structure function a conversion from DFTK to pymatgen structures is also possible. In the following we use this to create a silicon supercell and find its LDA ground state using direct minimisation.","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"First we setup the silicon lattice in DFTK.","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"using DFTK\n\na = 10.263141334305942 # Lattice constant in Bohr\nlattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]\nSi = ElementPsp(:Si, psp=load_psp(\"hgh/lda/Si-q4\"))\natoms = [Si => [ones(3)/8, -ones(3)/8]];\nnothing #hide","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"Next we make a [2, 2, 2] supercell using pymatgen","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"pystruct = pymatgen_structure(lattice, atoms)\npystruct.make_supercell([2, 2, 2])\nlattice = load_lattice(pystruct)\natoms = [Si => [s.frac_coords for s in pystruct.sites]];\nnothing #hide","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"Setup an LDA model and discretize using a single kpoint and a small Ecut of 5 Hartree.","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"model = model_LDA(lattice, atoms)\nbasis = PlaneWaveBasis(model, 5, kgrid=(1, 1, 1))","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"Find the ground state using direct minimisation (always using SCF is boring ...)","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"scfres = direct_minimization(basis, tol=1e-5);\nnothing #hide","category":"page"},{"location":"examples/pymatgen/","page":"Creating supercells with pymatgen","title":"Creating supercells with pymatgen","text":"scfres.energies","category":"page"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/arbitrary_floattype.jl\"","category":"page"},{"location":"examples/arbitrary_floattype/#Arbitrary-floating-point-types","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"","category":"section"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"Since DFTK is completely generic in the floating-point type in its routines, there is no reason to perform the computation using double-precision arithmetic (i.e.Float64). Other floating-point types such as Float32 (single precision) are readily supported as well. On top of that we already reported[HLC2020] calculations in DFTK using elevated precision from DoubleFloats.jl or interval arithmetic using IntervalArithmetic.jl. In this example, however, we will concentrate on single-precision computations with Float32. The setup of such a reduced-precision calculation is basically identical to the regular case, since Julia automatically compiles all routines of DFTK at the precision, which is used for the lattice vectors. Apart from setting up the model with an explicit cast of the lattice vectors to Float32, there is thus no change in user code required:","category":"page"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"[HLC2020]: M. F. Herbst, A. Levitt, E. Cancès. A posteriori error estimation for the non-self-consistent Kohn-Sham equations ArXiv 2004.13549","category":"page"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"using DFTK\n\n# Setup silicon lattice\na = 10.263141334305942 # lattice constant in Bohr\nlattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]\nSi = ElementPsp(:Si, psp=load_psp(:Si, functional=\"lda\"))\natoms = [Si => [ones(3)/8, -ones(3)/8]]\n\n# Cast to Float32, setup model and basis\nmodel = model_DFT(Array{Float32}(lattice), atoms, [:lda_x, :lda_c_vwn])\nEcut = 7\nbasis = PlaneWaveBasis(model, Ecut, kgrid=[4, 4, 4])\n\n# Run the SCF\nscfres = self_consistent_field(basis, tol=1e-4);\nnothing #hide","category":"page"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"To check the calculation has really run in Float32, we check the energies and density are expressed in this floating-point type:","category":"page"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"scfres.energies","category":"page"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"eltype(scfres.energies.total)","category":"page"},{"location":"examples/arbitrary_floattype/","page":"Arbitrary floating-point types","title":"Arbitrary floating-point types","text":"eltype(scfres.ρ.real)","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/gross_pitaevskii.jl\"","category":"page"},{"location":"examples/gross_pitaevskii/#Gross-Pitaevskii-equation-in-one-dimension","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"","category":"section"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"In this example we will use DFTK to solve the Gross-Pitaevskii equation, and use this opportunity to explore a few internals.","category":"page"},{"location":"examples/gross_pitaevskii/#The-model","page":"Gross-Pitaevskii equation in one dimension","title":"The model","text":"","category":"section"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"The Gross-Pitaevskii equation (GPE) is a simple non-linear equation used to model bosonic systems in a mean-field approach. Denoting by ψ the effective one-particle bosonic wave function, the time-independent GPE reads in atomic units:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":" H ψ = left(-frac12 Δ + V + 2 C ψ^2right) ψ = μ ψ qquad ψ_L^2 = 1","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"where C provides the strength of the boson-boson coupling. It's in particular a favorite model of applied mathematicians because it has a structure simpler than but similar to that of DFT, and displays interesting behavior (especially in higher dimensions with magnetic fields, see Gross-Pitaevskii equation with magnetism).","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"We wish to model this equation in 1D using DFTK. First we set up the lattice. For a 1D case we supply two zero lattice vectors,","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"a = 10\nlattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"which is special cased in DFTK to support 1D models.","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"For the potential term V we just pick a harmonic potential. The real-space grid is in 01) in fractional coordinates( see Lattices and lattice vectors), therefore:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"pot(x) = (x - a/2)^2;\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"We setup each energy term in sequence: kinetic, potential and nonlinear term. For the non-linearity we use the PowerNonlinearity(C, α) term of DFTK. This object introduces an energy term C ρ(r)^α dr to the total energy functional, thus a potential term α C ρ^α-1. In our case we thus need the parameters","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"C = 1.0\nα = 2;\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"... and with this build the model","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"using DFTK\nusing LinearAlgebra\n\nn_electrons = 1 # Increase this for fun\nterms = [Kinetic(),\n ExternalFromReal(r -> pot(r[1])),\n PowerNonlinearity(C, α),\n]\nmodel = Model(lattice; n_electrons=n_electrons, terms=terms,\n spin_polarization=:spinless); # use \"spinless electrons\"\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"We discretize using a moderate Ecut (For 1D values up to 5000 are completely fine) and run a direct minimization algorithm:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"Ecut = 500\nbasis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))\nscfres = direct_minimization(basis, tol=1e-8) # This is a constrained preconditioned LBFGS\nscfres.energies","category":"page"},{"location":"examples/gross_pitaevskii/#Internals","page":"Gross-Pitaevskii equation in one dimension","title":"Internals","text":"","category":"section"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"We use the opportunity to explore some of DFTK internals.","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"Extract the converged density and the obtained wave function:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"ρ = real(scfres.ρ.real)[:, 1, 1] # converged density\nψ_fourier = scfres.ψ[1][:, 1]; # first kpoint, all G components, first eigenvector\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"Transform the wave function to real space and fix the phase:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\nψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"Check whether ψ is normalised:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"x = a * vec(first.(DFTK.r_vectors(basis)))\nN = length(x)\ndx = a / N # real-space grid spacing\n@assert sum(abs2.(ψ)) * dx ≈ 1.0","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"The density is simply built from ψ:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"norm(scfres.ρ.real - abs2.(ψ))","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"We summarize the ground state in a nice plot:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"using Plots\n\np = plot(x, real.(ψ), label=\"real(ψ)\")\nplot!(p, x, imag.(ψ), label=\"imag(ψ)\")\nplot!(p, x, ρ, label=\"ρ\")","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"The energy_hamiltonian function can be used to get the energy and effective Hamiltonian (derivative of the energy with respect to the density matrix) of a particular state (ψ, occupation). The density ρ associated to this state is precomputed and passed to the routine as an optimization.","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"E, ham = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)\n@assert E.total == scfres.energies.total","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"Now the Hamiltonian contains all the blocks corresponding to kpoints. Here, we just have one kpoint:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"H = ham.blocks[1];\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"H can be used as a linear operator (efficiently using FFTs), or converted to a dense matrix:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"ψ11 = scfres.ψ[1][:, 1] # first kpoint, first eigenvector\nHmat = Array(H) # This is now just a plain Julia matrix,\n# which we can compute and store in this simple 1D example\n@assert norm(Hmat * ψ11 - H * ψ11) < 1e-10","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"Let's check that ψ11 is indeed an eigenstate:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"norm(H * ψ11 - dot(ψ11, H * ψ11) * ψ11)","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"Build a finite-differences version of the GPE operator H, as a sanity check:","category":"page"},{"location":"examples/gross_pitaevskii/","page":"Gross-Pitaevskii equation in one dimension","title":"Gross-Pitaevskii equation in one dimension","text":"A = Array(Tridiagonal(-ones(N - 1), 2ones(N), -ones(N - 1)))\nA[1, end] = A[end, 1] = -1\nK = A / dx^2 / 2\nV = Diagonal(pot.(x) + C .* α .* (ρ.^(α-1)))\nH_findiff = K + V;\nmaximum(abs.(H_findiff*ψ - (dot(ψ, H_findiff*ψ) / dot(ψ, ψ)) * ψ))","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/gross_pitaevskii_2D.jl\"","category":"page"},{"location":"examples/gross_pitaevskii_2D/#Gross-Pitaevskii-equation-with-magnetism","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"","category":"section"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"We solve the 2D Gross-Pitaevskii equation with a magnetic field. This is similar to the previous example (Gross-Pitaevskii equation in one dimension), but with an extra term for the magnetic field. We reproduce here the results of https://arxiv.org/pdf/1611.02045.pdf Fig. 10","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"using DFTK\nusing StaticArrays\nusing Plots","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"Unit cell. Having one of the lattice vectors as zero means a 2D system","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"a = 15\nlattice = a .* [[1 0 0.]; [0 1 0]; [0 0 0]];\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"Confining scalar potential, and magnetic vector potential","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"pot(x, y, z) = ((x - a/2)^2 + (y - a/2)^2)/2\nω = .6\nApot(x, y, z) = ω * @SVector [y - a/2, -(x - a/2), 0]\nApot(X) = Apot(X...);\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"Parameters","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"Ecut = 20 # Increase this for production\nη = 500\nC = η/2\nα = 2\nn_electrons = 1; # Increase this for fun\nnothing #hide","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"Collect all the terms, build and run the model","category":"page"},{"location":"examples/gross_pitaevskii_2D/","page":"Gross-Pitaevskii equation with magnetism","title":"Gross-Pitaevskii equation with magnetism","text":"terms = [Kinetic(),\n ExternalFromReal(X -> pot(X...)),\n PowerNonlinearity(C, α),\n Magnetic(Apot),\n]\nmodel = Model(lattice; n_electrons=n_electrons,\n terms=terms, spin_polarization=:spinless) # \"spinless electrons\"\nbasis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))\nscfres = direct_minimization(basis, tol=1e-5) # Reduce tol for production\nheatmap(scfres.ρ.real[:, :, 1], c=:blues)","category":"page"},{"location":"guide/installation/#Installation","page":"Installation","title":"Installation","text":"","category":"section"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"Similar to installing any other registered package in Julia just run from a Julia REPL:","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"] add DFTK","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"or if you like to be fully up to date:","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"] add DFTK#master","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"At least Julia 1.4 is required.","category":"page"},{"location":"guide/installation/#Python-dependencies","page":"Installation","title":"Python dependencies","text":"","category":"section"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"Some parts of the code require a working Python installation with the pymatgen module. Check out which version of python is used by the PyCall.jl package. You can do this for example with the Julia commands","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"using PyCall\nPyCall.python","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"Then use the corresponding package manager (usually apt, pip, pip3 or conda) to install aforementioned libraries, for example","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"pip install pymatgen","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"or","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"conda install -c conda-forge pymatgen","category":"page"},{"location":"guide/installation/","page":"Installation","title":"Installation","text":"Afterwards you're all set and should be able to run the code in the examples directory.","category":"page"},{"location":"api/#API-reference","page":"API reference","title":"API reference","text":"","category":"section"},{"location":"api/","page":"API reference","title":"API reference","text":"This page provides a plain list of all documented functions, structs, modules and macros in DFTK. Note that this list is neither structured, complete nor particularly clean, so it only provides rough orientation at the moment. The best reference is the code itself.","category":"page"},{"location":"api/","page":"API reference","title":"API reference","text":"Modules = [DFTK, DFTK.Smearing]","category":"page"},{"location":"api/#DFTK.DFTK","page":"API reference","title":"DFTK.DFTK","text":"DFTK –- The density-functional toolkit. Provides functionality for experimenting with plane-wave density-functional theory algorithms.\n\n\n\n\n\n","category":"module"},{"location":"api/#DFTK.DFTK_DATADIR","page":"API reference","title":"DFTK.DFTK_DATADIR","text":"The default search location for Pseudopotential data files\n\n\n\n\n\n","category":"constant"},{"location":"api/#DFTK.timer","page":"API reference","title":"DFTK.timer","text":"TimerOutput object used to store DFTK timings.\n\n\n\n\n\n","category":"constant"},{"location":"api/#DFTK.Applyχ0Model","page":"API reference","title":"DFTK.Applyχ0Model","text":"Full χ0 application, optionally dropping terms or disabling Sternheimer. All keyword arguments passed to apply_χ0.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.AtomicLocal","page":"API reference","title":"DFTK.AtomicLocal","text":"Atomic local potential defined by model.atoms.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.AtomicNonlocal","page":"API reference","title":"DFTK.AtomicNonlocal","text":"Nonlocal term coming from norm-conserving pseudopotentials in Kleinmann-Bylander form. textEnergy = sum_a sum_ij sum_n f_n ψ_np_ai D_ij p_ajψ_n\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.DielectricMixing","page":"API reference","title":"DFTK.DielectricMixing","text":"We use a simplification of the Resta model DOI 10.1103/physrevb.16.2717 and set χ_0(q) = fracC_0 G^24π (1 - C_0 G^2 k_TF^2) where C_0 = 1 - ε_r with ε_r being the macroscopic relative permittivity. We neglect K_textxc, such that J^-1 α frack_TF^2 - C_0 G^2ε_r k_TF^2 - C_0 G^2\n\nBy 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 ρ_textspin in the same way.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.DielectricModel","page":"API reference","title":"DFTK.DielectricModel","text":"A localised dielectric model for χ_0:\n\nsqrtL(x) textIFFT fracC_0 G^24π (1 - C_0 G^2 k_TF^2) textFFT sqrtL(x)\n\nwhere C_0 = 1 - ε_r, L(r) is a real-space localisation function and otherwise the same conventions are used as in DielectricMixing.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.ElementCohenBergstresser-Tuple{Any}","page":"API reference","title":"DFTK.ElementCohenBergstresser","text":"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).\n\nkey may be an element symbol (like :Si), an atomic number (e.g. 14) or an element name (e.g. \"silicon\")\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.ElementCoulomb-Tuple{Any}","page":"API reference","title":"DFTK.ElementCoulomb","text":"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\")\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.ElementPsp-Tuple{Any}","page":"API reference","title":"DFTK.ElementPsp","text":"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\")\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.Energies","page":"API reference","title":"DFTK.Energies","text":"A simple struct to contain a vector of energies, and utilities to print them in a nice format.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.Entropy","page":"API reference","title":"DFTK.Entropy","text":"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.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.EtsfFolder-Tuple{AbstractString}","page":"API reference","title":"DFTK.EtsfFolder","text":"Initialize a EtsfFolder from the path to the folder which contains the data in the ETSF Nanoquanta format.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.Ewald","page":"API reference","title":"DFTK.Ewald","text":"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.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.ExternalFromFourier","page":"API reference","title":"DFTK.ExternalFromFourier","text":"External potential from the (unnormalized) Fourier coefficients V(G) G is passed in cartesian coordinates\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.ExternalFromReal","page":"API reference","title":"DFTK.ExternalFromReal","text":"External potential from an analytic function V (in cartesian coordinates). No low-pass filtering is performed.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.FourierMultiplication","page":"API reference","title":"DFTK.FourierMultiplication","text":"Fourier space multiplication, like a kinetic energy term: (Hψ)(G) = multiplier(G) ψ(G)\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.Hartree","page":"API reference","title":"DFTK.Hartree","text":"Hartree term: for a decaying potential V the energy would be\n\n1/2 ∫ρ(x)ρ(y)V(x-y) dxdy\n\nwith 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\n\n1/2 ∫ρ(x)ρ(y) G(x-y) dx dy\n\nwhere 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).\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.KerkerDosMixing","page":"API reference","title":"DFTK.KerkerDosMixing","text":"The same as KerkerMixing, but the Thomas-Fermi wavevector is computed from the current density of states at the Fermi level.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.KerkerMixing","page":"API reference","title":"DFTK.KerkerMixing","text":"Kerker mixing: J^-1 fracα G^2k_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.\n\nNotes:\n\nAbinit calls 1k_TF the dielectric screening length (parameter dielng)\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.Kinetic","page":"API reference","title":"DFTK.Kinetic","text":"Kinetic energy: 1/2 sumn fn ∫ |∇ψn|^2.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.Kpoint","page":"API reference","title":"DFTK.Kpoint","text":"Discretization information for kpoint-dependent quantities such as orbitals. More generally, a kpoint is a block of the Hamiltonian; eg collinear spin is treated by doubling the number of kpoints.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.LdosModel","page":"API reference","title":"DFTK.LdosModel","text":"Represents the LDOS-based χ_0 model\n\nχ_0(r r) = (-D_textloc(r) δ(r r) + D_textloc(r) D_textloc(r) D)\n\nwhere D_textloc is the local density of states and D the density of states. For details see Herbst, Levitt 2020 arXiv:2009.01665\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.LibxcDensity","page":"API reference","title":"DFTK.LibxcDensity","text":"Compute density in real space and its derivatives starting from ρ\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.Magnetic","page":"API reference","title":"DFTK.Magnetic","text":"Magnetic term A(-i). It is assumed (but not checked) that A = 0.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.MagneticFieldOperator","page":"API reference","title":"DFTK.MagneticFieldOperator","text":"Magnetic field operator A⋅(-i∇).\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.Model-Union{Tuple{AbstractArray{T,2}}, Tuple{T}} where T<:Real","page":"API reference","title":"DFTK.Model","text":"Model(lattice; n_electrons, atoms, magnetic_moments, terms, temperature,\n smearing, spin_polarization, symmetry)\n\nCreates the physical specification of a model (without any discretization information).\n\nn_electrons is taken from atoms if not specified\n\nspin_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.\n\nmagnetic_moments is only used to determine the symmetry and the spin_polarization; it is not stored inside the datastructure.\n\nsmearing is Fermi-Dirac if temperature is non-zero, none otherwise\n\nThe 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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.NonlocalOperator","page":"API reference","title":"DFTK.NonlocalOperator","text":"Nonlocal operator in Fourier space in Kleinman-Bylander format, defined by its projectors P matrix and coupling terms D: Hψ = PDP' ψ\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.NoopOperator","page":"API reference","title":"DFTK.NoopOperator","text":"Noop operation: don't do anything. Useful for energy terms that don't depend on the orbitals at all (eg nuclei-nuclei interaction).\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.NoopTerm","page":"API reference","title":"DFTK.NoopTerm","text":"A term with a constant zero energy.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.PlaneWaveBasis","page":"API reference","title":"DFTK.PlaneWaveBasis","text":"A plane-wave discretized Model. Normalization conventions:\n\nThings 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 xsqrt(unit_cell_volume). This is so that, eg norm(ψ) = 1 gives the correct normalization. This also holds for the density and the potentials.\nQuantities expressed on the real-space grid are in actual values.\n\nG_to_r and r_to_G convert between these representations.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.PlaneWaveBasis-2","page":"API reference","title":"DFTK.PlaneWaveBasis","text":"Creates a new basis identical to basis, but with a different set of kpoints\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.PlaneWaveBasis-Tuple{Model,Number}","page":"API reference","title":"DFTK.PlaneWaveBasis","text":"Creates a PlaneWaveBasis using the kinetic energy cutoff Ecut and a Monkhorst-Pack kpoint 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_size_from_minimal_spacing with a minimal spacing of 2π * 0.022 per Bohr.\n\nIf use_symmetry is true (default) the symmetries of the crystal are used to reduce the number of k-Points which are treated explicitly. In this case all guess densities and potential functions must agree with the crystal symmetries or the result is undefined.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.PlaneWaveBasis-Tuple{PlaneWaveBasis}","page":"API reference","title":"DFTK.PlaneWaveBasis","text":"\" Convert a basis into one that uses or doesn't use BZ symmetrization Mainly useful for debug purposes (e.g. in cases we don't want to bother with symmetry)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.PowerNonlinearity","page":"API reference","title":"DFTK.PowerNonlinearity","text":"Power nonlinearity, with energy C ∫ρ^α where ρ is the total density\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.PreconditionerNone","page":"API reference","title":"DFTK.PreconditionerNone","text":"No preconditioning\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.PreconditionerTPA","page":"API reference","title":"DFTK.PreconditionerTPA","text":"(simplified version of) Tetter-Payne-Allan preconditioning ↑ M.P. Teter, M.C. Payne and D.C. Allan, Phys. Rev. B 40, 12255 (1989).\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.PspCorrection","page":"API reference","title":"DFTK.PspCorrection","text":"Pseudopotential correction energy. TODO discuss the need for this.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.PspHgh-Union{Tuple{T}, Tuple{Any,Any,Any,Any,Array{Array{T,2},1}}} where T","page":"API reference","title":"DFTK.PspHgh","text":"PspHgh(Zion::Number, rloc::Number, cloc::Vector, rp::Vector, h::Vector;\n identifier=\"\", description=\"\")\n\nConstruct 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).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.RealFourierArray","page":"API reference","title":"DFTK.RealFourierArray","text":"A structure to facilitate manipulations of an array of real-space type T, in both real and fourier space. Create with from_real or from_fourier, and access with A.real and A.fourier.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.RealFourierOperator","page":"API reference","title":"DFTK.RealFourierOperator","text":"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.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.RealSpaceMultiplication","page":"API reference","title":"DFTK.RealSpaceMultiplication","text":"Real space multiplication by a potential: (Hψ)(r) V(r) ψ(r)\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.SimpleMixing","page":"API reference","title":"DFTK.SimpleMixing","text":"Simple mixing: J^-1 α\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.Xc","page":"API reference","title":"DFTK.Xc","text":"Exchange-correlation term, defined by a list of functionals and usually evaluated through libxc.\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.χ0Mixing","page":"API reference","title":"DFTK.χ0Mixing","text":"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).\n\n\n\n\n\n","category":"type"},{"location":"api/#DFTK.CROP","page":"API reference","title":"DFTK.CROP","text":"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.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.DOS-Tuple{Any,Any,Any}","page":"API reference","title":"DFTK.DOS","text":"Total density of states at energy ε\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.G_to_r!-Tuple{AbstractArray{T,3} where T,PlaneWaveBasis,AbstractArray{T,3} where T}","page":"API reference","title":"DFTK.G_to_r!","text":"In-place version of G_to_r.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.G_to_r-Tuple{PlaneWaveBasis,AbstractArray{T,3} where T}","page":"API reference","title":"DFTK.G_to_r","text":"G_to_r(basis::PlaneWaveBasis, [kpt::Kpoint, ] f_fourier)\n\nPerform 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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.G_vectors-Tuple{Any}","page":"API reference","title":"DFTK.G_vectors","text":"Return the list of wave vectors (integer coordinates) for the cubic basis set.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.G_vectors-Tuple{Kpoint}","page":"API reference","title":"DFTK.G_vectors","text":"The list of G vectors of a given basis or kpoint, in reduced coordinates.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.G_vectors_cart-Tuple{Kpoint}","page":"API reference","title":"DFTK.G_vectors_cart","text":"The list of G vectors of a given basis or kpoint, in cartesian coordinates.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.HybridMixing-Tuple{}","page":"API reference","title":"DFTK.HybridMixing","text":"The model for the susceptibility is\n\nbeginaligned\n χ_0(r r) = (-D_textloc(r) δ(r r) + D_textloc(r) D_textloc(r) D) \n + sqrtL(x) textIFFT fracC_0 G^24π (1 - C_0 G^2 k_TF^2) textFFT sqrtL(x)\nendaligned\n\nwhere C_0 = 1 - ε_r, D_textloc 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 localisation function L(r). For details see Herbst, Levitt 2020 arXiv:2009.01665\n\nImportant kwargs passed on to χ0Mixing\n\nα: Damping parameter\nRPA: Is the random-phase approximation used for the kernel (i.e. only Hartree kernel is used and not XC kernel)\nverbose: Run the GMRES in verbose mode.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.LDOS-NTuple{4,Any}","page":"API reference","title":"DFTK.LDOS","text":"Local density of states, in real space\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.NOS-Tuple{Any,Any,Any}","page":"API reference","title":"DFTK.NOS","text":"NOS(ε, basis, eigenvalues; smearing=basis.model.smearing,\n temperature=basis.model.temperature)\n\nThe number of Kohn-Sham states in a temperature window of width temperature around the energy ε contributing to the DOS at temperature T.\n\nThis 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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.ScfConvergenceDensity-Tuple{Any}","page":"API reference","title":"DFTK.ScfConvergenceDensity","text":"Flag convergence by using the L2Norm of the change between input density and unpreconditioned output density (ρout)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.ScfConvergenceEnergy-Tuple{Any}","page":"API reference","title":"DFTK.ScfConvergenceEnergy","text":"Flag convergence as soon as total energy change drops below tolerance\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.ScfDefaultCallback-Tuple{}","page":"API reference","title":"DFTK.ScfDefaultCallback","text":"Default callback function for self_consistent_field, which prints a convergence table\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.ScfDiagtol-Tuple{}","page":"API reference","title":"DFTK.ScfDiagtol","text":"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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.ScfPlotTrace","page":"API reference","title":"DFTK.ScfPlotTrace","text":"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!.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.ScfSaveCheckpoints","page":"API reference","title":"DFTK.ScfSaveCheckpoints","text":"Adds simplistic checkpointing to a DFTK self-consistent field calculation.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.add_response_from_band!-NTuple{11,Any}","page":"API reference","title":"DFTK.add_response_from_band!","text":"Adds the term (f'ₙ δεₙ |ψₙ|² + 2Re fₙ ψₙ * δψₙ to δρ_{k} where δψₙ is computed from δV partly using the known, computed states and partly by solving the Sternheimer equation (if sternheimer_contribution=true).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.anderson","page":"API reference","title":"DFTK.anderson","text":"Anderson-accelerated root-finding iteration for finding a root of 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.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.apply_kernel","page":"API reference","title":"DFTK.apply_kernel","text":"apply_kernel(basis::PlaneWaveBasis, dρ, dρspin=nothing; kwargs...)\n\nComputes the potential response to a perturbation (dρ, dρspin) in real space. Returns the array [dV_α, dV_β] for collinear spin-polarized calculations, else the array [dV_{tot}].\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.apply_ksymop-Tuple{Any,Any,Any,Union{AbstractArray{T,1}, AbstractArray{T,2}} where T}","page":"API reference","title":"DFTK.apply_ksymop","text":"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 kpoint).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.apply_ksymop-Tuple{Any,RealFourierArray}","page":"API reference","title":"DFTK.apply_ksymop","text":"Apply a k-point symmetry operation (the tuple (S, τ)) to a partial density.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.apply_χ0","page":"API reference","title":"DFTK.apply_χ0","text":"Returns the change in density δρ for a given δV. Drop all non-diagonal terms with (f(εn)-f(εm))/(εn-εm) factor less than droptol. If sternheimer_contribution is false, only compute excitations inside the provided orbitals.\n\nNote: 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 before using this function.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.atom_decay_length-Tuple{Any,Any}","page":"API reference","title":"DFTK.atom_decay_length","text":"Get the lengthscale of the valence density for an atom with n_elec_core core and n_elec_valence valence electrons. ```\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.build_fft_plans-Tuple{Union{Type{Float32}, Type{Float64}},Any}","page":"API reference","title":"DFTK.build_fft_plans","text":"Plan a FFT of type T and size fft_size, spending some time on finding an optimal algorithm. Both an inplace and an out-of-place FFT plan are returned.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.build_form_factors-Tuple{Any,Any}","page":"API reference","title":"DFTK.build_form_factors","text":"Build form factors (Fourier transforms of projectors) for an atom centered at 0.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.build_projection_vectors_-Union{Tuple{T}, Tuple{PlaneWaveBasis{T},Any,Kpoint}} where T","page":"API reference","title":"DFTK.build_projection_vectors_","text":"Build projection vectors for a atoms array generated by term_nonlocal\n\nHat = sumij Cij |pi> <pj| Hper = sumR sumij Cij |pi(x-R)> <pj(x-R)| = sumR sum_ij Cij |pi(x-R)> <pj(x-R)|\n\n<ekG'|Hper|ekG> = ... = 1/Ω sumij Cij pihat(k+G') pjhat(k+G)^*\n\nwhere pihat(q) = ∫_R^3 pi(r) e^{-iqr} dr\n\nWe store 1/√Ω pihat(k+G) in proj_vectors.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.bzmesh_ir_wedge-Tuple{Any,Any}","page":"API reference","title":"DFTK.bzmesh_ir_wedge","text":" bzmesh_ir_wedge(kgrid_size, symmetries; kshift=[0, 0, 0])\n\nConstruct 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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.bzmesh_uniform-Tuple{Any}","page":"API reference","title":"DFTK.bzmesh_uniform","text":"bzmesh_uniform(kgrid_size; kshift=[0, 0, 0])\n\nConstruct 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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.charge_ionic-Tuple{DFTK.Element}","page":"API reference","title":"DFTK.charge_ionic","text":"Return the total ionic charge of an atom type (nuclear charge - core electrons)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.charge_nuclear-Tuple{DFTK.Element}","page":"API reference","title":"DFTK.charge_nuclear","text":"Return the total nuclear charge of an atom type\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.clear_without_conjugate!-Union{Tuple{AbstractArray{T,3}}, Tuple{T}} where T<:Complex","page":"API reference","title":"DFTK.clear_without_conjugate!","text":"Zero all elements of a Fourier space array which have no complex-conjugate partner and may thus lead to an imaginary component in real space (after an iFFT).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.compute_current-Tuple{PlaneWaveBasis,Any,Any}","page":"API reference","title":"DFTK.compute_current","text":"Computes the probability (not charge) current, ∑ fn Im(ψn* ∇ψn)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.compute_density-Tuple{PlaneWaveBasis,Any,Any}","page":"API reference","title":"DFTK.compute_density","text":"compute_density(basis::PlaneWaveBasis, ψ::AbstractVector, occupation::AbstractVector)\n\nCompute the density and spin 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. If the Model underlying the basis is not collinear the spin density is nothing.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.compute_kernel-Union{Tuple{PlaneWaveBasis{T}}, Tuple{T}} where T","page":"API reference","title":"DFTK.compute_kernel","text":"compute_kernel(basis::PlaneWaveBasis; kwargs...)\n\nComputes 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, which are:\n\nleft(beginarraycc\n K_α texttot K_α textspin\n K_β texttot K_β textspin\nendarrayright)\n\ni.e. corresponding to a mapping (ρ_texttot ρ_textspin)^T = (ρ_α + ρ_β ρ_α - ρ_β)^T (V_α V_β)^T.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.compute_partial_density-NTuple{4,Any}","page":"API reference","title":"DFTK.compute_partial_density","text":"Compute the partial density at the indicated k-Point and return it (in Fourier space).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.compute_χ0-Tuple{Any}","page":"API reference","title":"DFTK.compute_χ0","text":"Compute the independent-particle susceptibility. Will blow up for large systems. Drop all non-diagonal terms with (f(εn)-f(εm))/(εn-εm) factor less than droptol. 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:\n\nleft(beginarraycc\n (χ_0)_texttot α (χ_0)_texttot β \n (χ_0)_textspin α (χ_0)_textspin β\nendarrayright)\n\ni.e. corresponding to a mapping ``(Vα, Vβ)^T ↦ (ρ\\text{tot}, ρ\\text{spin})^T = (ρα + ρβ, ρα - ρβ)^T.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.datadir_psp-Tuple{}","page":"API reference","title":"DFTK.datadir_psp","text":"Return the data directory with pseudopotential files\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.default_spin_polarization-Tuple{Any}","page":"API reference","title":"DFTK.default_spin_polarization","text":":none if no element has a magnetic moment, else :collinear or :full\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.default_symmetries-NTuple{5,Any}","page":"API reference","title":"DFTK.default_symmetries","text":"Default logic to determine the symmetry operations to be used in the model.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.determine_fft_size-Union{Tuple{T}, Tuple{AbstractArray{T,2},Any}} where T","page":"API reference","title":"DFTK.determine_fft_size","text":"Determine the minimal grid size for the cubic basis set to be able to represent product of orbitals (with the default supersampling=2).\n\nOptionally optimize the grid afterwards for the FFT procedure by ensuring factorization into small primes.\n\nThe function will determine the smallest parallelepiped containing the wave vectors G^22 leq E_textcut textsupersampling^2. For an exact representation of the density resulting from wave functions represented in the spherical basis sets, supersampling should be at least 2.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.diagonalize_all_kblocks-Tuple{Any,Hamiltonian,Int64}","page":"API reference","title":"DFTK.diagonalize_all_kblocks","text":"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)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.direct_minimization-Tuple{PlaneWaveBasis}","page":"API reference","title":"DFTK.direct_minimization","text":"Computes the ground state by direct minimization. kwargs... are passed to Optim.Options(). Note that the resulting ψ are not necessarily eigenvectors of the Hamiltonian.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.divergence_real-Tuple{Any,Any}","page":"API reference","title":"DFTK.divergence_real","text":"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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.energy_ewald-Tuple{Any,Any,Any}","page":"API reference","title":"DFTK.energy_ewald","text":"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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.energy_psp_correction-Tuple{Model}","page":"API reference","title":"DFTK.energy_psp_correction","text":"energy_psp_correction(model)\n\nCompute the correction term for properly modelling the interaction of the pseudopotential core with the compensating background charge induced by the Ewald term.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.eval_psp_energy_correction-Tuple{Any,PspHgh,Any}","page":"API reference","title":"DFTK.eval_psp_energy_correction","text":"eval_psp_energy_correction([T=Float64,] psp, n_electrons)\n\nEvaluate 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.\n\nNotice: 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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.eval_psp_local_fourier-Union{Tuple{T}, Tuple{PspHgh,T}} where T<:Real","page":"API reference","title":"DFTK.eval_psp_local_fourier","text":"eval_psp_local_fourier(psp, q)\n\nEvaluate the local part of the pseudopotential in reciprocal space.\n\nThis function computes V(q) = ∫R^3 Vloc(r) e^{-iqr} dr = 4π ∫{R+} sin(qr)/q r e^{-iqr} dr\n\nGTH98 except they do it with plane waves normalized by 1/sqrt(Ω).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.eval_psp_local_real-Union{Tuple{T}, Tuple{PspHgh,T}} where T<:Real","page":"API reference","title":"DFTK.eval_psp_local_real","text":"eval_psp_local_real(psp, r)\n\nEvaluate the local part of the pseudopotential in real space.\n\nGTH98\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.eval_psp_projection_radial-Union{Tuple{T}, Tuple{PspHgh,Any,Any,T}} where T<:Real","page":"API reference","title":"DFTK.eval_psp_projection_radial","text":"eval_psp_projection_radial(psp::PspHgh, i, l, q::Number)\n\nEvaluate the radial part of the i-th projector for angular momentum l at the reciprocal vector with modulus q.\n\np(q) = ∫{R+} r^2 p(r) jl(q r) dr\n\nHGH98 except they do it with plane waves normalized by 1/sqrt(Ω).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.eval_psp_projection_radial_real-Union{Tuple{T}, Tuple{PspHgh,Any,Any,T}} where T<:Real","page":"API reference","title":"DFTK.eval_psp_projection_radial_real","text":"eval_psp_projection_radial_real(psp::PspHgh, i, l, q::Real)\n\nEvaluate the radial part of the i-th projector for angular momentum l in real-space at the vector with modulus r. HGH98\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.fermi_level-Tuple{Any,Any}","page":"API reference","title":"DFTK.fermi_level","text":"Find the Fermi level.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.filled_occupation-Tuple{Any}","page":"API reference","title":"DFTK.filled_occupation","text":"Maximal occupation of a state (2 for non-spin-polarized electrons, 1 otherwise).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.find_irreducible_kpoints-Tuple{Any,Any,Any}","page":"API reference","title":"DFTK.find_irreducible_kpoints","text":"Implements a primitive search to find an irreducible subset of kpoints amongst the provided kpoints.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.find_occupation-Union{Tuple{T}, Tuple{PlaneWaveBasis{T},Any}} where T","page":"API reference","title":"DFTK.find_occupation","text":"Find the occupation and Fermi level.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.find_occupation_bandgap-Tuple{Any,Any}","page":"API reference","title":"DFTK.find_occupation_bandgap","text":"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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.gaussian_superposition-Union{Tuple{T}, Tuple{PlaneWaveBasis{T},Any}} where T","page":"API reference","title":"DFTK.gaussian_superposition","text":"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\n\nhatρ(G) = textcoefficient expleft(-(2π textlength G)^2right)\n\nand are placed at position (in fractional coordinates).\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.guess_density-Tuple{PlaneWaveBasis}","page":"API reference","title":"DFTK.guess_density","text":"guess_density(basis)\n\nBuild a superposition of atomic densities (SAD) guess density.\n\nWe 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\n\nhatρ(G) = Z expleft(-(2π textlength G)^2right)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.guess_spin_density","page":"API reference","title":"DFTK.guess_spin_density","text":"guess_spin_density(basis, magnetic_moments)\n\nForm a spin density guess. The magnetic moments should be specified in units of μ_B.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.hamiltonian_with_total_potential-Tuple{Hamiltonian,Any}","page":"API reference","title":"DFTK.hamiltonian_with_total_potential","text":"Returns a new Hamiltonian with local potential replaced by the given one\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.index_G_vectors-Union{Tuple{T}, Tuple{PlaneWaveBasis,AbstractArray{T,1}}} where T<:Integer","page":"API reference","title":"DFTK.index_G_vectors","text":"Return the index tuple I such that G_vectors(basis)[I] == G or the index i such that G_vectors(kpoint)[i] == G. Returns nothing if outside the range of valid wave vectors.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.interpolate_blochwave-Tuple{Any,Any,Any}","page":"API reference","title":"DFTK.interpolate_blochwave","text":"Interpolate Bloch wave between two basis sets. Limited feature set. Currently only interpolation to a bigger grid (larger Ecut) on the same lattice supported.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.interpolate_density-Tuple{RealFourierArray,PlaneWaveBasis}","page":"API reference","title":"DFTK.interpolate_density","text":"Interpolate a function expressed in a basis b_in to a basis b_out This interpolation uses a very basic real-space algorithm, and makes a DWIM-y attempt to take into account the fact that bout can be a supercell of bin\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.interpolate_kpoint-Tuple{Union{AbstractArray{T,1}, AbstractArray{T,2}} where T,Kpoint,Kpoint}","page":"API reference","title":"DFTK.interpolate_kpoint","text":"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\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.is_metal","page":"API reference","title":"DFTK.is_metal","text":"is_metal(band_data, εF, tol)\n\nDetermine whether the provided bands indicate the material is a metal, i.e. where bands are cut by the Fermi level.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.kgrid_monkhorst_pack-Tuple{Any}","page":"API reference","title":"DFTK.kgrid_monkhorst_pack","text":"Construct the coordinates of the kpoints in a (shifted) Monkorst-Pack grid\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.kgrid_size_from_minimal_spacing","page":"API reference","title":"DFTK.kgrid_size_from_minimal_spacing","text":"Selects a kgrid_size to ensure a minimal spacing (in inverse Bohrs) between kpoints. Default is 2π * 004 AA^-1.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.krange_spin-Tuple{PlaneWaveBasis,Integer}","page":"API reference","title":"DFTK.krange_spin","text":"Return the index range of k-points that have a particular spin component.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.lda_c_vwn!-Tuple{Any}","page":"API reference","title":"DFTK.lda_c_vwn!","text":"LDA correlation according to Vosko Wilk,and Nusair, (DOI 10.1139/p80-159)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.lda_x!-Tuple{Any}","page":"API reference","title":"DFTK.lda_x!","text":"LDA Slater exchange (DOI: 10.1017/S0305004100016108 and 10.1007/BF01340281)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.libxc_symmetric_index-Tuple{Any,Any}","page":"API reference","title":"DFTK.libxc_symmetric_index","text":"Translation from a double index (i, j) of two symmetric axes to the linear index over the non-redundant components used in libxc\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.list_psp","page":"API reference","title":"DFTK.list_psp","text":"list_psp(element; functional, family, core, datadir_psp)\n\nList the pseudopotential files known to DFTK. Allows various ways to restrict the displayed files.\n\nExamples\n\njulia> list_psp(family=\"hgh\")\n\nwill list all HGH-type pseudopotentials and\n\njulia> list_psp(family=\"hgh\", functional=\"lda\")\n\nwill only list those for LDA (also known as Pade in this context) and\n\njulia> list_psp(:O, core=:semicore)\n\nwill list all oxygen semicore pseudopotentials known to DFTK.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.load_atoms-Tuple{Any,EtsfFolder}","page":"API reference","title":"DFTK.load_atoms","text":"Load a DFTK-compatible atoms object from the ETSF folder. Use the scalar type T to represent the data.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.load_atoms_pymatgen-Tuple{Any,PyCall.PyObject}","page":"API reference","title":"DFTK.load_atoms_pymatgen","text":"Load a DFTK-compatible atoms representation from a supported pymatgen object. All atoms are using a Coulomb model.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.load_basis-Tuple{Any,EtsfFolder}","page":"API reference","title":"DFTK.load_basis","text":"Load a DFTK-compatible basis object from the ETSF folder. Use the scalar type T to represent the data.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.load_density-Tuple{Any,EtsfFolder}","page":"API reference","title":"DFTK.load_density","text":"Load a DFTK-compatible density object from the ETSF folder. Use the scalar type T to represent the data.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.load_lattice-Tuple{Any,EtsfFolder}","page":"API reference","title":"DFTK.load_lattice","text":"Load a DFTK-compatible lattice object from the ETSF folder\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.load_lattice-Tuple{Any,PyCall.PyObject}","page":"API reference","title":"DFTK.load_lattice","text":"Load a DFTK-compatible lattice object from a supported python object (e.g. pymatgen or ASE)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.load_model-Tuple{Any,EtsfFolder}","page":"API reference","title":"DFTK.load_model","text":"Load a DFTK-compatible model object from the ETSF folder. Use the scalar type T to represent the data.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.load_psp-Tuple{AbstractString}","page":"API reference","title":"DFTK.load_psp","text":"load_psp(key; datadir_psp)\n\nLoad 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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.load_scfres","page":"API reference","title":"DFTK.load_scfres","text":"load_scfres(filename)\n\nLoad back an scfres, which has previously been stored with save_scfres. Note the warning in save_scfres.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.local_potential_fourier-Tuple{DFTK.Element,AbstractArray{T,1} where T}","page":"API reference","title":"DFTK.local_potential_fourier","text":"Radial local potential, in Fourier space: V(q) = int_{R^3} V(x) e^{-iqx} dx.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.local_potential_real-Tuple{DFTK.Element,AbstractArray{T,1} where T}","page":"API reference","title":"DFTK.local_potential_real","text":"Radial local potential, in real space.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.model_DFT-Tuple{AbstractArray{T,2} where T,Array{T,1} where T,Any}","page":"API reference","title":"DFTK.model_DFT","text":"Build a DFT model from the specified atoms, with the specified functionals.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.model_LDA-Tuple{AbstractArray{T,2} where T,Array{T,1} where T}","page":"API reference","title":"DFTK.model_LDA","text":"Build an LDA model (Teter93 parametrization) from the specified atoms.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.model_PBE-Tuple{AbstractArray{T,2} where T,Array{T,1} where T}","page":"API reference","title":"DFTK.model_PBE","text":"Build an PBE-GGA model from the specified atoms.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.model_atomic-Tuple{AbstractArray{T,2} where T,Array{T,1} where T}","page":"API reference","title":"DFTK.model_atomic","text":"Convenience constructor, which builds a standard atomic (kinetic + atomic potential) model. Use extra_terms to add additional terms.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.n_elec_core-Tuple{DFTK.Element}","page":"API reference","title":"DFTK.n_elec_core","text":"Return the number of core electrons\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.n_elec_valence-Tuple{DFTK.Element}","page":"API reference","title":"DFTK.n_elec_valence","text":"Return the number of valence electrons\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.next_density-Tuple{Hamiltonian}","page":"API reference","title":"DFTK.next_density","text":"Obtain new density ρ by diagonalizing ham.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.normalize_kpoint_coordinate-Tuple{Real}","page":"API reference","title":"DFTK.normalize_kpoint_coordinate","text":"Bring kpoint coordinates into the range [-0.5, 0.5)\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.parse_hgh_file-Tuple{Any}","page":"API reference","title":"DFTK.parse_hgh_file","text":"parse_hgh_file(path; identifier=\"\")\n\nParse an HGH pseudopotential file and construct the PspHgh object. If identifier is given, this identifier will be set.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.plot_bandstructure-NTuple{4,Any}","page":"API reference","title":"DFTK.plot_bandstructure","text":"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. Unlike the rest of DFTK bands energies are plotted in :eV unless a different unit is selected.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.plot_dos","page":"API reference","title":"DFTK.plot_dos","text":"Plot the density of states over a reasonable range\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.psp_local_polynomial","page":"API reference","title":"DFTK.psp_local_polynomial","text":"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_textloc q and Q is a polynomial of at most degree 8. This function returns Q.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.psp_projection_radial_polynomial","page":"API reference","title":"DFTK.psp_projection_radial_polynomial","text":"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.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.qcut_psp_local-Tuple{Any,PspHgh}","page":"API reference","title":"DFTK.qcut_psp_local","text":"Estimate an upper bound for the argument q after which abs(eval_psp_local_fourier(psp, q)) is a strictly decreasing function.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.qcut_psp_projection_radial-Tuple{Any,PspHgh,Any,Any}","page":"API reference","title":"DFTK.qcut_psp_projection_radial","text":"Estimate an upper bound for the argument q after which eval_psp_projection_radial(psp, q) is a strictly decreasing function.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.r_to_G!-Tuple{AbstractArray{T,3} where T,PlaneWaveBasis,AbstractArray{T,3} where T}","page":"API reference","title":"DFTK.r_to_G!","text":"In-place version of r_to_G!. NOTE: If kpt is given, not only f_fourier but also f_real is overwritten.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.r_to_G-Tuple{PlaneWaveBasis,AbstractArray{T,3} where T}","page":"API reference","title":"DFTK.r_to_G","text":"r_to_G(basis::PlaneWaveBasis, [kpt::Kpoint, ] f_real)\n\nPerform 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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.r_vectors-Union{Tuple{PlaneWaveBasis{T}}, Tuple{T}} where T","page":"API reference","title":"DFTK.r_vectors","text":"Return the list of r vectors, in reduced coordinates. By convention, this is in [0,1)^3.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.r_vectors_cart-Tuple{PlaneWaveBasis}","page":"API reference","title":"DFTK.r_vectors_cart","text":"Return the list of r vectors, in cartesian coordinates.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.run_abinit_scf-Tuple{Model,Any}","page":"API reference","title":"DFTK.run_abinit_scf","text":"Run an SCF in ABINIT starting from a DFTK Model and some extra parameters. Write the result to the output directory in ETSF Nanoquanta format and return the EtsfFolder object.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.run_abinit_scf-Tuple{PyCall.PyObject,Any}","page":"API reference","title":"DFTK.run_abinit_scf","text":"Run an SCF in ABINIT starting from the input file infile represented as a abipy.abilab.AbinitInput python object. Write the result to the output directory in ETSF Nanoquanta format and return the result as an EtsfFolder object.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.save_scfres","page":"API reference","title":"DFTK.save_scfres","text":"save_scfres(filename, scfres)\n\nSave an scfres obtained from self_consistent_field to a JLD2 file.\n\nwarning: No compatibility guarantees\nNo guarantees are made with respect to this function at this point. It may change incompatibly between DFTK versions or stop working / be removed in the future.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.scf_damping_solver","page":"API reference","title":"DFTK.scf_damping_solver","text":"Create a damped SCF solver updating the density as x = β * x_new + (1 - β) * x\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.scf_nlsolve_solver","page":"API reference","title":"DFTK.scf_nlsolve_solver","text":"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.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.select_eigenpairs_all_kblocks-Tuple{Any,Any}","page":"API reference","title":"DFTK.select_eigenpairs_all_kblocks","text":"Function to select a subset of eigenpairs on each k-Point. Works on the Tuple returned by diagonalize_all_kblocks.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.self_consistent_field-Tuple{PlaneWaveBasis}","page":"API reference","title":"DFTK.self_consistent_field","text":"Solve the Kohn-Sham equations with a SCF algorithm, starting at ρ.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.spglib_atoms","page":"API reference","title":"DFTK.spglib_atoms","text":"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.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.spin_components-Tuple{Symbol}","page":"API reference","title":"DFTK.spin_components","text":"Explicit spin components of the KS orbitals and the density\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.standardize_atoms","page":"API reference","title":"DFTK.standardize_atoms","text":"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.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.symmetries_preserving_kgrid-Tuple{Any,Any}","page":"API reference","title":"DFTK.symmetries_preserving_kgrid","text":"Filter out the symmetry operations that respect the symmetries of the discrete BZ grid\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.symmetrize-Tuple{RealFourierArray}","page":"API reference","title":"DFTK.symmetrize","text":"Symmetrize a RealFourierArray by applying all the model symmetries (by default) and forming the average.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.symmetry_operations","page":"API reference","title":"DFTK.symmetry_operations","text":"Return the k-point symmetry operations associated to a lattice and atoms.\n\n\n\n\n\n","category":"function"},{"location":"api/#DFTK.total_local_potential-Tuple{Hamiltonian}","page":"API reference","title":"DFTK.total_local_potential","text":"Get the total local potential of the given Hamiltonian, in real space.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.unit_to_au-Tuple{Symbol}","page":"API reference","title":"DFTK.unit_to_au","text":"unit_to_ao(symbol)\n\nGet the factor converting from the unit symbol to atomic units. E.g. unit_to_au(:eV) returns the conversion factor from electron volts to Hartree.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.ylm_real-Union{Tuple{T}, Tuple{Integer,Integer,AbstractArray{T,1}}} where T","page":"API reference","title":"DFTK.ylm_real","text":"Returns the (l,m) real spherical harmonic Ylm(r). Consistent with https://en.wikipedia.org/wiki/Tableofsphericalharmonics#Realsphericalharmonics\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.@timing-Tuple","page":"API reference","title":"DFTK.@timing","text":"Shortened version of the @timeit macro from TimerOutputs, which writes to the DFTK timer.\n\n\n\n\n\n","category":"macro"},{"location":"api/#DFTK.@timing_seq-Tuple","page":"API reference","title":"DFTK.@timing_seq","text":"Similar to @timing, but disabled in parallel runs. Should be used to time threaded regions, since TimerOutputs is not thread-safe and breaks otherwise.\n\n\n\n\n\n","category":"macro"},{"location":"api/#DFTK.Smearing.entropy-Tuple{DFTK.Smearing.SmearingFunction,Any}","page":"API reference","title":"DFTK.Smearing.entropy","text":"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.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.Smearing.occupation-Tuple{DFTK.Smearing.SmearingFunction,Any}","page":"API reference","title":"DFTK.Smearing.occupation","text":"Occupation at x, where in practice x = (ε - εF) / T.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.Smearing.occupation_derivative-Tuple{DFTK.Smearing.SmearingFunction,Any}","page":"API reference","title":"DFTK.Smearing.occupation_derivative","text":"Derivative of the occupation function, approximation to minus the delta function.\n\n\n\n\n\n","category":"method"},{"location":"api/#DFTK.Smearing.occupation_divided_difference-Tuple{DFTK.Smearing.SmearingFunction,Any,Any,Any,Any}","page":"API reference","title":"DFTK.Smearing.occupation_divided_difference","text":"(f(x) - f(y))/(x - y), computed stably in the case where x and y are close\n\n\n\n\n\n","category":"method"},{"location":"guide/parallelisation/#Timings-and-parallelisation","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"","category":"section"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"This section summarises the options DFTK offers to monitor and influence performance of the code.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"using DFTK\na = 10.26 # Silicon lattice constant in Bohr\nlattice = a / 2 * [[0 1 1.];\n [1 0 1.];\n [1 1 0.]]\nSi = ElementPsp(:Si, psp=load_psp(\"hgh/lda/Si-q4\"))\natoms = [Si => [ones(3)/8, -ones(3)/8]]\n\nmodel = model_LDA(lattice, atoms)\nkgrid = [2, 2, 2]\nEcut = 5\nbasis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\n\nimport TimerOutputs\nTimerOutputs.reset_timer!(DFTK.timer)\nscfres = self_consistent_field(basis, tol=1e-8)","category":"page"},{"location":"guide/parallelisation/#Built-in-timing-measurements","page":"Timings and parallelisation","title":"Built-in timing measurements","text":"","category":"section"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"By default DFTK uses TimerOutputs.jl to record timings, memory allocations and the number of calls for selected routines inside the code. These numbers are accessible in the object DFTK.timer. Since the timings are automatically accumulated inside this datastructure, any timing measurement should first reset this timer before running the calculation of interest.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"For example to measure the timing of an SCF:","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"import TimerOutputs\n\nTimerOutputs.reset_timer!(DFTK.timer)\nscfres = self_consistent_field(basis, tol=1e-8)\n\nDFTK.timer","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"The output produced when printing or displaying the DFTK.timer now shows a nice table summarising total time and allocations as well as a breakdown over individual routines.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"note: Timing measurements and stack traces\nTiming measurements have the unfortunate disadvantage that they alter the way stack traces look making it sometimes harder to find errors when debugging. For this reason timing measurements can be disabled completely (i.e. not even compiled into the code) by setting the environment variable DFTK_TIMING to \"0\" or \"false\". For this to take effect recompiling all DFTK (including the precompile cache) is needed.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"note: Timing measurements and threading\nUnfortunately measuring timings in TimerOutputs is not yet thread-safe. Therefore taking timings of threaded parts of the code will be disabled unless you set DFTK_TIMING to \"all\". In this case you must not use Julia threading (see section below) or otherwise undefined behaviour results.","category":"page"},{"location":"guide/parallelisation/#Options-to-influence-threading","page":"Timings and parallelisation","title":"Options to influence threading","text":"","category":"section"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"At the moment DFTK employs shared-memory parallelism using multiple levels of threading which distribute the workload over different k-Points, bands or within an FFT or BLAS call between processors. At its current stage our approach to threading is quite simple and pragmatic, such that we do not yet achieve a great scaling. This should be improved in future versions of the code.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"Finding a good sweet spot between the number of threads to use and the extra performance gained by each additional working core is not always easy, since starting, terminating and synchronising threads takes time as well. Most importantly the best settings wrt. threading depend on both hardware and the problem (e.g. number of bands, k-Points, FFT grid size).","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"For the moment DFTK does not offer an automated selection mechanism of thread and parallelisation threading and just uses the Julia defaults. Since these are rarely good, users are advised to use the timing capabilities of DFTK to experiment with threading for their particular use case before running larger calculations.","category":"page"},{"location":"guide/parallelisation/#FFTW-threads","page":"Timings and parallelisation","title":"FFTW threads","text":"","category":"section"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"For typical small to medium-size calculations in DFTK the largest part of time is spent doing discrete Fourier transforms (about 80 to 90%). For this reason parallelising FFTs can have a large effect on the runtime of larger calculation in DFTK. Unfortunately for scaling of FFT threading for smaller problem sizes and large numbers of threads is not great, such that by default threading in FFTW is even disabled.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"The recommended setting for FFT threading with DFTK is therefore to only use moderate number of FFT threads, something like 2 or 4 and for smaller calculations disable FFT threading completely. To enable parallelisation of FFTs (which is by default disabled), use","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"using FFTW\nFFTW.set_num_threads(N)","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"where N is the number of threads you desire.","category":"page"},{"location":"guide/parallelisation/#BLAS-threads","page":"Timings and parallelisation","title":"BLAS threads","text":"","category":"section"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"All BLAS calls in Julia go through a parallelised OpenBlas or MKL (with MKL.jl. Generally threading in BLAS calls is far from optimal and the default settings can be pretty bad. For example for CPUs with hyper threading enabled, the default number of threads seems to equal the number of virtual cores. Still, BLAS calls typically take second place in terms of the share of runtime they make up (between 10% and 20%). Of note many of these do not take place on matrices of the size of the full FFT grid, but rather only in a subspace (e.g. orthogonalisation, Rayleigh-Ritz, ...) such that parallelisation is either anyway disabled by the BLAS library or not very effective.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"The recommendation is therefore to use the same number of threads as for the FFT threads. You can set the number of BLAS threads by","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"using LinearAlgebra\nBLAS.set_num_threads(N)","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"where N is the number of threads you desire. To check the number of BLAS threads currently used, you can use","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"Int(ccall((BLAS.@blasfunc(openblas_get_num_threads), BLAS.libblas), Cint, ()))","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"or (from Julia 1.6) simply BLAS.get_num_threads().","category":"page"},{"location":"guide/parallelisation/#Julia-threads","page":"Timings and parallelisation","title":"Julia threads","text":"","category":"section"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"On top of FFT and BLAS threading DFTK uses Julia threads (Thread.@threads) in a couple of places to parallelise over k-Points (density computation) or bands (Hamiltonian application). The number of threads used for these aspects is controlled by the environment variable JULIA_NUM_THREADS. To influence the number of Julia threads used, set this variable before starting the Julia process.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"Notice, that Julia threading is applied on top of FFTW and BLAS threading in the sense that the regions parallelised by Julia threads again use parallelised FFT and BLAS calls, such that the effects are not orthogonal. Compared to FFT and BLAS threading the parallelisation implied by using Julia threads tends to scale better, but its effectiveness is limited by the number of bands and the number of irreducible k-Points used in the calculation. Therefore this good scaling quickly diminishes for small to medium systems.","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"The recommended setting is to stick to 2 Julia threads and to use 4 or more Julia threads only for large systems and/or many k-Points. To check the number of Julia threads use Threads.nthreads().","category":"page"},{"location":"guide/parallelisation/#Summary-of-recommended-settings","page":"Timings and parallelisation","title":"Summary of recommended settings","text":"","category":"section"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"Calculation size JULIA_NUM_THREADS FFTW.set_num_threads(N) BLAS.set_num_threads(N)\ntiny 1 1 1\nsmall 2 1 1\nmedium 2 2 2\nlarge 2 4 4","category":"page"},{"location":"guide/parallelisation/","page":"Timings and parallelisation","title":"Timings and parallelisation","text":"Note, that this picture is likely to change in future versions of DFTK and Julia as improvements to threading and parallelisation are made in the language or the code.","category":"page"},{"location":"advanced/conventions/#Notation-and-conventions","page":"Notation and conventions","title":"Notation and conventions","text":"","category":"section"},{"location":"advanced/conventions/#Usage-of-unicode-characters","page":"Notation and conventions","title":"Usage of unicode characters","text":"","category":"section"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"DFTK liberally uses unicode characters to represent Greek characters (e.g. ψ, ρ, ε...). Make sure you use the proper Julia plugins to simplify typing them.","category":"page"},{"location":"advanced/conventions/#symbol-conventions","page":"Notation and conventions","title":"Symbol conventions","text":"","category":"section"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"Reciprocal-space vectors: k for vectors in the Brillouin zone, G for vectors of the reciprocal lattice, q for general vectors\nReal-space vectors: R for lattice vectors, r and x are usually used for unit for vectors in the unit cell or general real-space vectors, respectively. This convention is, however, less consistently applied.\nOmega is the unit cell, and Omega (or sometimes just Omega) is its volume.\nA are the real-space lattice vectors (model.lattice) and B the Brillouin zone lattice vectors (model.recip_lattice).\nThe Bloch waves are\npsi_nk(x) = e^ikcdot x u_nk(x)\nwhere n is the band index and k the k-point. In the code we sometimes use psi and u interchangeably.\nvarepsilon are the eigenvalues, varepsilon_F is the Fermi level.\nrho is the density.\nIn the code we use normalized plane waves:\ne_G(r) = frac 1 sqrtOmega e^i G cdot r\nY^l_m are the complex spherical harmonics, and Y_lm the real ones.\nj_l are the Bessel functions. In particular, j_0(x) = fracsin xx.","category":"page"},{"location":"advanced/conventions/#Units","page":"Notation and conventions","title":"Units","text":"","category":"section"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"In DFTK, atomic units are used throughout, most importantly lengths are in Bohr and energies in Hartree. See wikipedia for a list of conversion factors. Useful conversion factors can also be found in DFTK.units and using DFTK.unit_to_au:","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"import DFTK.units: eV\n10eV # 10eV in Hartree","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"import DFTK.units: Å\n1.2 / Å # 1.2 Bohr in Ångström","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"warning: Differing unit conventions\nDifferent electronic-structure codes use different unit conventions. For example for lattice vectors the common length units are Bohr (used by DFTK) and Ångström (used e.g. by ASE, 1Å ≈ 1.80 Bohr). When setting up a calculation for DFTK one needs to ensure to convert to Bohr and atomic units. For some Python libraries (currently ASE, pymatgen and abipy) DFTK directly ships conversion tools in form of the load_lattice and load_atoms functions, which take care of such conversions. Examples which demonstrate this are Creating slabs with ASE and Creating supercells with pymatgen.","category":"page"},{"location":"advanced/conventions/#conventions-lattice","page":"Notation and conventions","title":"Lattices and lattice vectors","text":"","category":"section"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"Both the real-space lattice (i.e. model.lattice) and reciprocal-space lattice (model.recip_lattice) contain the lattice vectors in columns. If 1D or 2D problems are to be treated these arrays are still 3 times 3 matrices, but contain two or one zero-columns, respectively. The real-space lattice vectors are sometimes referred to by A and the reciprocal-space lattice vectors by B = 2pi A^-T.","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"warning: Row-major versus column-major storage order\nJulia stores matrices as column-major, but other languages (notably Python and C) use row-major ordering. Care therefore needs to be taken to properly transpose the unit cell matrices A before using it with DFTK. For the supported third-party packages load_lattice and load_atoms again handle such conversion automatically.","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"We use the convention that the unit cell in real space is 0 1)^3 in reduced coordinates and the unit cell in reciprocal space (the reducible Brillouin zone) is -12 12)^3.","category":"page"},{"location":"advanced/conventions/#Reduced-and-cartesian-coordinates","page":"Notation and conventions","title":"Reduced and cartesian coordinates","text":"","category":"section"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"Unless denoted otherwise the code uses reduced coordinates for reciprocal-space vectors such as k, G, q or real-space vectors like r and R (see Symbol conventions). One switches to Cartesian coordinates by","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"x_textcart = M x_textred","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"where M is either A / model.lattice (for real-space vectors) or B / model.recip_lattice (for reciprocal-space vectors). A useful relationship is","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"b_textcart cdot a_textcart=2pi b_textred cdot a_textred","category":"page"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"if a and b are real-space and reciprocal-space vectors respectively. Other names for reduced coordinates are integer coordinates (usually for G-vectors) or fractional coordinates (usually for k-points).","category":"page"},{"location":"advanced/conventions/#Normalization-conventions","page":"Notation and conventions","title":"Normalization conventions","text":"","category":"section"},{"location":"advanced/conventions/","page":"Notation and conventions","title":"Notation and conventions","text":"The normalization conventions used in the code is that quantities stored in reciprocal space are coefficients in the e_G basis, and quantities stored in real space use real physical values. This means for instance that wavefunctions in the real space grid are normalized as fracOmegaN sum_r psi(r)^2 = 1 where N is the number of grid points and in reciprocal space its coefficients are ell^2-normalized, see the discussion in section PlaneWaveBasis and plane-wave discretisations where this is demonstrated.","category":"page"},{"location":"advanced/symmetries/#Crystal-symmetries","page":"Crystal symmetries","title":"Crystal symmetries","text":"","category":"section"},{"location":"advanced/symmetries/#Theory","page":"Crystal symmetries","title":"Theory","text":"","category":"section"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"In this discussion we will only describe the situation for a monoatomic crystal mathcal C subset mathbb R^3, the extension being easy. A symmetry of the crystal is an orthogonal matrix widetildeS and a real-space vector widetildetau such that","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"widetildeS mathcalC + widetildetau = mathcalC","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"The symmetries where widetildeS = 1 and widetildetau is a lattice vector are always assumed and ignored in the following.","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"We can define a corresponding unitary operator U on L^2(mathbb R^3) with action","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":" (Uu)(x) = uleft( widetildeS x + widetildetau right)","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"We assume that the atomic potentials are radial and that any self-consistent potential also respects this symmetry, so that U commutes with the Hamiltonian.","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"This operator acts on a plane-wave as","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"beginaligned\n(U e^iqcdot x) (x) = e^iq cdot widetildetau e^i (widetildeS^T q) x\n= e^- i(S q) cdot tau e^i (S q) cdot x\nendaligned","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"where we set","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"beginaligned\nS = widetildeS^T\ntau = -widetildeS^-1widetildetau\nendaligned","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"(these equations being also valid in reduced coordinates).","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"It follows that the Fourier transform satisfies","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"widehatUu(q) = e^- iq cdot tau widehat u(S^-1 q)","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"In particular, if e^ikcdot x u_k(x) is an eigenfunction, then by decomposing u_k over plane-waves e^i G cdot x one can see that e^i(S^T k) cdot x (U u_k)(x) is also an eigenfunction: we can choose","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"u_Sk = U u_k","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"This is used to reduce the computations needed. For a uniform sampling of the Brillouin zone (the reducible k-points), one can find a reduced set of k-points (the irreducible k-points) such that the eigenvectors at the reducible k-points can be deduced from those at the irreducible k-points.","category":"page"},{"location":"advanced/symmetries/#Example","page":"Crystal symmetries","title":"Example","text":"","category":"section"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"using DFTK\na = 10.26\nlattice = a / 2 * [[0 1 1.];\n [1 0 1.];\n [1 1 0.]]\nSi = ElementPsp(:Si, psp=load_psp(\"hgh/lda/Si-q4\"))\natoms = [Si => [ones(3)/8, -ones(3)/8]]\nEcut = 5\nkgrid = [4, 4, 4]","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"Let us demonstrate this in practice. We consider silicon, setup appropriately in the lattice and atoms objects as in Tutorial and to reach a fast execution, we take a small Ecut of 5 and a [4, 4, 4] Monkhorst-Pack grid. First we perform the DFT calculation disabling symmetry handling","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"model = model_LDA(lattice, atoms)\nbasis_nosym = PlaneWaveBasis(model, Ecut; kgrid=kgrid, use_symmetry=false)\nscfres_nosym = @time self_consistent_field(basis_nosym, tol=1e-8)\nnothing # hide","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"and then redo it using symmetry (the default):","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"basis_sym = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\nscfres_sym = @time self_consistent_field(basis_sym, tol=1e-8)\nnothing # hide","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"Clearly both yield the same energy but the version employing symmetry is faster, since less k-points are explicitly treated:","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"(length(basis_sym.kpoints), length(basis_nosym.kpoints))","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"Both SCFs would even agree in the convergence history if exact diagonalization was used for the eigensolver in each step of both SCFs. But since DFTK adjusts this diagtol value adaptively during the SCF to increase performance, a slightly different history is obtained. Try adding the keyword argument determine_diagtol=(args...; kwargs...) -> 1e-8 in each SCF call to fix the diagonalization tolerance to be 1e-8 for all SCF steps, which will result in an almost identical convergence history.","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"We can also explicitly verify both methods to yield the same density:","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"using LinearAlgebra # hide\n(norm(scfres_sym.ρ.real - scfres_nosym.ρ.real),\n norm(values(scfres_sym.energies) .- values(scfres_nosym.energies)))","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"To demonstrate the mapping between k-points due to symmetry, we pick an arbitrary k-point in the irreducible Brillouin zone:","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"ikpt_irred = 2\nkpt_irred_coord = basis_sym.kpoints[ikpt_irred].coordinate\nbasis_sym.ksymops[ikpt_irred]","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"This is a list of all symmetries operations (S tau) that can be used to map this irreducible k-point to reducible k-points. Let's pick the third symmetry operation of this k-point and check.","category":"page"},{"location":"advanced/symmetries/","page":"Crystal symmetries","title":"Crystal symmetries","text":"S, τ = basis_sym.ksymops[ikpt_irred][3]\nkpt_red_coord = S * basis_sym.kpoints[ikpt_irred].coordinate\nikpt_red = findfirst(kcoord -> kcoord ≈ kpt_red_coord,\n [k.coordinate for k in basis_nosym.kpoints])\n[scfres_sym.eigenvalues[ikpt_irred] scfres_nosym.eigenvalues[ikpt_red]]","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/ase.jl\"","category":"page"},{"location":"examples/ase/#Creating-slabs-with-ASE","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"","category":"section"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"ASE is short for the atomistic simulation environment, a Python package to simplify the process of setting up, running and analysing results from atomistic simulations across different programs. Extremely powerful in this respect are the routines this code provides for setting up complicated systems (including surface-adsorption scenarios, defects, nanotubes, etc). See also the ASE installation instructions.","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"This example shows how to use ASE to setup a particular gallium arsenide surface and run the resulting calculation in DFTK. If you are less interested in having access to the full playground of options in DFTK, but more interested in performing analysis in ASE itself, have a look at asedftk. This package provides an ASE-compatible calculator class based on DFTK, such that one may write the usual Python scripts against ASE, but the calculations are still run in DFTK.","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"The particular example we consider the (1, 1, 0) GaAs surface separated by vacuum with the setup slightly adapted from [RCW2001].","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"[RCW2001]: D. Raczkowski, A. Canning, and L. W. Wang Thomas-Fermi charge mixing for obtaining self-consistency in density functional calculations Phys. Rev. B 64, 121101(R).","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"Parameters of the calculation. Since this surface is far from easy to converge, we made the problem simpler by choosing a smaller Ecut and smaller values for n_GaAs and n_vacuum. More interesting settings are Ecut = 15 and n_GaAs = n_vacuum = 20.","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"miller = (1, 1, 0) # Surface Miller indices\nn_GaAs = 2 # Number of GaAs layers\nn_vacuum = 4 # Number of vacuum layers\nEcut = 5 # Hartree\nkgrid = (4, 4, 1); # Monkhorst-Pack mesh\nnothing #hide","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"Use ASE to build the structure:","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"using PyCall\n\nase_build = pyimport(\"ase.build\")\na = 5.6537 # GaAs lattice parameter in Ångström (because ASE uses Å as length unit)\ngaas = ase_build.bulk(\"GaAs\", \"zincblende\", a=a)\nsurface = ase_build.surface(gaas, miller, n_GaAs, 0, periodic=true);\nnothing #hide","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"Get the amount of vacuum in Ångström we need to add","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"d_vacuum = maximum(maximum, surface.cell) / n_GaAs * n_vacuum\nsurface = ase_build.surface(gaas, miller, n_GaAs, d_vacuum, periodic=true);\nnothing #hide","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"Write an image of the surface and embed it as a nice illustration:","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"pyimport(\"ase.io\").write(\"surface.png\", surface * (3, 3, 1),\n rotation=\"-90x, 30y, -75z\")","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"<img src=\"../surface.png\" width=500 height=500 />","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"Use the load_atoms and load_lattice functions to convert to DFTK datastructures. These two functions not only support importing ASE atoms into DFTK, but a few more third-party datastructures as well. Typically the imported atoms use a bare Coulomb potential, such that appropriate pseudopotentials need to be attached in a post-step:","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"using DFTK\n\natoms = load_atoms(surface)\natoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional=\"pbe\")) => position\n for (el, position) in atoms]\nlattice = load_lattice(surface);\nnothing #hide","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"We model this surface with (quite large a) temperature of 0.01 Hartree to ease convergence. Try lowering the SCF convergence tolerance (tol or the temperature to see the full challenge of this system.","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"model = model_DFT(lattice, atoms, [:gga_x_pbe, :gga_c_pbe],\n temperature=0.001, smearing=DFTK.Smearing.Gaussian())\nbasis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\n\nscfres = self_consistent_field(basis, tol=1e-4, mixing=KerkerMixing());\nnothing #hide","category":"page"},{"location":"examples/ase/","page":"Creating slabs with ASE","title":"Creating slabs with ASE","text":"scfres.energies","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/docs/src/guide/tutorial.jl\"","category":"page"},{"location":"guide/tutorial/#Tutorial","page":"Tutorial","title":"Tutorial","text":"","category":"section"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"(Image: ) (Image: )","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"This document provides an overview of the structure of the code and how to access basic information about calculations. Basic familiarity with the concepts of plane-wave density functional theory is assumed throughout. Feel free to take a look at the density-functional theory chapter for some introductory lectures and introductory material on the topic.","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"note: Convergence parameters in the documentation\nWe use rough parameters in order to be able to automatically generate this documentation very quickly. Therefore results are far from converged. Tighter thresholds and larger grids should be used for more realistic results.","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"For our discussion we will use the classic example of computing the LDA ground state of the silicon crystal. Performing such a calculation roughly proceeds in three steps.","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"using DFTK\nusing Plots\n\n# 1. Define lattice and atomic positions\na = 10.26 # Silicon lattice constant in Bohr\nlattice = a / 2 * [[0 1 1.]; # Silicon lattice vectors\n [1 0 1.]; # specified column by column\n [1 1 0.]]\n\n# Load HGH pseudopotential for Silicon\nSi = ElementPsp(:Si, psp=load_psp(\"hgh/lda/Si-q4\"))\n\n# Specify type and positions of atoms\natoms = [Si => [ones(3)/8, -ones(3)/8]]\n\n# 2. Select model and basis\nmodel = model_LDA(lattice, atoms)\nkgrid = [4, 4, 4] # k-point grid (Regular Monkhorst-Pack grid)\nEcut = 7 # kinetic energy cutoff in Hartree\nbasis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\n\n# 3. Run the SCF procedure to obtain the ground state\nscfres = self_consistent_field(basis, tol=1e-8);\nnothing #hide","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"That's it! Now you can get various quantities from the result of the SCF. For instance, the different components of the energy:","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"scfres.energies","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"Eigenvalues:","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"hcat(scfres.eigenvalues...)","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"eigenvalues is an array (indexed by kpoints) of arrays (indexed by eigenvalue number). The \"splatting\" operation ... calls hcat with all the inner arrays as arguments, which collects them into a matrix.","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"The resulting matrix is 7 (number of computed eigenvalues) by 8 (number of kpoints). There are 7 eigenvalues per kpoint because there are 4 occupied states in the system (4 valence electrons per silicon atom, two atoms per unit cell, and paired spins), and the eigensolver gives itself some breathing room by computing some extra states (see n_ep_extra argument to self_consistent_field).","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"We can check the occupations:","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"hcat(scfres.occupation...)","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"And density:","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"rvecs = collect(r_vectors(basis))[:, 1, 1] # slice along the x axis\nx = [r[1] for r in rvecs] # only keep the x coordinate\nplot(x, scfres.ρ.real[:, 1, 1], label=\"\", xlabel=\"x\", ylabel=\"ρ\", marker=2)","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"We can also perform various postprocessing steps: for instance compute a band structure","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"plot_bandstructure(scfres, kline_density=5, unit=:eV)","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"or forces","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"forces(scfres)[1] # Select silicon forces","category":"page"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"The [1] extracts the forces for the first kind of atoms, i.e. Si (silicon) in the setup of the atoms list of step 1 above. As expected, they are almost zero in this highly symmetric configuration.","category":"page"},{"location":"guide/tutorial/#Where-to-go-from-here","page":"Tutorial","title":"Where to go from here","text":"","category":"section"},{"location":"guide/tutorial/","page":"Tutorial","title":"Tutorial","text":"Take a look at the example index to continue exploring DFTK.","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/polarizability.jl\"","category":"page"},{"location":"examples/polarizability/#Polarizability-by-linear-response","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"","category":"section"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"We compute the polarizability of a Helium atom. The polarizability is defined as the change in dipole moment","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"mu = int r ρ(r) dr","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"with respect to a small uniform electric field E = -x.","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"We compute this in two ways: first by finite differences (applying a finite electric field), then by linear response. Note that DFTK is not really adapted to isolated atoms because it uses periodic boundary conditions. Nevertheless we can simply embed the Helium atom in a large enough box (although this is computationally wasteful).","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"As in other tests, this is not fully converged, convergence parameters were simply selected for fast execution on CI,","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"using DFTK\nusing LinearAlgebra\n\na = 10.\nlattice = a * I(3) # cube of ``a`` bohrs\nHe = ElementPsp(:He, psp=load_psp(\"hgh/lda/He-q2\"))\natoms = [He => [[1/2; 1/2; 1/2]]] # Helium at the center of the box\n\nkgrid = [1, 1, 1] # no kpoint sampling for an isolated system\nEcut = 30\ntol = 1e-8\n\n# dipole moment of a given density (assuming the current geometry)\nfunction dipole(ρ)\n basis = ρ.basis\n rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]\n d = sum(rr .* ρ.real) * basis.model.unit_cell_volume / prod(basis.fft_size)\nend;\nnothing #hide","category":"page"},{"location":"examples/polarizability/#Polarizability-by-finite-differences","page":"Polarizability by linear response","title":"Polarizability by finite differences","text":"","category":"section"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"We first compute the polarizability by finite differences. First compute the dipole moment at rest:","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"model = model_LDA(lattice, atoms; symmetries=false)\nbasis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\nres = self_consistent_field(basis, tol=tol)\nμref = dipole(res.ρ)","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"Then in a small uniform field:","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"ε = .01\nmodel_ε = model_LDA(lattice, atoms; extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],\n symmetries=false)\nbasis_ε = PlaneWaveBasis(model_ε, Ecut; kgrid=kgrid)\nres_ε = self_consistent_field(basis_ε, tol=tol)\nμε = dipole(res_ε.ρ)","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"polarizability = (με - μref) / ε\n\nprintln(\"Reference dipole: $μref\")\nprintln(\"Displaced dipole: $με\")\nprintln(\"Polarizability : $polarizability\")","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"The result on more converged grids is very close to published results. For example DOI 10.1039/C8CP03569E quotes 1.65 with LSDA and 1.38 with CCSD(T).","category":"page"},{"location":"examples/polarizability/#Polarizability-by-linear-response-2","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"","category":"section"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"Now we use linear response to compute this analytically; we refer to standard textbooks for the formalism. In the following, chi_0 is the independent-particle polarizability, and K the Hartree-exchange-correlation kernel. We denote with dV_rm ext an external perturbing potential (like in this case the uniform electric field). Then:","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"drho = chi_0 dV = chi_0 (dV_rm ext + K drho)","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"which implies","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"drho = (1-chi_0 K)^-1 chi_0 dV_rm ext","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"From this we identify the polarizability operator to be chi = (1-chi_0 K)^-1 chi_0. Numerically, we apply chi to dV = -x by solving a linear equation (the Dyson equation) iteratively.","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"using KrylovKit\n\n# KrylovKit cannot deal with the density as a 3D array, so we need `vec` to vectorize\n# and `devec` to \"unvectorize\"\ndevec(arr) = from_real(basis, reshape(arr, size(res.ρ.real)))\n\n# Apply (1- χ0 K) to a vectorized dρ\nfunction dielectric_operator(dρ)\n dρ = devec(dρ)\n dv = apply_kernel(basis, dρ; ρ=res.ρ)\n χ0dv = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, dv...)[1]\n vec((dρ - χ0dv).real)\nend\n\n# dVext is the potential from a uniform field interacting with the dielectric dipole\n# of the density.\ndVext = from_real(basis, [-a * (r[1] - 1/2) for r in r_vectors(basis)])\n\n# Apply χ0 once to get non-interacting dipole\ndρ_nointeract = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, dVext)[1]\n\n# Solve Dyson equation to get interacting dipole\ndρ = devec(linsolve(dielectric_operator, vec(dρ_nointeract.real), verbosity=3)[1])\n\nprintln(\"Non-interacting polarizability: $(dipole(dρ_nointeract))\")\nprintln(\"Interacting polarizability: $(dipole(dρ))\")","category":"page"},{"location":"examples/polarizability/","page":"Polarizability by linear response","title":"Polarizability by linear response","text":"As expected, the interacting polarizability matches the finite difference result. The non-interacting polarizability is higher.","category":"page"},{"location":"advanced/data_structures/#Data-structures","page":"Data structures","title":"Data structures","text":"","category":"section"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"using DFTK\na = 10.26 # Silicon lattice constant in Bohr\nlattice = a / 2 * [[0 1 1.];\n [1 0 1.];\n [1 1 0.]]\nSi = ElementPsp(:Si, psp=load_psp(\"hgh/lda/Si-q4\"))\natoms = [Si => [ones(3)/8, -ones(3)/8]]\n\nmodel = model_LDA(lattice, atoms)\nkgrid = [4, 4, 4]\nEcut = 15\nbasis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\nscfres = self_consistent_field(basis, tol=1e-8);","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"In this section we assume a calculation of silicon LDA model in the setup described in Tutorial.","category":"page"},{"location":"advanced/data_structures/#Model-datastructure","page":"Data structures","title":"Model datastructure","text":"","category":"section"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"The physical model to be solved is defined by the Model datastructure. It contains the unit cell, number of electrons, atoms, type of spin polarization and temperature. Each atom has an atomic type (Element) specifying the number of valence electrons and the potential (or pseudopotential) it creates with respect to the electrons. The Model structure also contains the list of energy terms defining the energy functional to be minimised during the SCF. For the silicon example above, we used an LDA model, which consists of the following terms[2]:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"[2]: If you are not familiar with Julia syntax, typeof.(model.term_types) is equivalent to [typeof(t) for t in model.term_types].","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"typeof.(model.term_types)","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"DFTK computes energies for all terms of the model individually, which are available in scfres.energies:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"scfres.energies","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"For now the following energy terms are available in DFTK:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"Kinetic energy\nLocal potential energy, either given by analytic potentials or specified by the type of atoms.\nNonlocal potential energy, for norm-conserving pseudopotentials\nNuclei energies (Ewald or pseudopotential correction)\nHartree energy\nExchange-correlation energy\nPower nonlinearities (useful for Gross-Pitaevskii type models)\nMagnetic field energy\nEntropy term","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"Custom types can be added if needed. For examples see the definition of the above terms in the src/terms directory.","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"By mixing and matching these terms, the user can create custom models not limited to DFT. Convenience constructors are provided for common cases:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"model_LDA: LDA model using the Teter parametrisation\nmodel_DFT: Assemble a DFT model using any of the LDA or GGA functionals of the libxc library, for example: model_DFT(lattice, atoms, [:gga_x_pbe, :gga_c_pbe]) model_DFT(lattice, atoms, :lda_xc_teter93) where the latter is equivalent to model_LDA. Specifying no functional is the reduced Hartree-Fock model: model_DFT(lattice, atoms, [])\nmodel_atomic: A linear model, which contains no electron-electron interaction (neither Hartree nor XC term).","category":"page"},{"location":"advanced/data_structures/#PlaneWaveBasis-and-plane-wave-discretisations","page":"Data structures","title":"PlaneWaveBasis and plane-wave discretisations","text":"","category":"section"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"The PlaneWaveBasis datastructure handles the discretization of a given Model in a plane-wave basis. In plane-wave methods the discretization is twofold: Once the k-point grid, which determines the sampling inside the Brillouin zone and on top of that a finite plane-wave grid to discretise the lattice-periodic functions. The former aspect is controlled by the kgrid argument of PlaneWaveBasis, the latter is controlled by the cutoff energy parameter Ecut:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"PlaneWaveBasis(model, Ecut; kgrid=kgrid)","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"The PlaneWaveBasis by default uses symmetry to reduce the number of k-points explicitly treated. For details see Crystal symmetries.","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"As mentioned, the periodic parts of Bloch waves are expanded in a set of normalized plane waves e_G:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"beginaligned\n psi_k(x) = e^i k cdot x u_k(x)\n = sum_G in mathcal R^* c_G e^i k cdot x e_G(x)\nendaligned","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"where mathcal R^* is the set of reciprocal lattice vectors. The c_G are ell^2-normalized. The summation is truncated to a \"spherical\", k-dependent basis set","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":" S_k = leftG in mathcal R^* middle\n frac 1 2 k+ G^2 le E_textcutright","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"where E_textcut is the cutoff energy.","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"Densities involve terms like psi_k^2 = u_k^2 and therefore products e_-G e_G for G G in S_k. To represent these we use a \"cubic\", k-independent basis set large enough to contain the set G-G G G in S_k. We can obtain the coefficients of densities on the e_G basis by a convolution, which can be performed efficiently with FFTs (see G_to_r and r_to_G functions). Potentials are discretized on this same set.","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"The normalization conventions used in the code is that quantities stored in reciprocal space are coefficients in the e_G basis, and quantities stored in real space use real physical values. This means for instance that wavefunctions in the real space grid are normalized as fracOmegaN sum_r psi(r)^2 = 1 where N is the number of grid points.","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"For example let us check the normalization of the first eigenfunction at the first k-point in reciprocal space:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"ψtest = scfres.ψ[1][:, 1]\nsum(abs2.(ψtest))","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"We now perform an IFFT to get ψ in real space. The k-point has to be passed because ψ is expressed on the k-dependent basis. Again the function is normalised:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"ψreal = G_to_r(basis, basis.kpoints[1], ψtest)\nsum(abs2.(ψreal)) * model.unit_cell_volume / prod(basis.fft_size)","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"The list of k points of the basis can be obtained with basis.kpoints.","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"basis.kpoints","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"The G vectors of the \"spherical\", k-dependent grid can be obtained with G_vectors:","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"[length(G_vectors(kpoint)) for kpoint in basis.kpoints]","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"ik = 1\nG_vectors(basis.kpoints[ik])[1:4]","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"The list of G vectors (Fourier modes) of the \"cubic\", k-independent basis set can be obtained similarly with G_vectors(basis).","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"length(G_vectors(basis)), prod(basis.fft_size)","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"collect(G_vectors(basis))[1:4]","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"Analogously the list of r vectors (real-space grid) can be obtained with r_vectors(basis):","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"length(r_vectors(basis))","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"collect(r_vectors(basis))[1:4]","category":"page"},{"location":"advanced/data_structures/#Accessing-Bloch-waves-and-densities","page":"Data structures","title":"Accessing Bloch waves and densities","text":"","category":"section"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"Wavefunctions are stored in an array scfres.ψ as ψ[ik][iG, iband] where ik is the index of the kpoint (in basis.kpoints), iG is the index of the plane wave (in G_vectors(basis.kpoints[ik])) and iband is the index of the band. Densities are usually stored in a special type, RealFourierArray, from which the representation in real and reciprocal space can be accessed using ρ.real and ρ.fourier respectively.","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"using Plots # hide\nrvecs = collect(r_vectors(basis))[:, 1, 1] # slice along the x axis\nx = [r[1] for r in rvecs] # only keep the x coordinate\nplot(x, scfres.ρ.real[:, 1, 1], label=\"\", xlabel=\"x\", ylabel=\"ρ\", marker=2)","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"G_energies = [sum(abs2.(model.recip_lattice * G)) ./ 2 for G in G_vectors(basis)][:]\nscatter(G_energies, abs.(scfres.ρ.fourier[:]);\n yscale=:log10, ylims=(1e-12, 1), label=\"\", xlabel=\"Energy\", ylabel=\"|ρ|^2\")","category":"page"},{"location":"advanced/data_structures/","page":"Data structures","title":"Data structures","text":"Note that the density has no components on wavevectors above a certain energy, because the wavefunctions are limited to frac 1 2k+G^2 E_rm cut.","category":"page"},{"location":"guide/density_functional_theory/#density-functional-theory","page":"Density-functional theory","title":"Introductory resources on density-functional theory and DFTK","text":"","category":"section"},{"location":"guide/density_functional_theory/","page":"Density-functional theory","title":"Density-functional theory","text":"This page collects a bunch of lecture notes, textbooks and recordings related to density-functional theory (DFT) and DFTK. Since DFTK aims for an interdisciplinary audience the level and scope of the referenced works varies. They are roughly ordered from beginner to advanced. For a list of articles dealing with novel research aspects achieved using DFTK, see Publications.","category":"page"},{"location":"guide/density_functional_theory/","page":"Density-functional theory","title":"Density-functional theory","text":"DFTK: A Julian approach for simulating electrons in solids by M. F. Herbst: Pre-recorded talk for JuliaCon 2020. Assumes no knowledge about DFT and gives the broad picture of DFTK. Slides.\nIntroduction to plane-wave DFT and DFTK by M. F. Herbst: Jupyter notebooks a two-hour lecture introducing plane-wave DFT and DFTK to an interdisciplinary audience of undergrads. You can also run the notebooks interactively on binder.\nElectronic Structure: Basic theory and practical methods by R. M. Martin (Cambridge University Press, 2004): Standard textbook introducing most common methods of the field (lattices, pseudopotentials, DFT, ...) from the perspective of a physicist\nA Mathematical Introduction to Electronic Structure Theory by L. Lin and J. Lu (SIAM, 2019): Monograph attacking DFT from a mathematical angle. Covers topics such as DFT, pseudos, SCF, response, ...","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/scf_checkpoints.jl\"","category":"page"},{"location":"examples/scf_checkpoints/#Saving-SCF-results-on-disk-and-SCF-checkpoints","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"","category":"section"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"For longer DFT calculations it is pretty standard to run them on a cluster in advance and to perform postprocessing (band structure calculation, plotting of density, etc.) at a later point and potentially on a different machine.","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"To support such workflows DFTK offers the two functions save_scfres and load_scfres, which allow to save the data structure returned by self_consistent_field on disk or retrieve it back into memory, respectively. For this purpose DFTK uses the JLD2.jl file format and Julia package. For the moment this process is considered an experimental feature and has a number of caveats, see the warnings below.","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"warning: Saving `scfres` is experimental\nThe load_scfres and save_scfres pair of functions are experimental features. This means:The interface of these functions as well as the format in which the data is stored on disk can change incompatibly in the future. At this point we make no promises ...\nJLD2 is not yet completely matured and it is recommended to only use it for short-term storage and not to archive scientific results.\nIf you are using the functions to transfer data between different machines ensure that you use the same version of Julia, JLD2 and DFTK for saving and loading data.","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"To illustrate the use of the functions in practice we will compute the total energy of the O₂ molecule at PBE level. To get the triplet ground state we use a collinear spin polarisation (see Collinear spin and magnetic systems for details) and a bit of temperature to ease convergence:","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"using DFTK\nusing LinearAlgebra\nusing JLD2\n\nd = 2.079 # oxygen-oxygen bondlength\na = 9.0 # size of the simulation box\nlattice = diagm(a * ones(3))\nO = ElementPsp(:O, psp=load_psp(\"hgh/pbe/O-q6.hgh\"))\natoms = [O => d / 2a * [[0, 0, 1], [0, 0, -1]]]\nmagnetic_moments = [O => [1., 1.]]\n\nEcut = 10 # Far too small to be converged\nmodel = model_PBE(lattice, atoms, temperature=0.02, smearing=smearing=Smearing.Gaussian(),\n magnetic_moments=magnetic_moments)\nbasis = PlaneWaveBasis(model, Ecut; kgrid=[1, 1, 1])\n\nρspin = guess_spin_density(basis, magnetic_moments)\nscfres = self_consistent_field(basis, tol=1e-2, ρspin=ρspin)\nsave_scfres(\"scfres.jld2\", scfres);\nnothing #hide","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"scfres.energies","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"The scfres.jld2 file could now be transfered to a different computer, Where one could fire up a REPL to inspect the results of the above calculation:","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"using DFTK\nusing JLD2\nloaded = load_scfres(\"scfres.jld2\")\npropertynames(loaded)","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"loaded.energies","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"Since the loaded data contains exactly the same data as the scfres returned by the SCF calculation one could use it to plot a band structure, e.g. plot_bandstructure(load_scfres(\"scfres.jld2\")) directly from the stored data.","category":"page"},{"location":"examples/scf_checkpoints/#Checkpointing-of-SCF-calculations","page":"Saving SCF results on disk and SCF checkpoints","title":"Checkpointing of SCF calculations","text":"","category":"section"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"A related feature, which is very useful especially for longer calculations with DFTK is automatic checkpointing, where the state of the SCF is periodically written to disk. The advantage is that in case the calculation errors or gets aborted due to overrunning the walltime limit one does not need to start from scratch, but can continue the calculation from the last checkpoint.","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"To enable automatic checkpointing in DFTK one needs to pass the ScfSaveCheckpoints callback to self_consistent_field, for example:","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"callback = DFTK.ScfSaveCheckpoints()\nscfres = self_consistent_field(basis, tol=1e-2, ρspin=ρspin, callback=callback);\nnothing #hide","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"Notice that using this callback makes the SCF go silent since the passed callback parameter overwrites the default value (namely DefaultScfCallback()) which exactly gives the familiar printing of the SCF convergence. If you want to have both (printing and checkpointing) you need to chain both callbacks:","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"callback = DFTK.ScfDefaultCallback() ∘ DFTK.ScfSaveCheckpoints(keep=true)\nscfres = self_consistent_field(basis, tol=1e-2, ρspin=ρspin, callback=callback);\nnothing #hide","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"For more details on using callbacks with DFTK's self_consistent_field function see Monitoring self-consistent field calculations.","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"By default checkpoint is saved in the file dftk_scf_checkpoint.jld2, which is deleted automatically once the SCF completes successfully. If one wants to keep the file one needs to specify keep=true as has been done in the ultimate SCF for demonstration purposes: now we can continue the previous calculation from the last checkpoint as if the SCF had been aborted. For this one just loads the checkpoint with load_scfres:","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"oldstate = load_scfres(\"dftk_scf_checkpoint.jld2\")\nscfres = self_consistent_field(oldstate.basis, ρ=oldstate.ρ, ρspin=oldstate.ρspin,\n ψ=oldstate.ψ, tol=1e-3);\nnothing #hide","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"note: Availability of `load_scfres`, `save_scfres` and `ScfSaveCheckpoints`\nAs JLD2 is an optional dependency of DFTK these three functions are only available once one has both imported DFTK and JLD2 (using DFTK and using JLD2).","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"(Cleanup files generated by this notebook)","category":"page"},{"location":"examples/scf_checkpoints/","page":"Saving SCF results on disk and SCF checkpoints","title":"Saving SCF results on disk and SCF checkpoints","text":"rm(\"dftk_scf_checkpoint.jld2\")\nrm(\"scfres.jld2\")","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/cohen_bergstresser.jl\"","category":"page"},{"location":"examples/cohen_bergstresser/#Cohen-Bergstresser-model","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"","category":"section"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"This example considers the Cohen-Bergstresser model[CB1966], reproducing the results of the original paper. This model is particularly simple since its linear nature allows one to get away without any self-consistent field calculation.","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"[CB1966]: M. L. Cohen and T. K. Bergstresser Phys. Rev. 141, 789 (1966) DOI 10.1103/PhysRev.141.789","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"We build the lattice using the tabulated lattice constant from the original paper, stored in DFTK:","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"using DFTK\n\nSi = ElementCohenBergstresser(:Si)\nlattice = Si.lattice_constant / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]\natoms = [Si => [ones(3)/8, -ones(3)/8]];\nnothing #hide","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"Next we build the rather simple model and discretise it with moderate Ecut:","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"Ecut = 10.0\nmodel = Model(lattice; atoms=atoms, terms=[Kinetic(), AtomicLocal()])\nbasis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1));\nnothing #hide","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"We diagonalise at the Gamma point to find a Fermi level ...","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"ham = Hamiltonian(basis)\neigres = diagonalize_all_kblocks(DFTK.lobpcg_hyper, ham, 6)\nεF = fermi_level(basis, eigres.λ)","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"... and compute and plot 8 bands:","category":"page"},{"location":"examples/cohen_bergstresser/","page":"Cohen-Bergstresser model","title":"Cohen-Bergstresser model","text":"using Plots\n\nn_bands = 8\nρ0 = guess_density(basis) # Just dummy, has no meaning in this model\nρspin0 = nothing\np = plot_bandstructure(basis, ρ0, ρspin0, n_bands, εF=εF, kline_density=10)\nylims!(p, (-5, 6))","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/custom_potential.jl\"","category":"page"},{"location":"examples/custom_potential/#Custom-potential","page":"Custom potential","title":"Custom potential","text":"","category":"section"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"We solve the 1D Gross-Pitaevskii equation with a custom potential. This is similar to Gross-Pitaevskii equation in one dimension and we show how to define local potentials attached to atoms, which allows for instance to compute forces.","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"using DFTK\nusing LinearAlgebra","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"First, we define a new element which represents a nucleus generating a custom potential","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"struct ElementCustomPotential <: DFTK.Element\n pot_real::Function # Real potential\n pot_fourier::Function # Fourier potential\nend","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"We need to extend two methods to access the real and Fourier forms of the potential during the computations performed by DFTK","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"function DFTK.local_potential_fourier(el::ElementCustomPotential, q::Real)\n return el.pot_fourier(q)\nend\nfunction DFTK.local_potential_real(el::ElementCustomPotential, r::Real)\n return el.pot_real(r)\nend","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"We set up the lattice. For a 1D case we supply two zero lattice vectors","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"a = 10\nlattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];\nnothing #hide","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"In this example, we want to generate two Gaussian potentials generated by two nuclei localized at positions x_1 and x_2, that are expressed in 01) in fractional coordinates. x_1 - x_2 should be different from 05 to break symmetry and get nonzero forces.","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"x1 = 0.2\nx2 = 0.8;\nnothing #hide","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"We define the width of the Gaussian potential generated by one nucleus","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"L = 0.5;\nnothing #hide","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"We set the potential in its real and Fourier forms","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"pot_real(x) = exp(-(x/L)^2)\npot_fourier(q::T) where {T <: Real} = exp(- (q*L)^2 / 4);\nnothing #hide","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"And finally we build the elements and set their positions in the atoms array. Note that in this example pot_real is not required as all applications of local potentials are done in the Fourier space.","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"nucleus = ElementCustomPotential(pot_real, pot_fourier)\natoms = [nucleus => [x1*[1,0,0], x2*[1,0,0]]];\nnothing #hide","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"Setup the Gross-Pitaevskii model","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"C = 1.0\nα = 2;\nn_electrons = 1 # Increase this for fun\nterms = [Kinetic(),\n AtomicLocal(),\n PowerNonlinearity(C, α),\n]\nmodel = Model(lattice; atoms=atoms, n_electrons=n_electrons, terms=terms,\n spin_polarization=:spinless); # use \"spinless electrons\"\nnothing #hide","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"We discretize using a moderate Ecut and run a SCF algorithm to compute forces afterwards. As there is no ionic charge associated to nucleus we have to specify a starting density and we choose to start from a zero density.","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"Ecut = 500\nbasis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))\nρ = zeros(complex(eltype(basis)), basis.fft_size)\nscfres = self_consistent_field(basis, tol=1e-8, ρ=from_fourier(basis, ρ))\nscfres.energies","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"Computing the forces can then be done as usual:","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"hcat(forces(scfres)...)","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"Extract the converged total local potential","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1\nnothing #hide","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"Extract other quantities before plotting them","category":"page"},{"location":"examples/custom_potential/","page":"Custom potential","title":"Custom potential","text":"ρ = real(scfres.ρ.real)[:, 1, 1] # converged density\nψ_fourier = scfres.ψ[1][:, 1]; # first kpoint, all G components, first eigenvector\nψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\nψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));\n\nusing Plots\nx = a * vec(first.(DFTK.r_vectors(basis)))\np = plot(x, real.(ψ), label=\"real(ψ)\")\nplot!(p, x, imag.(ψ), label=\"imag(ψ)\")\nplot!(p, x, ρ, label=\"ρ\")\nplot!(p, x, tot_local_pot, label=\"tot local pot\")","category":"page"},{"location":"#DFTK.jl:-The-density-functional-toolkit.","page":"Home","title":"DFTK.jl: The density-functional toolkit.","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"The density-functional toolkit, DFTK for short, is a library of Julia routines for playing with plane-wave density-functional theory (DFT) algorithms. In its basic formulation it solves periodic Kohn-Sham equations. The unique feature of the code is its emphasis on simplicity and flexibility with the goal of facilitating methodological development and interdisciplinary collaboration. In about 5k lines of pure Julia code we already support a sizeable set of features, after just a good year of development. Our performance is of the same order of magnitude as much larger production codes such as Abinit, Quantum Espresso and VASP.","category":"page"},{"location":"#package-features","page":"Home","title":"Package features","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Methods and models:\nKohn-Sham-like models, with an emphasis on flexibility: compose your own model, from Cohen-Bergstresser linear eigenvalue equations to Gross-Pitaevskii equations and sophisticated LDA/GGA functionals (any functional from the libxc library)\nAnalytic potentials or Godecker norm-conserving pseudopotentials (GTH, HGH)\nBrillouin zone symmetry for k-Point sampling using spglib\nSmearing functions for metals\nCollinear spin polarization for magnetic systems\nSelf-consistent field approaches: Damping, Kerker mixing, Anderson/Pulay/DIIS mixing\nDirect minimization\nMulti-level threading (kpoints, eigenvectors, FFTs, linear algebra)\n1D / 2D / 3D systems\nMagnetic fields\nTreat systems beyond 500 electrons\nGround-state properties and post-processing:\nTotal energy\nForces\nDensity of states (DOS), local density of states (LDOS)\nBand structures\nEasy access to all intermediate quantities (e.g. density, Bloch waves)\nSupport for arbitrary floating point types, including Float32 (single precision) or Double64 (from DoubleFloats.jl). For DFT this is currently restricted to LDA (with Slater exchange and VWN correlation).\nThird-party integrations:\nUse structures prepared in pymatgen, ASE or abipy.\nasedftk: DFTK-based calculator implementation for ASE.\nRead data in ETSF Nanoquanta format.","category":"page"},{"location":"","page":"Home","title":"Home","text":"Missing a feature? Look for an open issue or create a new one. Want to contribute? See our contributing notes.","category":"page"},{"location":"#example-index","page":"Home","title":"Example index","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"First, new users should take a look at the Installation and Tutorial sections. More details about DFTK are explained in the examples as we go along:","category":"page"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"examples/metallic_systems.md\",\n \"examples/pymatgen.md\",\n \"examples/ase.md\",\n \"examples/collinear_magnetism.md\",\n \"examples/geometry_optimization.md\",\n \"examples/scf_callbacks.md\",\n \"examples/scf_checkpoints.md\",\n \"examples/polarizability.md\",\n \"examples/gross_pitaevskii.md\",\n \"examples/gross_pitaevskii_2D.md\",\n \"examples/cohen_bergstresser.md\",\n \"examples/arbitrary_floattype.md\",\n \"examples/custom_solvers.md\",\n \"examples/custom_potential.md\",\n]\nDepth = 1","category":"page"},{"location":"","page":"Home","title":"Home","text":"These and more examples can be found in the examples directory of the main code.","category":"page"},{"location":"","page":"Home","title":"Home","text":"If you have a great example you think would fit here, please open a pull request!","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"EditURL = \"https://github.com/JuliaMolSim/DFTK.jl/blob/master/examples/geometry_optimization.jl\"","category":"page"},{"location":"examples/geometry_optimization/#Geometry-optimization","page":"Geometry optimization","title":"Geometry optimization","text":"","category":"section"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"(Image: ) (Image: )","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"We use the DFTK and Optim packages in this example to find the minimal-energy bond length of the H_2 molecule. We setup H_2 in an LDA model just like in the Tutorial for silicon.","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"using DFTK\nusing Optim\nusing LinearAlgebra\nusing Printf\n\nkgrid = [1, 1, 1] # k-Point grid\nEcut = 5 # kinetic energy cutoff in Hartree\ntol = 1e-8 # tolerance for the optimization routine\na = 10 # lattice constant in Bohr\nlattice = a * Diagonal(ones(3))\nH = ElementPsp(:H, psp=load_psp(\"hgh/lda/h-q1\"));\nnothing #hide","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"We define a blochwave and a density to be used as global variables so that we can transfer the solution from one iteration to another and therefore reduce the optimization time.","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"ψ = nothing\nρ = nothing","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"First, we create a function that computes the solution associated to the position x in mathbbR^6 of the atoms in reduced coordinates (cf. Reduced and cartesian coordinates for more details on the coordinates system). They are stored as a vector: x[1:3] represents the position of the first atom and x[4:6] the position of the second. We also update ψ and ρ for the next iteration.","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"function compute_scfres(x)\n atoms = [H => [x[1:3], x[4:6]]]\n model = model_LDA(lattice, atoms)\n basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)\n global ψ, ρ\n if ρ === nothing\n ρ = guess_density(basis)\n end\n scfres = self_consistent_field(basis; ψ=ψ, ρ=ρ,\n tol=tol / 10, callback=info->nothing)\n ψ = scfres.ψ\n ρ = scfres.ρ\n scfres\nend;\nnothing #hide","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"Then, we create the function we want to optimize: fg! is used to update the value of the objective function F, namely the energy, and its gradient G, here computed with the forces (which are, by definition, the negative gradient of the energy).","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"function fg!(F, G, x)\n scfres = compute_scfres(x)\n if G != nothing\n grad = forces(scfres)\n G .= -[grad[1][1]; grad[1][2]]\n end\n scfres.energies.total\nend;\nnothing #hide","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"Now, we can optimize on the 6 parameters x = [x1, y1, z1, x2, y2, z2] in reduced coordinates, using LBFGS(), the default minimization algorithm in Optim. We start from x0, which is a first guess for the coordinates. By default, optimize traces the output of the optimization algorithm during the iterations. Once we have the minimizer xmin, we compute the bond length in cartesian coordinates.","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"x0 = vcat(lattice \\ [0., 0., 0.], lattice \\ [1.4, 0., 0.])\nxres = optimize(Optim.only_fg!(fg!), x0, LBFGS(),\n Optim.Options(show_trace=true, f_tol=tol))\nxmin = Optim.minimizer(xres)\ndmin = norm(lattice*xmin[1:3] - lattice*xmin[4:6])\n@printf \"\\nOptimal bond length for Ecut=%.2f: %.3f Bohr\\n\" Ecut dmin","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"We used here very rough parameters to generate the example and setting Ecut to 10 Ha yields a bond length of 1.523 Bohr, which agrees with ABINIT.","category":"page"},{"location":"examples/geometry_optimization/","page":"Geometry optimization","title":"Geometry optimization","text":"note: Degrees of freedom\nWe used here a very general setting where we optimized on the 6 variables representing the position of the 2 atoms and it can be easily extended to molecules with more atoms (such as H_2O). In the particular case of H_2, we could use only the internal degree of freedom which, in this case, is just the bond length.","category":"page"}]
}