My Community

Hi,

I want to calculate the band structure of a 2D material with an applied electric field in the z-direction. First, I aim to verify if I can reproduce the original band structure with zero electric field applied. To do this, I used the following INCAR file for scf calculation:

ISTART = 0
ICHARG = 2
LORBIT = 11

ISMEAR = 1
SIGMA = 0.03
ENCUT = 400

LSORBIT = .TRUE.
LMAXMIX = 4
ISPIN = 1
SAXIS = 0 0 1
MAGMOM = 9*0

ALGO = Normal
AMIX = 0.1
BMIX = 0.00001
EDIFF = 1E-7

EFIELD = 0.8
IDIPOL = 3
LDIPOL = .TRUE.
DIPOL = 0.00 0.00 0.50

However, this did not reproduce my original results. Could you please help me resolve this issue?

Finally, I would like to calculate the band structure with an applied electric field of 0.8 eV perpendicular to the monolayer.

Could you elaborate please? Can you provide the two sets of input files that you used to produce the two different outputs? If the calculation is very expensive to run, please try to reproduce the effect in a simple cell.

The INCAR file that you provided has already an electric field so I would not expect it to yield the same results as one without.

Hello,

Thank you for your quick response, and I apologize for the confusion. I have set 'EFIELD = 0.0' in the above INCAR file. Here's a detailed explanation of my process:

1. I started with the following INCAR file for the self-consistent field (SCF) calculation of a unit cell of 2D NbSe₂. This gave the correct band structure after non-SCF calculations:

ISTART = 0
ICHARG = 2
LORBIT = 11

ISMEAR = 1
SIGMA = 0.03
ENCUT = 400

LSORBIT = .TRUE.
LMAXMIX = 4
ISPIN = 1
SAXIS = 0 0 1
MAGMOM = 9*0

ALGO = Normal
AMIX = 0.1
BMIX = 0.00001
EDIFF = 1E-7

For the band structure calculation, I defined the desired K-path (M-K-Γ-K’-M’) and modified only the following two flags, leaving everything else unchanged from the initial INCAR file:

ISTART = 1
ICHARG = 11

2. Next, I aimed to calculate the band structure of 2D NbSe₂ with an applied electric field perpendicular to the plane. Before applying a finite electric field, I first attempted to reproduce the original band structure with zero electric field, using the same K-path and following INCAR file for SCF calculation:

ISTART = 0
ICHARG = 2
LORBIT = 11

ISMEAR = 1
SIGMA = 0.03
ENCUT = 400

LSORBIT = .TRUE.
LMAXMIX = 4
ISPIN = 1
SAXIS = 0 0 1
MAGMOM = 9*0

ALGO = Normal
AMIX = 0.1
BMIX = 0.00001
EDIFF = 1E-7

EFIELD = 0.0
IDIPOL = 3
LDIPOL = .TRUE.
DIPOL = 0.00 0.00 0.50

After the SCF calculation, I carried out the band structure calculation, modifying only the same two flags as in case 1, and keeping the rest of the parameters unchanged with the same K-path:

ISTART = 1
ICHARG = 11

However, the results were not same as case 1, and the output showed random curves that were clearly incorrect. Could you advise on what might be wrong with this approach?

3. Finally, I intend to calculate the band structure with an electric field of 0.8 eV/Å applied perpendicular to the monolayer. Could you guide me through the correct steps for this calculation?

Thank you in advance for your help.

If you compute a surface, you should (almost) always apply electrostatic corrections (IDIPOL/LDIPOL) even if you do not want to apply an electric field. Without these VASP will remove any dipole from your cell to achieve a continuous potential, i.e., any system that would have a dipole will be incorrectly described. By activating the corrections, VASP is free to include a discontinuity in the potential so that a dipole can form.

So for your calculation, I would predict that the first calculation has a different electronic configuration because VASP rearranged the electrons such that they cancel out any dipole the ions induce. In the second calculation, you allow for a nonzero dipole moment. You can run

grep dipolmoment OUTCAR

to check that. If this output is nonzero than only the second approach is correct.

Could you please outline the steps for applying an external electric field to a 2D material in the z-direction? I have found the following flags—are there any additional flags that I need to specify to incorporate this effect:

EFIELD = 0.8
IDIPOL = 3
LDIPOL = .TRUE.
DIPOL = 0.00 0.00 0.50

Also, do I need to perform relaxation calculations with these flags before the SCF calculations, or is that unnecessary?

