ionic relaxation cannot reach accuracy

[yujia_teng]

Dear admin,
I had problem with getting the reached-accuracy relaxed structure. I've did thousands of ionic steps but it just cannot relax to required accuracy. There's no problem with electronic step. For this kind of structure, if there're less number atoms of last 2 elements, the structure is relaxed to accuracy quickly.

The calculation of file 'input.zip' attached starts from the WAVECAR that comment out dipole correction tags. The 'input-IBRION-1.zip' is the calculation where I used ADDGRID=.True. and IBRION=1. For these two calculations, EDIFFG=-0.005, which is a little bit high.

Looking at the end of OUTCAR, the 'TOTAL-FORCE (eV/Angst)' column, the value is still larger than 0.01. So even if I set EDIFFG=-0.01, it would still not converge. The job with EDIFFG=-0.01 is still running and I don't see any sign of convergence.

So I'm wondering is there a good way to solve this issue.

[michael_wolloch]

Dear Yujia Teng,

Thanks for your question and for taking care to attach all of your files.

I think your conversion problems are mainly due to the use of dipole corrections. Your system is nearly symmetric, and in your OUTCAR files, you can see that your dipole is very small (less than 0.01 electrons x Angstrom):
e.g. first iteration of the OUTCAR file you provide from the input folder:

Code: Select all

DIPCOR: dipole corrections for dipol
 direction  3 min pos     1,
 dipolmoment           0.000000      0.000000      0.005787 electrons x Angstroem
 Tr[quadrupol]      -696.843705

 energy correction for charged system         0.000000 eV
 dipol+quadrupol energy correction           -0.000002 eV
 added-field ion interaction         -0.001074 eV  (added to PSCEN)
 

As you can see, the resulting energy corrections are vanishingly small as well.

Generally, we recommend relaxing a system (even if it has a sizable dipole moment) without dipole corrections and add them only after convergence is reached. This information is not well documented yet on the wiki, but we are currently working on preparing a summary page for surface calculation with much advice and many tips, and this will be among them.

Forces below 0.01 eV/A are very small and hard to distinguish from noise. It is very seldom useful to go below that threshold.

It is furthermore beneficial to clamp one ion using the selective dynamics POSCAR option. This helps to reduce possible oscillations in the layers of the slab without reducing any degrees of freedom. I chose the central Se ion, but in principle this is arbitrary, and you could also e.g., use the bottommost ion.

Here is the modified POSCAR:

Code: Select all

Zr10 In2 Te20 Se3
1.0
   4.1030924467899776    0.0000000000000845    0.0000000000000000
  -2.0515462233430743    3.5533822925326985    0.0000000000000000
   0.0000000000000000    0.0000000000000000   92.8828312293978513
In Se Zr Te
2 3 10 20
Selective Dynamics 
direct
   0.3333333329999988    0.6666666670000013    0.4770187862906877 T T T
   0.6666666670000012    0.3333333329999988    0.5224750845172059 T T T
   0.3333333329999988    0.6666666670000013    0.5043641580049560 F F F
   0.0000000000000000    0.0000000000000000    0.4633063782846610 T T T
   0.0000000000000000    0.0000000000000000    0.5372603260841553 T T T
   0.6666666670000012    0.3333333329999988    0.4098347377971909 T T T
   0.3333333329999988    0.6666666670000013    0.5900451425926466 T T T
   0.3333333329999988    0.6666666670000013    0.6607853240054014 T T T
   0.3333333329999988    0.6666666670000013    0.7315255054181562 T T T
   0.6666666670000012    0.3333333329999988    0.3390945563844361 T T T
   0.3333333329999988    0.6666666670000013    0.8730058682436657 T T T
   0.6666666670000012    0.3333333329999988    0.1976141935589267 T T T
   0.6666666670000012    0.3333333329999988    0.2683543749716814 T T T
   0.3333333329999988    0.6666666670000013    0.8022656868309109 T T T
   0.6666666670000012    0.3333333329999988    0.1268740121461719 T T T
   0.0000000000000000    0.0000000000000000    0.6093767508487827 T T T
   0.0000000000000000    0.0000000000000000    0.7508571136742922 T T T
   0.6666666659999977    0.3333333340000025    0.6416755078988421 T T T
   0.0000000000000000    0.0000000000000000    0.6801169322615374 T T T
   0.6666666659999977    0.3333333340000025    0.5709353264860875 T T T
   0.0000000000000000    0.0000000000000000    0.3906232491512173 T T T
   0.0000000000000000    0.0000000000000000    0.8215972950870469 T T T
   0.6666666659999977    0.3333333340000025    0.7831558707243517 T T T
   0.6666666659999977    0.3333333340000025    0.7124156893115969 T T T
   0.6666666659999977    0.3333333340000025    0.8538960521371064 T T T
   0.0000000000000000    0.0000000000000000    0.8923374764998017 T T T
   0.0000000000000000    0.0000000000000000    0.3198830677384625 T T T
   0.3333333340000025    0.6666666659999977    0.3582961570195209 T T T
   0.3333333340000025    0.6666666659999977    0.2875559756067662 T T T
   0.0000000000000000    0.0000000000000000    0.2491428863257078 T T T
   0.3333333340000025    0.6666666659999977    0.2168157941940114 T T T
   0.3333333340000025    0.6666666659999977    0.4290363384322756 T T T
   0.0000000000000000    0.0000000000000000    0.1784027049129530 T T T
   0.3333333340000025    0.6666666659999977    0.1460756127812567 T T T
   0.0000000000000000    0.0000000000000000    0.1076625235001983 T T T

