Example problem


This page provides a small example problem, which is not exactly electronic structure theory, but that structurally resembles the challenges we will face in density-functional theory (DFT). The tasks are deliberately formulated in a very open fashion to give you the flexibility to tackle them in a form you feel comfortable with.

A large part of the exercise is to implement a discretisation and solution procedure for the problem discussed below. Please keep an eye on writing your code in a readable and understandable fashion as I will not manage to spend a lot of time looking at it. Ideally the code should be unit tested and work beyond the scenario we are considering here, but for the sake of this exercise this is optional as you likely won't have enough time for this.

Both your code as well as your reasoning during the exercise will be discussed in the interview. Don't worry if you don't manage to complete the full exercise. The idea is that you spend a couple of hours on this, such that we have a basis for discussion and I get an idea how you approach such a problem.

It is completely natural when getting into a new topic to get stuck somewhere. I very much encourage you to send me an email with questions if this is the case.

The Gross-Pitaevskii equation

In this example we will deal with the numerical modelling of the (time-independent) Gross-Pitaevskii equation (GPE). The GPE is simpler than DFT, but since it still displays interesting behaviour, which makes it a favourite model of applied mathematicians, that study electronic structure problems. Of note the GPE itself has physical significance as well: it is a mean-field model for bosonic systems.

In one dimension the GPE reads in atomic units:

Hψ=(12Δ+V+2Cψ2)ψ=λ1ψψL2=1, H ψ = \left(-\frac12 Δ + V + 2 C |ψ|^2\right) ψ = λ_1 ψ \qquad \|ψ\|_{L^2} = 1,

where ψL2(R)\psi \in L^2(\mathbb R) is the effective one-particle bosonic wave function, CC provides the strength of the boson-boson coupling. The application of the external potential V:RRV : \mathbb R \mapsto \mathbb R is point-wise (in real space), i.e.

[Vψ](x)=V(x)ψ(x). \left[V\psi\right](x) = V(x) \psi(x).

Of note HH is an eigenvector-dependent operator. Overall the GPE is thus a non-linear eigenvalue problem, where we seek the lowest eigenvalue λ1\lambda_1 and associated eigenfunction ψ\psi.

In this example we will consider external potentials of the form

Gα,L,R(x)=α2πLexp[12(xRL)2]V(x)=Gα,L,R+Gα,L,R+βx2, \begin{align*} G_{\alpha, L, R}(x) &= -\frac{\alpha}{\sqrt{2\pi} \cdot L} \exp\left[- \frac12 \Big(\frac{x - R}{L}\Big)^2\right] \\ V(x) &= G_{\alpha, L, R} + G_{\alpha', L', R'} + \beta x^2, \end{align*}

i.e. the sum of two inverted Gaussians plus a confining parabola. Due to the confinement ψ\psi decays exponentially as x\left| x \right| \to \infty.

Task 1: Implement a numerical method for GPE problems

Consider the setting C=1C = 1, β=1/10\beta = 1/10, L=L=1L = L' = 1, α=α=5\alpha = \alpha' = 5, R=3R = -3, R=3R = 3 and find the corresponding λ1\lambda_1 and ψ\psi (lowest eigenpair of HH) numerically. You can use any discretisation method and programming language of your choice. Ensure that your method converges λ1\lambda_1 to at least 44 digits.

The remainder of this task description is basically some hints that help you to achieve this in one particularly simple way. You don't have to follow any of this and in fact there are plenty of better methods out there, which you are welcome to follow if you want.

Tip on the algorithm: Since this is a non-linear eigenvalue problem you will need another iterative procedure around the actual diagonalisation to find the eigenfunction ψ\psi. Since the phase of the eigenvectors is ill-defined, it is usually more convenient to iterate on the bosonic density ρ=ψ2\rho = |ψ|^2 instead of the eigenfunction ψ\psi. One approach we will consider here is the self-consistent field (SCF) procedure, which is a fixed-point method.

It can be implemented as follows:

Tip on the implementation: If you are unsure how to set up and discretise the problem, you can take a look at the DFTK code developed in our group. It uses plane waves and periodic boundary conditions, see this introduction in the documentation. As a result of this you will put the GPE problem above into a large box with artificial periodic boundary conditions, see this example for details. Don't forget to converge wrt. the box size. Custom potentials (like the two Gaussians) are discussed in this example. Once you have constructed a PlaneWaveBasis object, you can easily obtain a representation of the discretised Hamiltonian HH at given bosonic density ρ by constructing a Hamiltonian object, e.g.

# Make sure to use kgrid=(1,1,1) when setting up the PlaneWaveBasis
basis = PlaneWaveBasis( #= ... more arguments ... =#, kgrid=(1,1,1) )
ρ = zeros(basis.fft_size..., 1)
H = Hamiltonian(basis; ρ=ρ).blocks[1]   # [1] Select 1st (and only) k-Point

The resulting Hamiltonian object H can be directly passed to an eigensolver as it supports e.g. the * operator to perform a matrix-vector product. E.g. here we use an LOBPCG built into DFTK to find the smallest 3 eigenpairs starting from a random guess:

X    = randn(ComplexF64, size(H, 2), 3)     # ramdom guess
prec = PreconditionerTPA(H.basis, H.kpoint)  # Setup preconditioner

# Solve and get eigenpairs back
(; λ, X) = lobpcg_hyper(H, X; prec, tol=1e-6, maxiter=300)

It is generally a good idea to use the eigenvectors X from a previous SCF step as an initial guess to save on the number of required iterations of the diagonalisation procedure.

Note that DFTK has a self_consistent_field function to do a (more sophisticated) SCF. You can use it for testing, but I would still ask you to code up your own SCF as part of this exercise.

Task 2: Investigating the eigenspectrum of H

Once your numerical SCF scheme from the previous section is up and running, play with the model's parameters, in particular CC and the Gaussian parameters α\alpha, α\alpha', RR and RR'. For each parameter choice first run the SCF to obtain the solution ψ\psi and then investigate the total potential V+2Cψ2V + 2 C |ψ|^2, the density ρ=ψ2\rho = |\psi|^2 as well as the five lowest eigenvalues of the GPE operator HH. How do the eigenvalues depend on your parameter choices?

Task 3: Inverse design

For C=1C = 1 and β=1/10\beta = 1/10 find values for the Gaussian parameters such that the gaps between the first three eigenvalues have a ratio of 1:101 : 10. In other words if we denote the lowest eigenvalues of HH by λ1,λ2,λ3\lambda_1, \lambda_2, \lambda_3, then we desire

λ3λ2=10λ2λ1. |\lambda_3 - \lambda_2| = 10 |\lambda_2 - \lambda_1|.

Tip: One solution can be obtained by restricting L=L=1L = L' = 1, α=α=5\alpha = \alpha' = 5 and R=RR = -R'.

Edit this page on
Last modified: May 27, 2024. Website built with Franklin.jl and the Julia programming language.