Constrained occupation not respected in VASP6
Posted: Fri Dec 02, 2022 2:49 pm
Dear developers,
Testing constrained occupation calculations in versions 6.2 and 6.3, we notice that structural relaxation sometimes does not keep the occupation as specified, even changing the total number of electrons if the supercell is charged.
We set the constrained occupation calculations with ISMEAR=-2 and specify the occupation with FERWE/FERDO. When we use LDIAG=.FALSE. together with a hybrid functional (HSE06 in our case), the initially constrained occupation resets to the neutral ground state after the first ionic step. In contrast, in versions of VASP before 6 (we test 5.4.4), as well as in VASP 6.2-6.3 with LDIAG=.TRUE. and a GGA functional, the specified occupation is kept throughout the full structural relaxation as one may expect.
We include a minimal test case: where we calculate the first excited state of the silicon vacancy in 4H-SiC, which has been previously studied with VASP, with a negatively charged supercell of 32 atom size. The test case includes the ground state calculation and constrained occupation runs in VASP versions 5.4.4 and 6.3.0. One may, in particular, investigate the occupations seen in the first and second ionic iterations written in the OUTCAR.
Near-bandgap states for vasp 6.3.2:
End of Ionic step 1:
After first iteration of ionic step 2:
It is worth noting the second ionic step not only resets the excited electron occupation in the second spin component but also changes the number of electrons in the calculation, completely ignoring the FERWE/FERDO tags reflected in the occupation of the first step. Performing these tests for supercells of larger size (as necessary for total energy calculations of defects, we go up to 576 atoms) removes the partial occupation, but the occupation reset still occurs.
As for the source of the issue:
We know that this is normally not an issue in VASP version 5.4.4 (apart from problems mentioned in https://www.vasp.at/forum/viewtopic.php ... ion#p21231, which we patch in our version of 5.4.4), and a run with that version is included and performs as expected. Comparing the source of 6.3.0 and 5.4.4, we noticed that the occupation is reset in the part of the main.F commented as the regeneration of Fermi-weights following the "prediction of wavefunctions" (line ~4580 in 6.3.0). This is also in line with the argument that the issue appears when using LDIAG=.FALSE. as doing otherwise would trigger a false result in the "pre_subrot" if-statement if we have interpreted our parameters correctly. There was also a check in 5.4.4 for if wavefunctions were predicted or not (PRED%IPRE>1), which appears moved in 6.3.0 such that it no longer encapsulates this weight-regenerating code. Moving it back using the patch file "pre_patch.txt" between 5.4.4 and 6.3.0, and running again seems to behave identically to the 5.4.4 runs in practice.
Patch-file for main.F for 6.3.0:
This occupation behaviour and difference in source seems present in versions 6.2.1-6.3.2 and possibly earlier vasp6 versions we did not test.
We hope this is the source of the issue, but it would be best if anyone on this forum could confirm and/or check if this does not break anything else along the way.
Your attention to this issue would be highly appreciated.
Testing constrained occupation calculations in versions 6.2 and 6.3, we notice that structural relaxation sometimes does not keep the occupation as specified, even changing the total number of electrons if the supercell is charged.
We set the constrained occupation calculations with ISMEAR=-2 and specify the occupation with FERWE/FERDO. When we use LDIAG=.FALSE. together with a hybrid functional (HSE06 in our case), the initially constrained occupation resets to the neutral ground state after the first ionic step. In contrast, in versions of VASP before 6 (we test 5.4.4), as well as in VASP 6.2-6.3 with LDIAG=.TRUE. and a GGA functional, the specified occupation is kept throughout the full structural relaxation as one may expect.
We include a minimal test case: where we calculate the first excited state of the silicon vacancy in 4H-SiC, which has been previously studied with VASP, with a negatively charged supercell of 32 atom size. The test case includes the ground state calculation and constrained occupation runs in VASP versions 5.4.4 and 6.3.0. One may, in particular, investigate the occupations seen in the first and second ionic iterations written in the OUTCAR.
Near-bandgap states for vasp 6.3.2:
End of Ionic step 1:
Code: Select all
Spin-component 1
...
59 6.1461 1.00000
60 6.4252 1.00000
61 6.4253 1.00000
62 7.6738 1.00000
63 7.7898 1.00000
64 7.7900 1.00000
65 9.9708 0.00000
66 9.9708 0.00000
67 10.4378 0.00000
68 10.6074 0.00000
69 10.7180 0.00000
70 10.7181 0.00000
...
Spin-component 2
...
59 6.4948 1.00000
60 6.4949 1.00000
61 7.0620 0.00000
62 8.2484 1.00000
63 9.9884 0.00000
64 9.9884 0.00000
65 10.4848 0.00000
66 10.6524 0.00000
67 10.6524 0.00000
68 10.6791 0.00000
69 11.1340 0.00000
70 11.1342 0.00000
Code: Select all
Spin-component 1
...
59 5.9909 1.00000
60 5.8544 1.00000
61 5.8572 1.00000
62 7.0095 1.00000
63 8.7679 0.56788
64 8.9009 0.43212
65 10.2737 0.00000
66 10.2986 0.00000
67 10.3802 0.00000
68 11.0945 0.00000
69 10.8691 0.00000
70 10.8920 0.00000
...
Spin-component 2
...
59 5.9176 1.00000
60 5.9367 1.00000
61 6.5199 1.00000
62 8.0579 0.00000
63 10.2320 0.00000
64 10.2440 0.00000
65 10.4294 0.00000
66 10.6985 0.00000
67 10.7002 0.00000
68 11.0881 0.00000
69 11.0111 0.00000
70 10.9885 0.00000
As for the source of the issue:
We know that this is normally not an issue in VASP version 5.4.4 (apart from problems mentioned in https://www.vasp.at/forum/viewtopic.php ... ion#p21231, which we patch in our version of 5.4.4), and a run with that version is included and performs as expected. Comparing the source of 6.3.0 and 5.4.4, we noticed that the occupation is reset in the part of the main.F commented as the regeneration of Fermi-weights following the "prediction of wavefunctions" (line ~4580 in 6.3.0). This is also in line with the argument that the issue appears when using LDIAG=.FALSE. as doing otherwise would trigger a false result in the "pre_subrot" if-statement if we have interpreted our parameters correctly. There was also a check in 5.4.4 for if wavefunctions were predicted or not (PRED%IPRE>1), which appears moved in 6.3.0 such that it no longer encapsulates this weight-regenerating code. Moving it back using the patch file "pre_patch.txt" between 5.4.4 and 6.3.0, and running again seems to behave identically to the 5.4.4 runs in practice.
Patch-file for main.F for 6.3.0:
Code: Select all
@@ -4575,7 +4575,9 @@
! wavefunctions are not diagonal, so if they are written to the file
! or if we do not diagonalize before the optimization
! rotate them now
- pre_subrot: IF (.NOT.INFO%LDIAG .OR. INFO%LONESW .OR. &
+ pre_done: IF (PRED%IPRE>1) THEN
+
+ pre_subrot: IF (.NOT.INFO%LDIAG .OR. INFO%LONESW .OR. &
& INFO%LSTOP .OR. (MOD(NSTEP,10)==0 &
& .AND. PRED%IWAVPR >= 1 .AND. PRED%IWAVPR < 10) ) THEN
@@ -4621,6 +4623,8 @@
ENDIF pre_subrot
+ ENDIF pre_done
+
IF (INFO%LONESW .AND. INFO%NELMDL == 0 ) THEN
CALL START_TIMING("G")We hope this is the source of the issue, but it would be best if anyone on this forum could confirm and/or check if this does not break anything else along the way.
Your attention to this issue would be highly appreciated.