Constrained occupation not respected in VASP6

[oscarlb]

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:

vasp6_occup_test_out.tar.gz

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

After first iteration of ionic step 2:

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

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:

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")

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.

[henrique_miranda]

Apart from the post, you indicated above we've recently received a very similar user report about this issue:
https://www.vasp.at/forum/viewtopic.php?f=4&t=18725

My proposal is to modify the densta routine in dos.F:
https://www.vasp.at/forum/viewtopic.php ... 725#p22989

Let me know if this also solves the issue for you.
We are planning to include this fix in the next release of VASP.