Yes these flags are correct assuming your surface normal is in z direction. DIPOL is not required, VASP will determine automatically where to put the divergence but you can avoid possible issues by setting it manually. Just make sure to have the correct position.

Additionally you may want to look into WRT_POTENTIAL to visualize the potential or LVACPOTAV to get a better plot of the averaged potential.

Regarding relaxation it depends what you want to simulate. If it is a short field pulse like light where the atoms would not have time to adjust not relaxing the structure would be appropriate. If the field is permanent the atoms would relax in reality then computing without relaxation is an approximation that may or may not be justified.

martin.schlipf wrote: ↑Tue Sep 24, 2024 9:06 pm

If you compute a surface, you should (almost) always apply electrostatic corrections (IDIPOL/LDIPOL) even if you do not want to apply an electric field. Without these VASP will remove any dipole from your cell to achieve a continuous potential, i.e., any system that would have a dipole will be incorrectly described. By activating the corrections, VASP is free to include a discontinuity in the potential so that a dipole can form.

So for your calculation, I would predict that the first calculation has a different electronic configuration because VASP rearranged the electrons such that they cancel out any dipole the ions induce. In the second calculation, you allow for a nonzero dipole moment. You can run

Code: Select all

grep dipolmoment OUTCAR

to check that. If this output is nonzero than only the second approach is correct.

You mentioned that when performing calculations for surfaces or 2D slabs, I should always use the LDIPOL tag. I applied this tag in the SCF calculations for the 2D MoSe₂ structure, but the calculation did not converge. However, when I removed the IDIPOL and LDIPOL tags, the calculations converged successfully.

I have attached my INCAR and POSCAR files for your reference. Could you please review them and let me know which flags I should modify to achieve convergence? I have already tried adjusting the mixing parameter, but that did not resolve the issue.

Thank you for your guidance.

=================================
INCAR

ISTART = 0
ICHARG = 2
LORBIT = 11

ISMEAR = 0
SIGMA = 0.03
ENCUT = 400
LMAXMIX = 4

ALGO = NORMAL
AMIX = 0.1
BMIX = 0.00001
EDIFF = 1E-7

IDIPOL = 3
LDIPOL = .TRUE.

=================================

POSCAR

0
1.0
3.3157939911 0.0000000000 0.0000000000
-1.6578969955 2.8715618300 0.0000000000
0.0000000000 0.0000000000 23.3586500000
Mo Se
1 2
Direct
0.333333343 0.666666687 0.500014544
0.000000000 0.000000000 0.571561337
0.000000000 0.000000000 0.428467691

=================================

Hi,

I just noticed your post and I am currently answering a similar problem for the same structure: https://vasp.at/forum/viewtopic.php?t=19952 . I found that VASP does not place the dipole appropriately in this structure (I think at least). Setting manually:

DIPOL = 0.0 0.0 0.5

the dipole between two of the monolayers made the calculations converge in that example. Maybe give that a try.

Best regards,
Alex

Hi,

Thank you for your quick response. Could you please explain why VASP does not correctly assume the center (DIPOL) for this structure? More generally, what types of structures require a manual selection of the DIPOL center?
Additionally, why did you choose this center (0 0 1) instead of something else, such as (0 0 0.5)?

Thank you.

Well I am by no means an expert, so if anyone sees this and is more experienced please do correct me in case it is wrong, but I read a bit about it here: https://christoph-wolf.at/2018/05/02/di ... -espresso/ . This is for QE but this of course holds the same way for any other DFT code. It is more about the fact that the blog post explains where the dipole should be placed. Since your mono-layer is in the middle of the cell I thought it is best to place the dipole in the vacuum, as explained in the blog post. Hence, my choice of 0.0 0.0 0.5 . But of course this should be tested. I just observed that the calculation seems to converge then w/o problems, and I understood from the VASP manual that it will probably not place the dipole for the correction there. I hope that helps.

Best,
Alex H.

If you run into convergence problems, you can converge the structure without dipole corrections first and then afterwards activate them. Keep in mind that these corrections are small in particular if you do not have a large field or if you have a large amount of vacuum.

Alex already mentioned that the position of the dipole correction should be placed in the vacuum. You should do this because the dipole correction is a discontinuity in the electrostatic potential so it need to be somewhere where there is no electronic density. Otherwise the integral over potential * density does not converge.

Finally, if you are using VASP.6.5, you can also try the new kernel truncation method. This should allow you to use a smaller cell while getting the same results.

