Geometry optimization for n-layers in VASP

What is the best way to optimize monolayer geometry in VASP? Should the same method be used for an n-layer (n=2,3,4,5)?

What is the best way to optimize monolayer geometry in VASP?

For the geometric optimization of the monolayer in VASP, you should use the following key tags:

 ISIF=4             % firstly using 4 then 2
IBRION=2
NSW=300
EDIFFG=-0.005


You can search the explanation for each tag in VASPWIKI. For completeness, I give an INCAR template for geometric optimization in VASP.

System=Monolayer
ISTART=0       !startjob: 0-new 1-cont 2-samecut
ICHARG=2       !charge: 1-file 2-atom 10-const
ENCUT=500      !energy cutoff in eV
EDIFF=1E-6     !stopping-criterion for electronic upd.
NELM=300       !nr. of electronic steps
ISMEAR=0       !part. occupancies: -5 Blochl -4-tet -1-fermi 0-gaus 0 MP
SIGMA=0.05     !broadening in eV -4-tet -1-fermi 0-gaus
IALGO=38       !algorithm: use only 8 (CG) or 48 (RMM-DIIS), default CG algorithm (IALGO=38)

Dynamic:
ISIF=4         !2:relax ions only; 3:also relax volume and cell shape; 4:relax ions+cellshape, volume=fixed
IBRION=2       !ionic relaxation: 0-MD 1-quasi-New 2-CG
NSW=300        !number of steps for ionic upd
EDIFFG=-0.005  !stopping-criterion for ionic upd

Output:
LCHARG=.FALSE. !don't create CHGCAR
LWAVE=.FALSE.  !don't create WAVECAR


I assume that you can generate POTCAR and KPOINTS file (see another answer) for your calculation. Note the lattice constant in POSCAR of your monolayer should take the experimental lattice constant if exists. Or you can take the other answer's strategy. After all these input files are prepared, you can perform your calculation.

Should the same method be used for an n-layer (n=2,3,4,5)?

Almost you can use the previous tags. However, you should add one more tag to consider van der Waals interaction between layers, which is important to the simulation of n-layers 2D materials. There are three main strategies to consider van der Waals interaction.

#Strategy A:
IVDW = 11

#Strategy B:
LUSE_VDW = .TRUE.
GGA = MK
PARAM1 = 0.1234
PARAM2 = 1.0000
LUSE_VDW = .TRUE.
AGGAC = 0.0000

#Strategy C:
LUSE_VDW = .TRUE.
GGA = BO
PARAM1 = 0.1833333333
PARAM2 = 0.2200000000
LUSE_VDW = .TRUE.
AGGAC = 0.0000


For more strong interlayer interaction, you should use the scan+rvv10 method (VASP 5.4.4 or more recent version):

METAGGA = SCAN
LASPH = T
LUSE_VDW = T
BPARAM = 15.7


In addition, if you POSCAR contains lots of atoms with n-layer structure, larger than 10, you should add:

 LREAL=auto.


May it helps.

• Could you give more information? Commented Aug 17, 2020 at 13:59
• Do you need what? INCAR/KPOINTS/POTCAR/POSCAR are basic inputs,you should generate them easily.
– Jack
Commented Aug 17, 2020 at 14:04
• @Jack I think Cody is maybe looking for some explanation of the keywords and why these might be good options to use. As is the answer is a a little short.
– Tyberius
Commented Aug 17, 2020 at 14:56
• This also ignores some possible key details that depend on the system. For some multilayer systems, neglecting to consider dispersion corrections will give you very wrong answers. Commented Aug 18, 2020 at 14:46
– Jack
Commented Aug 19, 2020 at 8:27

I highly recommend reading: Efficient creation and convergence of surface slabs

The following answer will assume a reasonable level of VASP knowledge (where keywords can be looked up at the VASP wiki).

The best way to optimize a monolayer or surface in VASP follows:

• First, optimize your bulk structure. This will give you a reasonable estimation.
• From the optimized bulk structure, form your monolayer or surface. There are many codes that can do this for you. I recommend pymatgen.
• Introduce a vacuum layer of about 15 A, to limit interactions between periodic images.
• You will now want to run the same INCAR file you used to optimize your bulk structure with the difference: ISIF = 2.
• You should also change your KPOINT file to k k 1; where k is equal to the number of points used to optimize your bulk structure and 1 is set in the direction of the vacuum.

The ionic relaxation of your INCAR file should take the form:

IBRION = 2
NSW = 200
EDIFFG = -1E-02
ISIF = 2


Your KPOINT file should look like:

Automatic mesh
0
Gamma
k   k   1
0.  0.  0.


Note: This is a gamma centred mesh, which is often advantageous. If you are doing any kind of surface calculations, I also recommend the use of the revised for solids PBE (PBEsol) functional. This has been proven to give better results than PBE and other GGA functionals.

If you wish to deal with magnetism, then this is a lot harder with a few pitfalls. To understand these pitfalls, I would recommend asking this as a separate question. However, the paper 'Noncollinear Relativistic DFT + U Calculations of Actinide Dioxide Surfaces' offers a detailed explanation.

There is a patch of vasp used to fix any axis. For a monolayer materials, to fix c-axis is a good choice.

https://github.com/Chengcheng-Xiao/VASP_OPT_AXIS

When you want to optimize n-layer materials (n=2,3,4,5), Van der Waals (vdW) correction needs to be added in INCAR to condider the interaction between two layers. Usually I use IVDW = 11 to describe van der Waals interactions in my calculation.

• Note that the reason to use IVDW is that IVDW can work when spin-orbit coupling is considered. You can also choose other vdW correction function. Commented Apr 5, 2021 at 14:01
• Useful comment. I am using IVDW=11 for MXenes (layered material) for elastic property calculations without LSORBIT=.TRUE because there are no such heavy elements. What can go wrong there? Commented May 7 at 13:09