Multiple GW Solutions: Limitations of GW methods

In this part we will discuss how to get deeper insight into the solution behavior of GW and address cases with multiple solutions. In most cases, a unique quasiparticle (QP) solution (\(\epsilon_i^{QP}\)) exists. However, if multiple solutions occur, iterating the QP equation,

\[ \epsilon_i^{QP} = \epsilon_i^{KS} + \mathrm{Re} \Sigma^{\mathrm c}_i(\epsilon_i^{QP}) + \Sigma^{\mathrm x}_i - V^{xc}_i \quad, \]

will still yield only one solution without any information on its spectral weight. Multiple solutions are typically a sign that the \(GW\) level of theory is insufficient for the system under investigation and, e.g., including some level of self-consistency might be required. A possible indication for a multiple solution behavior is that the iteration of the QP equations does not converge or that different \(GW\) codes yield different results, even though the same numerical approximations are used.

We can obtain all possible solutions by solving the QP equation graphically. The graphical solution is obtained by plotting the real part of the correlation part of the self-energy matrix elements \(\Sigma^c_i(\omega)\) over a frequency range of several eV around the expected solution. The following figure shows the self-energy matrix elements for the HOMO of a single water molecule.

graphical_solution_figure Figure adapted from Golze et al. ¹

The intersections with the red straight line \(\omega + V^{xc}_i -\Sigma^x_i - \epsilon_i^{KS}\) are solutions of the QP equation. An intersection with a small slope of \(\mathrm{Re}\Sigma_i^c\) indicates that the solution has a large spectral weight, whereas intersections with a large slope of \(\mathrm{Re}\Sigma_i^c\) correspond to solutions with very low spectral weight (satellites); see the \(GW\) compendium ¹ for more information. In our example, we have a clear unique QP solution at around -12 eV, whereas the intersections at around -30 eV are satellite features.

Graphical solution for the C₂H₄ molecule

In this part of the exercise, we want to evaluate the QP energy for the HOMO of C\(_2\)H\(_4\) graphically at the \(G_0W_0@\)PBE level.

Calculate the self-energy matrix elements for the HOMO, which corresponds to state 8. Modify the control.in that you used for the \(G_0W_0@\)PBE calculation with the tier2 basis sets by adding

print_self_energy 8

You will obtain the file self_energy_analytic_state_8.dat, which contains the frequency in eV (1st column) and \(\mathrm{Re}\Sigma_{\mathrm{HOMO}}^c\) in Hartree (2nd column). Plot the self-energy matrix elements in the frequency range -30 to 20 eV. Don't forget to convert the self-energy also to eV.
Plot in the same figure the straight line \(\omega + V^{xc}_{\mathrm{HOMO}}-\Sigma^x_{\mathrm{HOMO}} - \epsilon_{\mathrm{HOMO}}^{KS}\). You get \(V^{xc}_{\mathrm{HOMO}}\), \(\Sigma^x_{\mathrm{HOMO}}\) and \(\epsilon_{\mathrm{HOMO}}^{KS}\) from the output file of the calculation; see Part-1.1 (GW eigenvalues and Photoemission).
Find the frequency/frequencies where the self-energy and the straight line intersect. You should find that the graphical solution agrees with the iterative solution of the QP solution, which you obtained before. You should find exactly one solution.

Multiple solution case: Boron nitride dimer

In this exercise, we will recompute a multisolution case from the famous \(GW\)100 benchmark set.² Our example is the HOMO excitation of the BN dimer computed with \(G_0W_0@\)PBE.

Create a geometry.in file with the following header

atom 0.0000 0.0000 0.0000 B
atom 0.0000 0.0000 1.281  N

Copy the control.in from the C\(_2\)H\(_4\) calculation above and add/modify the following lines

  n_anacon_par 128
  frequency_points 400
  print_self_energy 6

  relativistic none

We perform here an analytic continuation with 128 parameters (default is 16) and thus increase the number of frequency points to 400 (default is 100). The self-energy elements are printed for the HOMO (= state 6). The \(GW100\) benchmark results are entirely non-relativistic. To (exactly) reproduce the results of the benchmark, we perform a non-relativistic calculation.

Remove the tier basis sets from the control.in and use the def2-QZVP basis sets, which were used for the \(GW100\) benchmark study. You find them in the FHIaims directory under

species_defaults/non-standard/gaussian_tight_770/def2-QZVP/

Proceed as before and plot the self-energy matrix elements for the HOMO and the straight line \(\omega + V^{xc}_{\mathrm{HOMO}}-\Sigma^x_{\mathrm{HOMO}} - \epsilon_{\mathrm{HOMO}}^{KS}\) in the range from -15 to -8 eV. Find the intersections. You should obtain two solutions, one at -11.69 eV and another at -11.02 eV. Compare your results to the GW100 benchmark¹ (check also the SI).
Optional: Repeat the whole procedure, but use

qpe_calc           ev_scgw

This triggers a partially self-consistent GW scheme instead of \(G_0W_0\), also known as eigenvalue-selfconsistent GW (evGW). You should find a unique solution now at -12.03 eV. Hint: Make sure you take the GW table of the last iteration to extract \(V^{xc}_{\mathrm{HOMO}}\), \(\Sigma^x_{\mathrm{HOMO}}\) and \(\epsilon_{\mathrm{HOMO}}^{KS}\). Note that the calculation might take up to 30 min on a PC or laptop.