SCF loop diverging with Spin-orbit on

[paulfons]

I am running a calculation using Vasp 5.4.4 on a coincidence lattice system with a total of 182 ions. The system consists of Sb2Te3 and a TMDC compound. There are three different TMDC compounds in the study. All three systems were structurally relaxed (without spin-orbit coupling) without problem. Similarly, I was able to calculate the band structure for each system without trouble (again without spin-orbit coupling). For two of the systems, the SOC calculations converge without problem, but for the third (MoS2), the system fails to converge. There are a total of 182 atoms in the system consisting of 2 ML of MoS2 and 7 QL of Sb2Te3. The computational cell has a 1:10 aspect ratio with the z-axis being significantly longer.

I then took the relaxed structures and attempted to include the effects of spin-orbit coupling. For two of the systems, the system converged fine while the third system did not. I have tried many variations of input and k points, but the system always has diverged in energy. This occurred even when ICHARG=12 (fix charges) was set. There are 1120 electrons in the system and enlarged the number of bands to 1500 to avoid numerical troubles. I have no idea what to try next hence my inquiry here.


N E dE d eps ncg rms rms(c)
DAV: 1 0.649191374646E+04 0.64919E+04 -0.38243E+05 50816 0.106E+03
DAV: 2 -0.100053666920E+03 -0.65920E+04 -0.62718E+04 57696 0.209E+02
DAV: 3 -0.800304317535E+03 -0.70025E+03 -0.67965E+03 55920 0.740E+01
DAV: 4 -0.829168447772E+03 -0.28864E+02 -0.28265E+02 59760 0.146E+01
DAV: 5 -0.830553120111E+03 -0.13847E+01 -0.13781E+01 72432 0.314E+00 0.171E+01
DAV: 6 -0.830390999775E+03 0.16212E+00 -0.61298E-01 65032 0.105E+00 0.103E+01
DAV: 7 -0.830371226976E+03 0.19773E-01 -0.17680E-01 57976 0.441E-01 0.351E+00
DAV: 8 -0.830367238150E+03 0.39888E-02 -0.56942E-02 58064 0.335E-01 0.114E+00
DAV: 9 -0.298739524471E+04 -0.21570E+04 -0.23792E+02 57368 0.325E-01 0.137E+01
DAV: 10 -0.133848642162E+05 -0.10397E+05 -0.28436E+04 55488 0.279E+00 0.530E+01
DAV: 11 -0.190743253401E+06 -0.17736E+06 -0.46255E+05 56256 0.343E+01 0.985E+01
DAV: 12 -0.903883872139E+06 -0.71314E+06 -0.60175E+06 54964 0.431E+01 0.185E+02
DAV: 13 -0.242854531447E+06 0.66103E+06 -0.11228E+06 57220 0.102E+02 0.274E+02
DAV: 14 -0.517376273239E+06 -0.27452E+06 -0.95459E+03 55588 0.295E+02 0.200E+02
DAV: 15 -0.108644484194E+09 -0.10813E+09 -0.18256E+08 52828 0.457E+02 0.171E+02
DAV: 16 -0.412246177950E+08 0.67420E+08 -0.51783E+07 54864 0.589E+02 0.154E+02
DAV: 17 -0.565707635239E+09 -0.52448E+09 -0.24748E+05 58112 0.605E+02 0.147E+02
DAV: 18 -0.433209920258E+08 0.52239E+09 -0.59808E+05 59696 0.710E+02 0.141E+02


A typical input file is attached below. Note that I have also tried ISMEAR = -1 (Gaussian smearing) and -5 (tetrahedron method with Blöchl corrections), increasing the kpoint density, both IBRION = 1 (quasi-Newton) and 2 (conjugate gradient) and more. I am using a 4x4x1 Gamma centered mesh (the system is hexagonal) although I have increased the mesh density further only to encounter a similar divergence. I should also add that after the tutorial slide on convergence problems, I also set ICHARGE=12 (to fix the charge) only to find the SCF loop diverges like above. Note that the same system with Te and Se converges fine). Any suggestions as to what to try next would be gratefully received.



ALGO = Normal
EDIFF = 1E-5
ENCUT = 420
IBRION = 2
ICHARG = 1
ISIF = 3
ISMEAR = -5
ISPIN = 2
ISYM = -1
LAECHG = True
LCHARG = True
LORBIT = 11
LREAL = False
LSORBIT = T
LVHAR = True
LWAVE = False
NELM = 100
NSW = 0
PREC = Accurate
SAXIS = 0 0 1
SIGMA = 0.05

NCORE = 9
AMIN = 0.01
NBANDS = 1500

[admin]

The solution to your problem is found in steps.
In the following description the vasp_ncl version was used in every step, even for collinear calculations, in order to have compatible WAVECARs.

If convergence problems appear for SOC/NCL calculations,
we usually try to obtain a solution without SOC/NCL at the Gamma point only and use this WAVECAR for the actual SOC/NCL calculation.

So following steps were done to obtain the solution for the full k-point grid with SOC turned on in the non-collinear spin orientation.

1. Gamma-only calculation for LSORBIT = .FALSE. ; LNONCOLLINEAR = .FALSE. ; keep WAVECAR
2. use WAVECAR of previous step with LSORBIT = .FALSE. ; LNONCOLLINEAR = .TRUE. ; keep WAVECAR
3. use WAVECAR of previous step with LSORBIT = .TRUE. ; LNONCOLLINEAR = .TRUE. ; keep WAVECAR
4. Calculation with full k-point grid (gamma-point must be included in grid) based on WAVECAR of previous step

Very often, step 2 can be skipped.
This does not apply to your system, but for magnetic systems in general the MAGMOM tag should be set in all steps to study the desired spin-orientation.
This means, the user should set MAGMOM in step 1 (and 2) with LORBIT = 11, which results in an table of the form in the OUTCAR

Code: Select all

 magnetization (z)
 
# of ion       s       p       d       tot
------------------------------------------
    1       -0.000   0.000  -0.000  -0.000
    2       -0.000   0.000  -0.000  -3.753
    3       -0.000   0.000  -0.000  -0.406
    4       -0.000   0.000  -0.000  -0.000

The values in the last column can be used for MAGMOM in the INCAR for the non-collinear setup in step 3 and 4.
Note, however, that here 3 values for each ion should be supplied, as written on our wiki site:
https://cms.mpi.univie.ac.at/wiki/index ... NCOLLINEAR