Geometry optimization for n-layers in VASP

Question

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)?

Answer 1

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
ADDGRID = 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.

Answer 2

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.

Answer 3

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.