Does the plot of the workfunction correspond to the INCAR file that you used? If you use dipole corrections, you should see a discontinuity in the vacuum region. We have a test in the testsuite that should specifically test this case (surface_AlC_VACPOTAV). Perhaps you could run that and see if you see a workfunction that you expect.
Generally, we recommend to start with a simple system and then gradually move towards more complex ones. In your case, you have a 2x2 cell in the plane and 6 layers of material. You could try using less layers and a 1x1 cell first. Furthermore, MAPI should not be magnetic; if you use
you will see that the magnetization disappears over the iterations. Calculating without ISPIN = 2 will be about twice faster. I am also not sure whether the mixing tags you set will really help you.
So here is a step-by-step guide:
Run the example from the testsuite, check the workfunction e.g. with
Code: Select all
import py4vasp
py4vasp.calculation.workfunction.plot()
Switch to a minimal example of your cell, e.g., 1x1 cell with 3 layers and a bit smaller vacuum. A 4x4x1 k-point mesh should probably be sufficient to get a qualitative picture. Use the same INCAR file as the one in the testsuite (except for the DIPOL position).
Switch to the INCAR setting that you want to use for your production run, i.e., add vdW and other tags.
Extend the cell to the real cell that you want to use, reduce the k-point mesh to 2x2x1.
Converge the k-point mesh until your property of interest does not change anymore.
If you run into problems at any of these steps, please get back to me.
Finally, please remember to check your results with SOC activated. SOC has a strong impact on the electronic structure in MAPI so you need to make sure that it does not influence the property that you are interested in.