I cleaned up your INCAR file a bit and removed the dipole corrections. Then the calculation converged in 81 steps, which is reasonable for a vdW-bonded slab.
Note that we recommend keeping the INCAR file concise without too many comments to avoid careless mistakes.

Code: Select all

 SYSTEM=Zr10In2Te20Se3
  

   ISPIN  =  1
   LORBIT =  11
   LMAXMIX = 4
   EDIFF  = 1.0E-06
   IVDW = 11
   NSW    = 200
   IBRION =    2
   ISIF   =    2
   EDIFFG = -0.01
   MAXMIX = 120 
   ALGO = Normal 
   ENCUT = 250
   ISMEAR = 0 
   SIGMA =  0.1
   NEDOS  = 2000
   EMIN   = -15
   EMAX  = 15
   AMIN = 0.01
   LREAL= Auto
   KPAR=10
   NCORE=4

I did reduce ENCUT by 50% compared to your original value of 500 eV. Since the largest ENMAX in your POTCAR file is 239 eV, this seemed a bit excessive for a relaxation without cell shape changes. To compare total energies with other calculations, you will have to take care that ENCUT is equivalent, of course.
I also added MAXMIX and gave it a positive value of around 3 times the electronic steps of the first ionic iteration (I restarted the relaxation from the WAVECAR file of an initial static run.)

I also did a follow-up relaxation with the dipole tags added back in, and IBRION=1. It converged after the first ionic step, and the dipole was essentially 0, as expected.

The OUTCAR file of my relaxation is attached. Please let me know if this resolves your issue, or if you have follow-up questions.
Michael

[yujia_teng]

Dear Michael,
Thanks very much for your detailed explanation and effort. I was not aware that it's recommended to relax system without dipole tag, and this is the key problem why I couldn't get convergence before. I followed your suggestion, comment out dipole tag to relax first, and it got relaxed soon (I have a few more steps, probably because it's different on different machine and different number of cores). Then I add dipole tag back, and it converged in 1-2 ionic steps in my several tests. (IBRION=1 and 2 doesn't make much difference here. But from wiki, looks like using 1 is more reasonable because it's very close to ground state.)

Force below 0.01 is indeed too small, so I switched to EDIFFG = -0.01, not -0.005 before. All other tips you provided could further improve the convergence (converge in fewer ionic steps), like clamping one ion and using MAXMIX.

I have one question about ENCUT here. For relaxation of atomic positions only with cell fixed, a ENCUT slightly larger than ENCUT is enough? I remember that in wiki it says we need to at least use ENCUT = 1.3 * ENMAX. Also, from my convergence test below (with fewer Zr and Te atoms), it converges at around 500 eV. But in a paper of this material, they also just use 250 eV ENCUT. So for the following scf and nscf calculations, I'm not sure which one is best to be used.

[michael_wolloch]

Dear Yujia Teng,

Convergence of absolute energies with respect to ENCUT (and really anything) is usually very slow. But it is also unnecessary, since absolute energies are entirely meaningless in DFT. You want to converge energy differences with respect to ENCUT!
A good example is to perform static calculations of a small bulk cell of your material at different volumes and fit an equation of state. This fit will give you an equilibrium volume (and also info about the bulk modulus), and it depends only on the energy differences between the different volumes. You will see that your equilibrium volume will converge quite fast with respect to the plane-wave cutoff.

Generally, the maximum ENMAX of your POTCAR file is a fine value for ENCUT. Of course, you should converge the parameter of interest with respect to ENCUT, but for a relaxation with a fixed cell, ENMAX is usually OK.
If the cell volume changes, the number of plane waves changes, since ENCUT is constant. This can lead to problems if ENCUT is too small, see the section on Pulay stress. For this reason, it is recommended to increase it by 30% if the cell is allowed to change.

I hope this answers your question. Let me know if it does, so I can lock the thread, or if it does not, so I can follow up with more details.
Michael

[yujia_teng]

Dear Michael,
Thanks for the explanation, it is quite clear and answers most of my questions. The problem of ENCUT is a little bit off-topic of this thread, so I think I should open another one to ask about this. And this thread can be closed.

Best,
Yujia