RIM and cutoff

[alsaidi]

Hello-

I am trying to understand more the convergence of the GW bandgap or eigenstates with slab thickness using the Si surface example file provided in the tutorail.
I only modified the kgird to 8x8x1 for both SCF and for generating the KSS file.

Using the RIM method using:

RandQpts=300000 # [RIM] Number of random q-points in the BZ
RandGvec=100

I noticed that the bandgap converges rather slowly with slab thickness.
For acell_z= 20,25,30,35,40, and 50 bohr, the bandgap is 1.46, 1.56,1.61,1.65,1.69, and 1.78 eV without the RIM method, and is 1.40, 1.52, 1.6, 1.66, 1.71 and 1.81 eV with the RIM method.

Is the expected? The RIM approach does not seem to be helping.

From the forum, it is best to use the RIM in conjunction with a Coloumb cutoff.
I tested this but I was not sure how to choose the cutoff for my slab.
How do I choose the cutoff for the box for a supercell of size (a,b,c)?
and do I have to center the system in the cell in this case.

Thank you.

Best,

Wissam Saidi
University of Pitt.

[Daniele Varsano]

Dear Wissam,
I've not well understood what you are trying to do:
are you varying the slab thickness, or the empty space surrounding the slab?

From the forum, it is best to use the RIM in conjunction with a Coloumb cutoff.

The RIM is mandatory for a box shaped cutoff, whereas is useless for other shapes, in the sense that it
has not effect.

If you are varying the empty space ( as I suppose) you should use a box with infinte xy direction,

CUTGeo= "box z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere
% CUTBox
0.000 | 0.000 | zcut | # [CUT] [au] Box sides
%

the value of the zcut depend on the size of your system and the side of the supercell, and have to be chosen in order
that all the electrons in your supercell do interact without interacting with the replica of your system, so big supercell
are also needed.
Remember that the RIM is needed before constructing this kind of cutoff potential.
You can also si how it looks in real space removing the flag "#CUTCol_test" that appears in the input
when generating with yambo -c -RL, and may be it is a good idea to look at it as the cutoff potential
for this kind of geometry has not be well tested.

Hope it helps.

Daniele

[alsaidi]

Dear Daniele,

Just to clarify, I am trying to increase the supercell size and not the slab thickness.
I understood from the Yambo paper that the RIM approach helps in convergence for reduced dimensional systems, but I have not seen this for the slab case. Is this correct?

For the Couloumb cutoff method applied to a slab: Assuming that the "periodic" cell lattice vectors are "a" and "b" and the physical extension along the non-periodic z-direction is "d". What is the best way to choose the supercell lattice constant "c" and the box cutoff along z?
I am trying to get this from your PRB but so far I am still confused .

Wissam
Univ of Pitt.

[Daniele Varsano]

Dear Wissam,

the Zcut has to be bigger than your system size (s), but not so big
to permit interaction with the system replica:
s < Zcut < d-s
you can have a look at this paper where also exact condition are wrote:
# An exact Coulomb cutoff technique for supercell calculations.
# C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio,
# Phys. Rev. B 73, 205119 (2006).

even if here the case of the box it is not treated, it relies on the same philosophy.
The only important difference is that for constructing the box potential you need
and integration over the Bz, and as your sampling it is not 3d, you also have to use
the RIM method.
The zcut in the input is the total, length of the size of the box , i.e. V(z) = 1/z for -zc/2 < z < zc/2,
you can check it visualizing it in real space (this is an FFT, so the more Gvec you include the more
the shpe of the potential look like to a box).

Cheers,

Daniele

[alsaidi]

Dear Daniele,

I am afraid I still have some questions. I am using a pretty big supercell , 50 bohrs along the z-axis. The system's extension is around 17 bohrs and there are 33 bohr vaccuum.
I believe that you recommend using in your PRB paper (2006) a box cutoff of L/2=25 bohr. Doing so, I get a bandgap of 2.5 eV. If I reduce the box cutoff to 20 bohrs I get 2.05 eV?
I would have thought that both should be okay? In both cases the potential projected on the xz and yz planes seems okay.

I did not center the system in the supercell? Is this relevant?

Thanks again.

Best, Wissam Saidi
Univ. of Pitt.


