The bug of surface construction

[cql-cxl]

The primitive cell is shown in the following:
Primitive Cell
1.000000
2.06900038067880 1.19453792673834 9.54666700000005
-2.06900038067880 1.19453792673834 9.54666700000005
0.00000000000000 -2.38907585347668 9.54666700000005
Bi Se
2 3
DIRECT
0.8989999999999998 0.8989999999999998 0.8989999999999999 Bi1
0.1010000000000002 0.1010000000000001 0.1010000000000001 Bi2
0.5000000000000000 0.5000000000000000 0.5000000000000000 Se1
0.7059999999999997 0.7059999999999997 0.7059999999999997 Se2
0.2940000000000003 0.2940000000000002 0.2940000000000002 Se3

However, when I use the tag MILLER_INDICES 1 1 1, we obtan the wrong slab structure by running wannier_tools-2.5.1:
POSCAR for slab defined by SURFACE card and Nslab and Vacuum_thickness_in_Angstrom in wt.in by WannierTools
1.0
-2.069000 -3.583614 0.000000
2.069000 -3.583614 0.000000
4.669699 2.696052 21.546667
Bi Se
2 3
Cartesian
1.300349 0.750757 31.747360 Bi
1.300349 0.750757 8.892640 Bi
1.300349 0.750757 20.320000 Se
1.300349 0.750757 26.219840 Se
1.300349 0.750757 14.420160 Se

Is there any problem for the surface construction?

[ElAlderson]

Hello @cql-cxl,

There is indeed a bug in the surface construction. In my experience (with WannierTools v2.5.1), this occurs if the surface to be studied is other than (001) surface of the supplied unit cell (or (0001) surface for hexagonal unit cells). Although the generated POSCAR-rotated is correct, the POSCAR-slab will be wrong (one may also check these in VESTA).

From what I have understood, POSCAR-slab is generated by repeating POSCAR-rotated a finite number of times (NSLAB) and then adding vacuum on both sides. But in v2.5.1, the top few and bottom few layers are wrong. The only way I could find to bypass this bug is taking the POSCAR-rotated as the starting point for the DFT calculations and Wannierization. Then the SURFACE card would simply be identity matrix (or MILLER_INDICES 0 0 1) and no unitary transformation is done on the supplied unit cell in WT.

One other thing to keep in mind is the termination. One has to shift the POSCAR (or POSCAR-rotated) before performing DFT such that in WT, while creating the slab, the correct termination is obtained. For Bi2Se3, the termination should be at the Te-layer of the full quintuple layer.

So, first generate POSCAR-rotated from your supplied standard rhombohedral POSCAR of Bi2Se3 by setting the following SURFACE card:

SURFACE
1  -1   0
0   1  -1
0   0   1

MILLER_INDICES 1 1 1 will have the same effect. But it is better to avoid it since it is a legacy feature. The generated POSCAR-rotated may not have the correct termination, so give a shift to the z-coordinate to obtain the correct termination (this can be done in VESTA). Then, perform DFT calculations and Wannierization using the shifted POSCAR-rotated. Finally, in WannierTools, set the following SURFACE card (MILLER_INDICES 0 0 1 will have the same effect):

SURFACE
1   0   0
0   1   0
0   0   1

This should generate the correct POSCAR-slab and correspondingly the correct surface states of Bi2Se3.