I am calculating the normal electric field-induced band structure (BS) of monolayer MoS2. When I explicitly specify the dipole location as mentioned in the previous discussion:

DIPOL = 0.0 0.0 0.5

it produces a drastically different band structure compared to the one calculated without the DIPOL flag. How should I determine which band structure is correct? When should I explicitly specify the DIPOL parameter, and when should I let VASP handle it automatically?

Here is my SCF file (with the DIPOL flag):
ISTART = 0
ICHARG = 2
LORBIT = 11

ISMEAR = 0
SIGMA = 0.03
ENCUT = 400

LSORBIT = .TRUE.
LMAXMIX = 4
ISPIN = 1
SAXIS = 0 0 1
MAGMOM = 9*0

ALGO = NORMAL
AMIX = 0.1
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
EDIFF = 1E-7

EFIELD = 0.2
IDIPOL = 3
LDIPOL = .TRUE.
DIPOL = 0.0 0.0 0.5

If you want to activate dipole corrections, you obtain a discontinuity in the electric potential. The discontinuity is a plane normal to the IDIPOL direction and goes through the point DIPOL. It is important that this plane does not cut through the crystal otherwise the results will be incorrect. VASP will usually do a reasonable job determining this cut but sometimes it may fail to do so. Therefore we always recommend that experienced users take this choice in their own hand.

This is a general approach in VASP: Make a reasonable default choice that will work most of the time so that if you forget to set one tag you can usually still use the results. Give the user the option to overwrite this choice to address the systems where the default does not work.

So to your question, which one is more appropriate: Whichever setup leads to the planes not touching the 2d layer.

Edit: It is not correct that the cut goes through DIPOL, see below.

Thank you for your response. I still have some doubts and would be truly grateful if you could help clarify them.

My structure is centered around z = 0.25, and based on your explanation, all DIPOL centers should be valid as long as they do not cross the 2D structure. I have tried several DIPOL values (z = 0.5, 0.75, 1.0) while ensuring that the discontinuity plane does not cut through the structure. However, these choices result in significantly different band structures.

Could you kindly take a look at my structure and INCAR file to help me understand where I might be making an error? Additionally, I also allowed VASP to determine the DIPOL center automatically, which led to yet another different result. Where can I check what DIPOL value VASP assumed for the calculation?

I have attached the reference structure, INCAR file, and band structure results for your review:

Initial structure (Before relaxing)
Mo S2
1.0
3.1689999104 0.0000000000 0.0000000000
-1.5844999552 2.7444344270 0.0000000000
0.0000000000 0.0000000000 12.3240003586
Mo S
1 2
Direct
0.333333343 0.666666687 0.250000000
0.666666627 0.333333313 0.376999974
0.666666627 0.333333313 0.123000026

INCAR FILE
ISTART = 0
ICHARG = 2
LORBIT = 11

ISMEAR = 0
SIGMA = 0.03
ENCUT = 400

LSORBIT = .TRUE.
LMAXMIX = 4
ISPIN = 1
SAXIS = 0 0 1
MAGMOM = 9*0

ALGO = NORMAL
AMIX = 0.1
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
EDIFF = 1E-7

EFIELD = 0.2
IDIPOL = 3
LDIPOL = .TRUE.
DIPOL = 0.0 0.0 (0.5/0.75/1.0)

Sorry, I misremembered how this feature worked. DIPOL should not specify the plane of the cut but the center of your slab.

Here is the algorithm that VASP uses to create the cut:

1. Get the number of points N in the grid along the direction IDIPOL. For IDIPOL = 3, N is the value NGZF in the OUTCAR file.

2. Computing the center of the cell by multiplying N with the selected value in DIPOL. Only the IDIPOL's component of DIPOL is actually (exception IDIPOL = 4).

3. Add N/2 to the center of the cell and allow for a potential discontinuity there.

If you do not set DIPOL, VASP will determine the minimum of the plane-averaged density and put the cut there. It will then compute DIPOL by adding N/2.

VASP reports the position of the cut in a somewhat convoluted manner. You will find the line direction 3 min pos X in the OUTCAR file. The X is the index of the center + N/2. If you set DIPOL in the INCAR file that should be a constant, otherwise it will oscillate a bit depending where VASP determines the minimum of the density to be. So to get the equivalent DIPOL setting, you need to take this value, divide it by N and then subtract or add 0.5 depending on how you put your lattice.