Here is the input file I used for the setup:
rim_cut # [R RIM CUT] Coulomb interaction
setup # [R INI] Initialization
RandQpts= 3000000 # [RIM] Number of random q-points in the BZ
RandGvec= 111 RL # [RIM] Coulomb interaction RS components
CUTGeo= "box z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere
% CUTBox
0.00000 | 0.00000 | 20.00000 | # [CUT] [au] Box sides
%
CUTRadius= 0.000000 # [CUT] [au] Sphere/Cylinder radius
CUTCylLen= 0.000000 # [CUT] [au] Cylinder length

and for the GW part:

rim_cut # [R RIM CUT] Coulomb interaction
HF_and_locXC # [R XX] Hartree-Fock Self-energy and Vxc
ppa # [R Xp] Plasmon Pole Approximation
gw0 # [R GW] GoWo Quasiparticle energy levels
em1d # [R Xd] Dynamical Inverse Dielectric Matrix
RandQpts= 3000000 # [RIM] Number of random q-points in the BZ
RandGvec= 111 RL # [RIM] Coulomb interaction RS components
CUTGeo= "box z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere
% CUTBox
0.00000 | 0.00000 | 20.00000 | # [CUT] [au] Box sides
%
CUTRadius= 0.000000 # [CUT] [au] Sphere/Cylinder radius
CUTCylLen= 0.000000 # [CUT] [au] Cylinder length
EXXRLvcs= 3089 RL # [XX] Exchange RL components
% QpntsRXp
1 | 25 | # [Xp] Transferred momenta
%
% BndsRnXp
1 | 100 | # [Xp] Polarization function bands
%
NGsBlkXp= 100 RL # [Xp] Response block size
% LongDrXp
1.000000 | 0.000000 | 0.000000 % GbndRnge
1 | 100 | # [GW] G[W] bands range
%
dScStep= 0.10000 eV # [GW] Energy step to evalute Z factors
% GDmRnge
0.10000 | 0.10000 | eV # [GW] G_gw damping range
%
%QPkrange # [GW] QP generalized Kpoint/Band indices
1| 1| 1|100|
%
%QPerange # [GW] QP generalized Kpoint/Energy indices
1| 25| 0.0|-1.0|
|
PPAPntXp= 27.21138 eV # [Xp] PPA imaginary energy

===================== Abinit input file


kptopt1 1
ngkpt1 8 8 1
shiftk1 0 0 0

acell 12.467903 7.1983 50.
xcart 5.63153281207509E+00 -3.59915000000000E+00 -6.51556659141272E-01
5.94978424629170E-01 3.99585919907963E-16 -4.03050523190062E-01
6.08969181251996E+00 -3.59915000000000E+00 2.15014461258978E+00
1.30822391984638E-01 3.99585919907963E-16 2.40042103807968E+00
2.09017759849385E+00 -3.59915000000000E+00 4.03001654213921E+00
-4.16312784010154E+00 3.99585919907963E-16 3.31388779850620E+00
2.08937918542687E+00 -3.59915000000000E+00 8.47747342177409E+00
-4.14872787258864E+00 3.99585919907963E-16 7.66996576372150E+00
-2.26223250633845E+00 -3.59915000000000E+00 9.38816468640326E+00
4.36177504392197E+00 3.99585919907963E-16 9.67172605366154E+00
-2.92762191101085E-02 -3.59915000000000E+00 1.60166759685954E+01
2.06568224183842E+00 3.99585919907963E-16 1.50655779346659E+01
-3.72624715164578E+00 -3.59915000000000E+00 1.35923712310877E+01
6.19304955132109E+00 3.99585919907963E-16 1.37734236641877E+01

[Daniele Varsano]

Dear Wissam,

It looks you have a huge difference of the gap for not so big change in the cutoff parameter.
In this moment I'm not sure that the value in the input file refers to zc/2 or to zc,
most probably to zc and (in this case your zc is quite small, as zc/2 have to be at least long
as your system).
Moreover this kind of geometry has not been carefully tested.
What I can suggest to you, if your calculation is not too heavy, is to scan different values
of zcut. In a good situation the gap should reach a plateau, as you can see from the paper.

In order to speed-up your calculation I suggest you to calculate the correction in correspondence
of the gap using the variable:

%QPkrange # [GW] QP generalized Kpoint/Band indices

Hope it is clear,

Daniele

[alsaidi]

Hello-

In my investigations of the surface, I have noticed some sensitivity of the GW eigenstates to the location of the system with respect to the simulation box.

