Imag Dielectric function of gold not matching with experimental (RPA)

[muhammadhasan]

Hi,

I want to calculate dielectric function for Au (my focus on low energy ~ 1 ev). For this purpose, I have used 4*4*4 supercell, 729 k points and 800 bands in first step. I compare my results (imag dielectric function) with experimental data (https://pmc.ncbi.nlm.nih.gov/articles/PMC4875142/)

compare_with_experimental.png

and It doesn't match. I thought the problem might be in lower k-points. After that, I followed a double grid method example (https://www.attaccalite.com/speed-up-di ... ith-yambo/) and applied it in my case. I run a new nscf calculation using 4096 k points (no shift) and then p2y -w (new) and finally ypp -m (old). However, when I compared the both results, I have found double grid method shows very weird result. I have attached the figure of this plot

double_grid_Au.png

. I have also attached necessary files if you able to help me to find the problem.

double_Au.tar.gz

Thank you so much for your patience and time.

Best
Md J Hasan
PhD student
Mechanical Engineering
University of Maine

[Daniele Varsano]

Dear Hasan,

please consider you are dealing with a metal and the long wavelength limit is not included by default, as you can see from the warning in the report file. You can add it empirically, considering a Drude model (DrudeWXd). You can also have a look at threads in the forum dealing with the topic, as for instance: viewtopic.php?t=1487

Best,
Daniele

[muhammadhasan]

Hi Professor,

I have gone through some threads and specially equation 12 of one of your papers (https://doi.org/10.1103/PhysRevB.107.155130). I am just sharing my understanding to check I got it correctly and not going to wrong direction.
Using -V resp, I can turn the option on at RPA level.

Code: Select all

DrudeWXd= ( 0.000000 , 0.000000 )  eV    # [Xd] Drude plasmon

here, first part is the Drude frequency and then second part is the damping.
Other words,
The form of the Drude model is: Y_D(w) = (w_D)^2/(w^2 + i*w*gamma_D)
where w_D and gamma_D are the real and the imaginary frequencies given by input (ev).

I have found a good paper (https://journals.aps.org/prb/abstract/1 ... .86.235147) to choose the frequency values for gold, and I am following the paper data to verify whether I am correct.
Within the Drude-Sommerfeld free electron model, they mentioned in equation 2,
epsilon_r(w) = 1- (w_p)^2/(w^2 + i*w*gamma)
where , gamma = 1/tau_D the electron relaxation rate, w_p= plasma frequency

I believe, these are the values ( w_p= plasma frequency, gamma = 1/tau_D the electron relaxation rate) ev, I need to pick up from the paper and use it in yambo.

They considered slightly different values for different gold samples as they mention in Figure 4, as follows,

FIG. 4. Dielectric function of Au (imaginary part epsilon_2) in the visible spectral region for evaporated (EV), template-stripped (TS), and single-crystal (SC) gold samples as well as the
dielectric function calculated for a Drude free electron metal with (hbar*ω_p) = 8.5 eV and gamma = 0.048 ev.

I want to run my model using the single-crystal (SC) gold sample ( 7.9 ev, 0.05 ev) and Drude free electron metal (8.5 ev, 0.048 ev) to check this time my result converges or not. Finally, my input will be:

Code: Select all

optics                           # [R] Linear Response optical properties
infver                           # [R] Input file variables verbosity
chi                              # [R][CHI] Dyson equation for Chi.
dipoles                          # [R] Oscillator strenghts (or dipoles)
FFTGvecs=  50           Ry    # [FFT] Plane-waves
Nelectro=  1216.00               # Electrons number
ElecTemp= 0.025852         eV    # Electronic Temperature
BoseTemp=-1.000000         eV    # Bosonic Temperature
OccTresh= 0.100000E-4            # Occupation treshold (metallic bands)
Chimod= "IP"                     # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% QpntsRXd
    1 |  1 |                       # [Xd] Transferred momenta
%
% BndsRnXd
    275 |  800 |                       # [Xd] Polarization function bands
%
% EnRngeXd
  0.00000 | 1.00000 |         eV    # [Xd] Energy range
%
% DmRngeXd
 0.04000 | 0.04000 |         eV    # [Xd] Damping range
%
ETStpsXd= 101                    # [Xd] Total Energy steps
DrudeWXd= ( 8.500000 , 0.048000 )  eV    # [Xd] Drude plasmon
% LongDrXd
 1.000000 | 0.000000 | 0.000000 |        # [Xd] [cc] Electric Field
%

Would you please check, professor, my thinking is going right direction or not? Thank you so much.

Best
Md J Hasan
PhD student
Mechanical Engineering
University of Maine

[muhammadhasan]

Hi Professor,

After using the DrudeWXd, I am facing the same problem as before. Here is the plot after employing Double grid (k 729 >> k 4096)

no_rigid_shift.png

However, this time I notice one interesting thing. If I somehow apply rigid shift backward on k 729 plot (~ $omega - 0.46ev), my result match with the experimental. As you can see here (Please ignore the negative values as I shift it manually):

rigid_shift.png

What I remember, some months ago I was working on temperature dependent dielectric function of Silicon (using Quantum Espresso). My result was very good agreement with experimental when I apply rigid shift (forward) on my silicon data (~ $omega+0.77 ev). I used LDA pseudopotential for silicon.

Here for the gold, I am using GGA type pseudopotential. But, I am subtracting (back shift) 0.46ev, Is that OK to you? Do you have any suggestion I can follow, professor? I am sharing the neccessary files here:

eps_Au.tar.gz

Best
Md J Hasan
PhD student
Mechanical Engineering
University of Maine

[andrea.ferretti]

Dear Md J Hasan,

here you can find a recent work where the dielectric function of gold has been computed:
https://www.nature.com/articles/s41524-019-0266-0

Despite calculations have not been done with Yambo, this could be a good reference to compare with.
As alluded to by Daniele Varsano, the metallic intraband contribution is crucial to describe the dielectric function of gold.
This is not explicitly included in Yambo at q=0, while it is at any finite q.
You have therefore 2 options:
* either you include the Drude model in the calculation at q=0
* or you compute the dielectric function for a small but finite q

in both cases you need a very dense k/q grid (this is also very crucial).

One question:
why do you need a 4x4x4 super cell ? Of course, if you can work with the unit cell (1 atom/cell) that would require larger k grids but would result in much less demanding calculations

[muhammadhasan]

Hi professor Andrea Ferretti,

Thank you so much for your reply.

I have included temperature dependent atomic positions in scf and nscf (using Quantum Espresso) calculation to include phonon effect. According to the paper https://doi.org/10.1103/PhysRevB.94.075125, we will get more accurate temperature dependent atomic positions when we use large supercell (e.g, 4x4x4, 6x6x6, 8x8x8 etc). That's why, I have run my simulation considering 4x4x4 (64 atoms) and I got very good agreement with phonon dispersion curve of Au.

That being said, I need the dielectric function of Au to all of the q points (not only q=0). Because, I need to solve the heat transmission function equation as follows:

eq.PNG

Here, the more q points I use, I will get more converge result. q =wave vector, d= vacuum separation distance of two Au, A=Area, epsilon^-1 is inverse dielectric function (each q point).

I prepared two SAVE files considering k-grid 729 points and 3375 points respectively. Another problem is, I cannot use symmetry option in nscf calculation because when I apply temperature, it breaks the symmetry (symmetry doesn't obey), so I also cannot get any help from symmetry.

My target was, I would consider the 3345 points (that means 3345 q dependent epsilon files) to solve the transmission equation. I got the result for q=0 after 2 Days of running (3 nodes, 10 cpus at each node, 50 gb memory at each cpu). I totally feeling bad to realize that if each q point takes 2 Days, then 3375 q points will take how many days!!!. What's happened, I am running the calculation (3345) for 2nd q points and I got into problem out of memory at same computer setup(but at q1=0 gives result)). I increase the memory to 70 gb at each cpu and it's also showing out of memory. Now I set it to 100 gb per cpu, I don't know when I will get the results.

Considering that, would you suggest me what should I do? Are 729 k points not enough for large 4x4x4 supercell. Then at least I will go for it and ignore the 3375 q points. I am totally hopeless.

Thanks

Best
Md J Hasan
Mechanical Engineering, PhD student
University of Maine

[Daniele Varsano]

Dear Hasan,

You can optimize the parallelization strategy to distribute memory among processors. The best option is to parallelize over bands.
You need to add in your input the following variables:

Code: Select all

DIP_CPU= "1 8 12"                      # [PARALLEL] CPUs for each role
DIP_ROLEs= "k c v"                    # [PARALLEL] CPUs roles (k,c,v)
X_CPU= "1 1 1 8 12"                        # [PARALLEL] CPUs for each role
X_ROLEs= "q g k c v"                      # [PARALLEL] CPUs roles (q,g,k,c,v)

This is an example considering 96 cpu. If you run with more or less cpu you need to change these assignments accordingly such the product of these two number matches the number of cpu.

You have a lot of k points, but this is an IP calculation and should not last so much time. I'm wondering if, because of the memory requirements, your nodes are swapping and hence running less than 100%.

Try to define this parallelization and see if the calculation run. Otherwise, it could be useful to recompile the code using the option --enable-memory-profile in the configure procedure. In this way, you will have the memory footprint in the log files.

About the mismatch with experiments, two possible reasons are 1) You need more k points (this can be verified looking to a trend with respect to the k grid). 2) Missing many body effect in the band structure, e.g. a GW calculation that will provide a stretching of your QP bands in respect to the DFT ones.

Best,

Daniele

[muhammadhasan]

Hi Professor Daniele,

Thank you for your suggested parallelization strategy. I am following it. I have some confusion if you give me some idea it would be great.

1) At finite q points (not q=0), Do I have to use the Drude parameters?

2) The most time consuming part, I think, to calculate the dipole matrix elements. Once I calculate it for q=0 and in the next run, when I calculate dielectric function for finite q points (keeping the same -J), Do I expect to get faster results as it doesn't need to recalculate dipole matrix again? If yes, I am confused why I couldn't get result for 2nd q point while keeping the same computer set up as for q = 0. ( it showed me out of memory).

