depolarization in thin films

[Christian Koenig]

Hi all,

I'm using the plasmon pole approximation in my GW calculations for thin films. To find a suitable value for PPAPntXp I plot the real and imaginary part of 1/epsilon_macroscopic with gnuplot, i.e. the o-...eps_q1_... file with

Re = ( ( $3) / ($2**3+$3**2) )
Im = ( (-$2) / ($2**3+$3**2) )

This is basically the negative of the electron energy loss spectrum (both components).

The local fields tutorial explains that depolarization happens if the electric field used for the q -> 0 limit in the screening is not in the plane of the film.
I can reproduce the results (yambo -o c -k hartree) as shown in the attached picture. If the electric field is perpendicular to the film and local fields are considered, then the peaks of Re and Im have a smaller amplitude.

My own calculations for thin films, e.g. for Si(111), behave somehow differently:
1) The amplitude with local fields is much smaller (Im ~ 0 + order 10e-7; Re has even smaller deviation from 1).
2) I get the same behaviour even if the electric field is in the plane.

Important things to note:
1) I have used parameters that seemed reasonable just like in the tutorial, i.e. convergence was not tested any further. The idea is to check the general behaviour.
2) A vacuum cutoff and random integration was used to account for the vacuum: yambo -o c -k hartree -r

Can someone please give me a hint want I'm doing wrong?
A plot of the results and the input file are attached.


Best,

Christian

[Daniele Varsano]

Dear Christian,

first, note that the electron energy loss spectrum is provided by the Yambo output o.eels.
Now the problem you are facing is due to the use of the coulomb cutoff potential. When computing the dielectric function
eps^-1=1+vX
here v is the cutoff potential which does not have the 1/q^2 behaviour while X behaves as q^2 in limit q->0, so you have the real part essentially equal to 1 and the Imag part with very low intensity.
This happens as the macroscopic dielectric function it is not well defined for a non-3D system and it will go to 0 for infinite volume.
What is well defined is the polarizability alfa (which is also proportional to X/q^2, as 1/eps_M i a bulk) which is also reported in output.
This is documented in the FAQ, which essentially reminds to this thread.

Hope it helps,
Daniele

[Christian Koenig]

Dear Daniele,

Thanks for your response. So that's why the polarizability is only printed out when the Coulomb cutoff is used. Then I'd best look for a peak in the polarizability to check that I'm well above the plasmon frequency. The imaginary part is smooth just like expected and shows well defined peaks but the real part isn't. That still doesn't look quite right.

Presumably the PPA still makes sense because the components of the inverse microscopic function are fitted and they still have a pole each? Can I somehow access the individual components?


Best,

Christian

[Daniele Varsano]

Dear Christian,

The imaginary part is smooth just like expected and shows well-defined peaks but the real part isn't. That still doesn't look quite right.

Well, actually the real part is really strange, it is seems something weird happened. Tomorrow I will try to check if this is reproducible in some test systems and try to understand what is going on.

Presumably the PPA still makes sense because the components of the inverse microscopic function are fitted and they still have a pole each?

I think so, but still, it is important to understand if the fitted curve make sense.

Can I somehow access the individual components?

You mean poles and residuals, or the component of the screening calculated at zero frequency and at imaginary plasmon pole frequency?
Yambo store in the ndb.pp all the components of eps^-1(w=0) and eps^-1(w=wp), next residues and poles are calculated internally.
Actually what it is stored is vX and not epw^-1=delta_gg'+1vX
You can access the content of the database by using the ncdump command, or better by using some python script reading the netcdf database.
Let me know if you need a script example,

Best,
Daniele

[Christian Koenig]

Hi,

Not sure where Davide's post went but increasing g_format_length and f_format_length didn't change things.
Without local fields the real part of alpha is also very smooth and I get a pole with the peaks of the order of 10^9.
I have tested with a different film which has the same behaviour.


Best,

Christian

[Daniele Varsano]

Hi Christian,
most probably it is a matter of precision, let me do some tests and we will be back to you soon.

Daniele

[Daniele Varsano]

Dear Christian,

please apply the attached patch and recompile:
this should solve both problem (1) precision issue in the real part; (2) large value found in the IP polarizability that was indeed a bug.

Best,
Daniele

[Christian Koenig]

Dear Daniele,

Thanks very much for the patch, both the real part and imaginary part look reasonable now.

Regarding the PPA fit I was referring to the components epsilon_G,G'^-1 (q, omega). I take it the code prints
the 0,0 component of the inverted matrix epsilon for different q as function of omega, i.e. the macroscopic dielectric
function in the files o- … .eps_q1_inv_rpa_dyson, o- … .eps_q2_inv_rpa_dyson, etc.
Therefore, the individual components epsilon_G,G'^-1 (q, omega) can not really be checked, can they?

I don't have a database ndb.pp, only ndb.dip_iR_and_P (or the static database ns.kb_pp_pwscf which shouldn't know anything about if the PPA is used).


Best,

Christian

[Daniele Varsano]

Dear Christian,

Therefore, the individual components epsilon_G,G'^-1 (q, omega) can not really be checked, can they?

Yes, they can be checked by looking inside the ndb.pp database, which is produced when doing a ppa calcualtion (yambo -p p), but not when an optic linear response (yambo -o c..) is run.

Note that what is inside the ndb.pp it is not exactly the epsilon_G,G'^-1, but vX_gg', i.e. there is a delta_gg' to be added int he diagonal components.

Best,
Daniele

[Christian Koenig]

Dear Daniele,

Sorry, it should have been obvious that I get ndb.pp only when doing the GW calculation.

I assume I get one file ndb.pp_fragment_X for each q point X. What is the data format?
In a test calulation I have 7675362 entries with two numbers which presumably are vX(w=0) and vX(w=wp) as you mentioned before. That doesn't seem to correspond to a square matrix (or half a square matrix).
From there I could probably obtain the fit parameters but I still can't compare the fit to the full function vX(w).

Why is the fit actually done for an imaginary frequency?


Best,

Christian