Here are the results at different values of K-point in the BZ.

DFT_Bandgap: 1.09 1.14 1.35 1.66 1.83 2.01 2.08 2.06 2.01
GW_Bandgap 1.78 1.85 2.08 2.41 2.56 2.75 2.86 2.8 2.7
HOMO_GW: 0.08 0.04 -0.09 -0.28 -0.28 -0.10 -0.26 -0.28 -0.23
LUMO_GW: 1.86 1.89 1.99 2.13 2.28 2.65 2.60 2.52 2.47
HOMO_DFT: -0.22 -0.25 -0.37 -0.55 -0.63 -0.43 -0.55 -0.60 -0.58
LUMO_DFT: 0.87 0.89 0.98 1.11 1.20 1.58 1.53 1.46 1.43

and this is for the same system after shifting it to the center of the simulation box:

FT_bandgap: 1.09 1.14 1.35 1.66 1.83 2.01 2.08 2.06 2.01
GW_bandgap: 1.7 1.76 1.99 2.38 1.96 2.41 2.57 2.62 2.44
HOMO_GW: 0.21 0.17 0.05 -0.15 0.42 0.54 0.27 0.13 0.29
LUMO_GW: 1.91 1.93 2.04 2.23 2.38 2.95 2.84 2.75 2.73
HOMO_DFT: -0.22 -0.25 -0.37 -0.55 -0.63 -0.43 -0.55 -0.60 -0.58
LUMO_DFT: 0.87 0.89 0.98 1.11 1.20 1.58 1.53 1.46 1.43

Note that the DFT results are the same (as they should be), but not the case with the GW results.
Also, note that there is a "larger" shift in the eigenstates than in the bandgaps.
In both cases, I am using the same parameters. Also, I am not using a cutoff or doing a RIM.

Best, Wissam
Univ. of Pitt.

[Daniele Varsano]

Dear Wissam,
did you check the convergences of all the variables in both cases?
in particular of the k-point sampling? Displacing your system in the
simulation box, could results in a difference of the symmetries of
the system that Yambo take into account, and consequently a different
number of the total k-point sampling in the full Bz. You can check it in the
report files. And this can affect the GW results if, for instance,
the Brillouin zone it is not well sampled. I don't know if it is the case,
and may be someone can you more about surfaces.

Cheers,
Daniele

[andrea marini]

Dear Wissam, I agree with Daniele's remark. At difference with DFT Yambo has to treat a lot of non-local summations. This is quite different from DFT and symmetries are crucial. TO be more precise, there have been several cases that I remember where QP energies were strongly affected by the use of symmetries. I remember in particular the case of graphene nanoribbons were the final bands did not even have zero derivaties at the high symmerty k-points if the ribbon was not placed in a symmetric position.

The reason is, as Daniele, explained in a simple fact. If you write down the HF self-energy, for example, you can easily prove that ONLY by summing on a symmetric group of RL vectors the final result is symmetric. Please try yourself to prove this and you will be easily convinced that this explains your results.

[alsaidi]

Dear Daniele and Andrea,

Thanks for explaining this.

Back to the cutoff issue.
First, I did several calculations with acell_z ranging from 20 to 300 bohr. I am using a small planewave cutoff so these calculations are not expensive.
In any case, the bandgap is extrapolated to be 2.4 eV using a 1/Volume linear fit. The bandgap with acell_z=200 bohr is 2.1eV.

Now, using the truncated Couloumb interaction, I have the following results
(for each unit cell case, I show the value of the cutoff and the bandgap).

acell_z=30 bohr
40 2.66
50 2.94
55 2.3
60 3.82

acellz=40 bohr
40 2.38
50 2.57
55 2.44
60 2.6

acellz=50 bohr
40 2.46
45 2.43
50 2.39
55 2.39
60 2.4

acellz=100 bohr
40 2.35
50 2.51
70 2.49
100 2.52

Now, the systems extension is only ~18 bohr so I would have thought that a box of 40 bohr should be enough with the truncated Couloumb but the results are not very satisfactory. Of course the differences between the naked Couloumb interaction and the cutoff one are huge and most the size effects are eliminated with the cutoff case. Still the differences between acell=50 and acell=100 are noticeable.

I have the checked the convergence of most of these results. I am still checking them but I think they are converged.

By the way, it seems to me that the implementation of the box cutoff in the code is different from that of the PRB 2006 paper.

Thanks.

Best,