3) I have seen some papers to use rigid shift to match with experimental, in my case, to match the results if I do rigid backward shift of frequency on my calculated result, Does it a permissible practice for DFT?

Thank you so much.


Best
Md J Hasan
PhD student
Mechanical Engineering
University of Maine

[Daniele Varsano]

Dear Md J Hasan,

1) No, at finite q all the intraband transition are taken into account
2) At finite q you do not need dipoles but <vk|e^iqr|ck-q> terms, these are not stored but calculated on the fly (see eq.2 in the yambo cheatsheet: https://wiki.yambo-code.eu/wiki/images/ ... et-5.0.pdf).
Moreover, at q=0 and finite q you have different symmetry operation, so the two calculations are not equivalent.
Anyway, if you have a problem with dipoles, you can speed up the calculations by neglecting the [Vnl.r]. This is an approximation, it needs to be verified, but it can be reasonable. In order to do that, just rename the ns.Vnl file in your SAVE directory.
3) I do not have an answer for that, you can of course apply an empirical shift, but it would be in my opinion you should have an argument for the disagreement. As mentioned earlier, possible reasons are (k grid not completely converged, calculation done on KS band structure instead of GW).

My suggestion here is to llok at the convergence trend with respect to k points, without using double grid, you can also consider grids with fewer k points and see the trend.

Best,
Daniele

[muhammadhasan]

Hi Professor,

Thank you so much for all of your suggestions.

This time, instead using 4x4x4 (64 Au atoms), I have used 3x3x3 (27 Au atoms) and this reduced the computational load significantly. In DFT run, firstly I have used 3375 k-points (energy cut 50 Ry) and 500 bands (Filled bands up to 253, Empty bands 261 to 500). Then after convergence test in yambo, I selected the number of bands 234-280. Based on the converged bands, in my next DFT nscf run of 4913 k-points, 5832 k-points, I have used 280 bands (instead of 500), which also reduced the computational load. I have also used the energy cutoff in yambo 20 Ry to reduce the memory load.

One of your suggestions was for speed up the calculations by neglecting the [Vnl.r]. I have seen [Vnl.r] is very significant in my calculation and decided to keep it. Here is the result of imag epsilon (Au).

Au_imag_eps.png

As you can see from the graph, I have found a very good agreement with experimental results (Imag epsilon) by selecting DmRngXd 0.018 ev (Without doing GW calculation and emperical shift). Also you can see, the blue solid line (3375 k-points, DmRngXd 0.1 ev) is away far from the experimental values. I have also compared the results with the suggested paper by professor Andrea Ferretti. This paper uses IP approximation (red circle), but not using yambo. However, using yambo (IP RPA), we got far better result than that. One important thing I have noticed that, if I change the DmRngXd to higher value (say 0.5 ev, 0.6 ev, etc.), the curve is very smooth (but going away from experiment result). I have also shown it (blue solid curve) for comparison.
So my question is:
using DmRngXd 0.018 as I have found good matched, however, the curve is not smooth, Is that OK? or do we need to change the value until the curve is reasonably smooth?

Next, I want to show you the result of real epsilon:

Au_real_eps.png

Even though, the imag part of epsilon is matched, real part is not matching. Also, the curves at higher k-points are nearly fall off in very lower energy range. I am predicting that if I can change the DmRngXd, we can also get a good match, however, in that case two different DmRngXd are unfeasible.
Would you please have any suggestions regarding that?

Thank you so much again.

Best
Md J Hasan
PhD student
Mechanical Engineering
University of Maine