Consultation on the zero filed splitting(ZFS) calculation and LDMATRIX tag

[xin_zhao]

Dear VASP developers,

Hello, I frequently perform calculations to obtain zero-field splitting (ZFS) matrices, previously using a custom program developed by our research group that processes VASP outputs to compute the ZFS matrices. Recently, I noticed on a GitHub discussion page that VASP seems to have an internal feature for calculating ZFS matrices(https://github.com/alejandrogallo/dmatrix/issues/1), marked by the tag LDMATRIX. I tried searching for this tag within the VASP source files and eventually found it in src/dmatrix.F. Although I am not familiar with Fortran, I have a general understanding of how ZFS calculations are implemented, and this file appears to compute the ZFS matrices while also accounting for PAW corrections.

However, it is quite strange that there is no mention of LDMATRIX or related information on ZFS matrices in the VASP wiki. I have also been unable to find any further documentation on this tag beyond the GitHub page mentioned earlier. Despite this, I decided to test LDMATRIX to see if it could produce accurate ZFS matrices.

For testing, I selected systems that have been widely studied: NV centers in diamond and the hh and kk types of NV defects in 4H-SiC. Unfortunately, the test results have been disappointing. There is a notable deviation from both experimental values and previous calculations. For example, for NV centers in diamond, the experimental D value is around 2.8–2.9 GHz, with previous literature (also using PAW methods) confirming values within this range. However, VASP outputs around 3.6 GHz, significantly higher than expected. Similarly, for NV centers in 4H-SiC, both experiments and other calculations suggest a D value around 1.4 GHz, while VASP calculates approximately 1.1 GHz, which is lower than expected.

Additionally, all results mentioned above were obtained with symmetry enabled (ISYM=2). In this setting, VASP outputs a ZFS matrix that respects the system's symmetry. However, with the same structure but ISYM=0, the resulting ZFS matrix differs from the one obtained with ISYM=2 and lacks the expected symmetry. In theory, physical quantities should remain consistent when structural and precision parameters are unchanged.

Given these observations and the lack of ZFS information in the VASP wiki, could there be limitations in the current implementation of this feature, which may impact its practical use? Or could there be other factors affecting the results?

I have attached the input files used for my ZFS calculations with VASP. The zero-field splitting matrix for spin-spin interactions can be found using grep MHz OUTCAR.

Thank you very much for your assistance.

[henrique_miranda]

Thank you for your report.

Could you point out the literature references you are referring too where these values are reported?
In particular, the ones where PAW was used too.

[henrique_miranda]

I also noticed two problems in your OUTCAR that might lead to issues in the caluclations:

1. You did not set MAGMOM in the INCAR file.

   Code: Select all

    -----------------------------------------------------------------------------
   |                                                                             |
   |           W    W    AA    RRRRR   N    N  II  N    N   GGGG   !!!           |
   |           W    W   A  A   R    R  NN   N  II  NN   N  G    G  !!!           |
   |           W    W  A    A  R    R  N N  N  II  N N  N  G       !!!           |
   |           W WW W  AAAAAA  RRRRR   N  N N  II  N  N N  G  GGG   !            |
   |           WW  WW  A    A  R   R   N   NN  II  N   NN  G    G                |
   |           W    W  A    A  R    R  N    N  II  N    N   GGGG   !!!           |
   |                                                                             |
   |     You use a magnetic or noncollinear calculation, but did not specify     |
   |     the initial magnetic moment with the MAGMOM tag. Note that a            |
   |     default of 9.302e-3 will be used for all atoms. This ferromagnetic      |
   |     setup may break the symmetry of the crystal, in particular it may       |
   |     rule out finding an antiferromagnetic solution. Thence, we              |
   |     recommend setting the initial magnetic moment manually or verifying     |
   |     carefully that this magnetic setup is desired.                          |
   |                                                                             |
    -----------------------------------------------------------------------------
    

2. You set ISYM=2 but LDMATRIX=.TRUE. activates a Hartree Fock calculation and leads to the following error:

   Code: Select all

    -----------------------------------------------------------------------------
   |                                                                             |
   |           W    W    AA    RRRRR   N    N  II  N    N   GGGG   !!!           |
   |           W    W   A  A   R    R  NN   N  II  NN   N  G    G  !!!           |
   |           W    W  A    A  R    R  N N  N  II  N N  N  G       !!!           |
   |           W WW W  AAAAAA  RRRRR   N  N N  II  N  N N  G  GGG   !            |
   |           WW  WW  A    A  R   R   N   NN  II  N   NN  G    G                |
   |           W    W  A    A  R    R  N    N  II  N    N   GGGG   !!!           |
   |                                                                             |
   |     You have selected ISYM=2 for a HF type calculation. This will           |
   |     symmetrize the charge density but not the exchange potential. I         |
   |     suggest to use ISYM=3 instead. This uses symmetry to obtain the         |
   |     orbitals at all k-points in the Brillouin zone, but does not apply      |
   |     symmetry to the Hartree potential directly. The resultant charge        |
   |     might have lower symmetry than the crystal but, at least, Hartree       |
   |     and exchange are fully compatible.                                      |
   |                                                                             |
    -----------------------------------------------------------------------------

Can you try setting MAGMOM and not setting ISYM=2 in the INCAR file and checking if at least the discrepancy between ISYM=2 and ISYM=0 disappears?

[xin_zhao]

Thank you for your report.

Could you point out the literature references you are referring too where these values are reported?
In particular, the ones where PAW was used too.

Thank you very much for your response.

Regarding the hh and kk NV centers in 4H-SiC, I mainly referred to the article https://doi.org/10.1103/PhysRevMaterials.5.074602. For NV centers in diamond, I referred to https://doi.org/10.1103/PhysRevB.90.235205. Additionally, it is worth noting a review article on NV centers in diamond (https://doi.org/10.1515/nanoph-2019-0154) mentions that "The all-electron version was then implemented by Martijn Marsman into VASP code in 2014." Could you confirm if this refers to the code I am currently using in VASP?

[xin_zhao]

I also noticed two problems in your OUTCAR that might lead to issues in the caluclations:
1. You did not set MAGMOM in the INCAR file.

Code: Select all

 -----------------------------------------------------------------------------
|                                                                             |
|           W    W    AA    RRRRR   N    N  II  N    N   GGGG   !!!           |
|           W    W   A  A   R    R  NN   N  II  NN   N  G    G  !!!           |
|           W    W  A    A  R    R  N N  N  II  N N  N  G       !!!           |
|           W WW W  AAAAAA  RRRRR   N  N N  II  N  N N  G  GGG   !            |
|           WW  WW  A    A  R   R   N   NN  II  N   NN  G    G                |
|           W    W  A    A  R    R  N    N  II  N    N   GGGG   !!!           |
|                                                                             |
|     You use a magnetic or noncollinear calculation, but did not specify     |
|     the initial magnetic moment with the MAGMOM tag. Note that a            |
|     default of 9.302e-3 will be used for all atoms. This ferromagnetic      |
|     setup may break the symmetry of the crystal, in particular it may       |
|     rule out finding an antiferromagnetic solution. Thence, we              |
|     recommend setting the initial magnetic moment manually or verifying     |
|     carefully that this magnetic setup is desired.                          |
|                                                                             |
 -----------------------------------------------------------------------------
 

2. You set ISYM=2 but LDMATRIX=.TRUE. activates a Hartree Fock calculation and leads to the following error:

Code: Select all

 -----------------------------------------------------------------------------
|                                                                             |
|           W    W    AA    RRRRR   N    N  II  N    N   GGGG   !!!           |
|           W    W   A  A   R    R  NN   N  II  NN   N  G    G  !!!           |
|           W    W  A    A  R    R  N N  N  II  N N  N  G       !!!           |
|           W WW W  AAAAAA  RRRRR   N  N N  II  N  N N  G  GGG   !            |
|           WW  WW  A    A  R   R   N   NN  II  N   NN  G    G                |
|           W    W  A    A  R    R  N    N  II  N    N   GGGG   !!!           |
|                                                                             |
|     You have selected ISYM=2 for a HF type calculation. This will           |
|     symmetrize the charge density but not the exchange potential. I         |
|     suggest to use ISYM=3 instead. This uses symmetry to obtain the         |
|     orbitals at all k-points in the Brillouin zone, but does not apply      |
|     symmetry to the Hartree potential directly. The resultant charge        |
|     might have lower symmetry than the crystal but, at least, Hartree       |
|     and exchange are fully compatible.                                      |
|                                                                             |
 -----------------------------------------------------------------------------

Can you try setting MAGMOM and not setting ISYM=2 in the INCAR file and checking if at least the discrepancy between ISYM=2 and ISYM=0 disappears?

Thank you very much for your suggestion.

Based on your advice, I performed tests on the NV center in diamond as a case study, examining various combinations of the MAGMOM parameter and the ISYM tag. Referring to the content on the wiki, I selected ISYM values of -1, 0, 1, 2, and 3, and configured three different MAGMOM settings:

• Without specifying MAGMOM, labeled as "nomag".

• Setting MAGMOM such that all C atoms are 0 and the N atom is 2, labeled as "mag1".

• Setting MAGMOM based on the electronic configuration of elements in the periodic table, where C atoms are 2 and the N atom is 4, labeled as "mag2".

In total, I performed 5×3=15 test calculations. The corresponding INCAR and OUTCAR files are attached, and I summarized the key results in the table below (three values correspond to the eigenvalues obtained from diagonalizing the zero-field splitting matrix).

Based on my results, the effects of the two tags are as follows:

1. Without enabling symmetry (ISYM = -1, 0), the values of the zero-field splitting matrix do not change with different MAGMOM settings. However, these computed values fail to reflect the inherent symmetry of the system (the NV center in diamond exhibits \(C_{3v}\) symmetry, meaning the two negative eigenvalues should be degenerate).

2. With ISYM = 1, the zero-field splitting matrix exhibits the desired symmetry, and the values remain unchanged across different MAGMOM settings. Moreover, the calculated values appear more reasonable compared to those obtained with ISYM = -1, 0.

3. For ISYM = 2 and ISYM = 3, the zero-field splitting matrix values vary with the MAGMOM settings, particularly for ISYM = 3, where the differences are pronounced. When MAGMOM is not specified, the results for these two settings differ significantly, but when MAGMOM is set, the differences become minimal.

4. Under the mag2 condition, enabling symmetry (ISYM = 1, 2, 3) produces identical results.

Summary and Recommendations:
Considering the above, for calculations of the zero-field splitting matrix using VASP, enabling ISYM = 1 appears to be an excellent choice. Alternatively, enabling symmetry with ISYM = 2 or ISYM = 3 and setting MAGMOM based on the number of unpaired electrons from the periodic table also yields reasonable results consistent with ISYM = 1.

However, from a theoretical perspective, I believe the zero-field splitting caused by spin-spin interactions should not change when the computational precision parameters are kept constant, regardless of other tag settings. This has been confirmed by our group's in-house program, which calculates zero-field splitting matrices using VASP's WAVECAR. At least for ISYM = -1, 2, 3 with MAGMOM unset, our program consistently produces identical results. Thus, I suspect the lack of robustness in VASP's zero-field calculation might be due to certain approximations made to save computational resources when handling symmetry internally.

I plan to test the NV center in SiC to verify whether the conclusions drawn here are consistent.

Finally, thank you again for your suggestion, which has been immensely helpful to me.

[henrique_miranda]

Yes, the implementation of this feature (LDMATRIX=.TRUE.) was done by Martijn Marsman in VASP.

Thank you for your detailed tests.
I should have been clearer about my suggestion to not set ISYM in the INCAR.
The point is that some settings of ISYM will not produce the correct results when used in combination with LHFCALC which is internally set to TRUE when LDMATRIX=.TRUE.
That being said, I do find it very strange that you obtain different values for ISYM=-1 and ISYM=3 for these calculations.
I would not expect that, so I tried running the same calculations as you locally and got agreement between ISYM=-1 and ISYM=3.

Can it be that something has been modified in your version of VASP?
Can you try re-runing the calculations with ISYM=-1 or without ISYM in the INCAR (which is the same as ISYM=3) and removing ISTART=1?

[xin_zhao]

Thank you very much for your response.

I will continue to test the setup, but I would like to highlight an issue I recently discovered. The computational platform I am currently using consists of two types of CPUs: Intel Xeon Scale 8358 and AMD EPYC 9654. To optimize performance, I use the OneAPI compiler to build VASP on the Intel platform and AOCC+AOCL compilers for the AMD platform. However, a significant issue has arisen: the calculation results differ between the two platforms under identical tag inputs. More concerning, when running the same calculations repeatedly on the AMD platform (with WAVECAR and CHGCAR removed to ensure identical initial conditions), the results are inconsistent across runs.

I noticed that your tests seem to be based on VASP compiled with the GNU toolchain. Could differences in the compilation environment and methodology be a potential reason for these discrepancies?

Additionally, I have made some modifications to VASP. Specifically, I am incorporating spin contamination corrections to the zero-field matrix as described in https://doi.org/10.1103/PhysRevResearch.2.022024. To achieve this, I calculate the zero-field matrix for the state with \(S=0\). Since VASP’s default settings skip this state, I modified dmatrix.F as follows:

from

Code: Select all

IF (.NOT.LDO_DMATRIX) RETURN

STOT=ABS(RHO0(GRIDC,CHTOT(1,2)))/2._q
IF (STOT<0.5_q) THEN
   IF (IU0>=0) THEN
      WRITE(IU0,'(A,F14.7,A)') 'DMATRIX_CALC: ERROR: S=',STOT,' < 1/2, skipping calculation of D-matrix ...'
   ENDIF
   RETURN
ENDIF

DMAT=0._q

! CALL DMATRIX_CALCULATE_JIJ_ALT(GRIDC,LATT_CUR,CHTOT(1,2),IU6,IU0)
CALL DMATRIX_CALCULATE_JIJ(WDES,W,P,LATT_CUR,LMDIM,NONL_S,NONLR_S,GRID,IU6,IU0)
CALL DMATRIX_CALCULATE_KIJ(WDES,W,P,LATT_CUR,NONL_S,NONLR_S,GRID,IU6,IU0)
CALL SET_DAB_PAW(WDES,T_INFO,P,LMDIM,CRHODE,IU6,IU0)

DMAT=DMAT/STOT/(2*STOT-1)

to

Code: Select all

REAL(q) STOT
REAL(q), EXTERNAL :: RHO0

IF (.NOT.LDO_DMATRIX) RETURN

! STOT=ABS(RHO0(GRIDC,CHTOT(1,2)))/2._q
! IF (STOT<0.5_q) THEN
!    IF (IU0>=0) THEN
!       WRITE(IU0,'(A,F14.7,A)') 'DMATRIX_CALC: ERROR: S=',STOT,' < 1/2, skipping calculation of D-matrix ...'
!    ENDIF
!    RETURN
! ENDIF

DMAT=0._q

! CALL DMATRIX_CALCULATE_JIJ_ALT(GRIDC,LATT_CUR,CHTOT(1,2),IU6,IU0)
CALL DMATRIX_CALCULATE_JIJ(WDES,W,P,LATT_CUR,LMDIM,NONL_S,NONLR_S,GRID,IU6,IU0)
CALL DMATRIX_CALCULATE_KIJ(WDES,W,P,LATT_CUR,NONL_S,NONLR_S,GRID,IU6,IU0)
CALL SET_DAB_PAW(WDES,T_INFO,P,LMDIM,CRHODE,IU6,IU0)

DMAT=DMAT

These changes should not affect the zero-field matrix for states with \(S=1\)

[henrique_miranda]

Ok, I tried comparing your OUTCARs with mine, and at first I did not see anything that could explain the differences you observe for ZFS.
I also tested with a version compiled using OneAPI.

Today I realized that you are using the vasp_gam version, while I used vasp_std.

I tried running with vasp_gam and sure enough I reproduce the issue you observed!
This is likely a bug, but first we need to check.

In the meantime, I would suggest that you try using the vasp_std version.

[xin_zhao]

Yes, you are right. I also find that vasp_std gives reasonable results, which is not depend on ISYM and MAGMOM tag. I checked the log file and found that vasp_gam skips calculating one integral compared to vasp_std during computation, which might be the source of the problem.

the vasp_gam log contains

Code: Select all

DAV:  13    -0.204926783124E+04   -0.92596E-04   -0.15218E-04  2656   0.215E-02    0.315E-02
DAV:  14    -0.204926785913E+04   -0.27891E-04   -0.31293E-05  2656   0.100E-02    0.129E-02
DAV:  15    -0.204926785909E+04    0.43518E-07   -0.55381E-06  2464   0.395E-03
 isp1:  1  isp2:  1  IDO:     2985984
 isp1:  1  isp2:  2  IDO:     2972160
 isp1:  2  isp2:  2  IDO:     2958400
 Jij:  2110.7986225   126.0136583 -2236.8121450  1370.3395847   964.6389002   770.6071336
 isp1:  1  isp2:  1  IDO:     2985984
 isp1:  1  isp2:  2  IDO:     2972160
 isp1:  2  isp2:  2  IDO:     2958400
-Kij: -1981.4570479  -217.4454315  2198.9022454 -2285.6204957    58.0288852    99.8881261
 D1c:    -0.0000121    -0.0000047     0.0000168    11.4488848    11.4487759    11.4488373
   1 F= -.20492679E+04 E0= -.20492679E+04  d E =-.805274E-14  mag=     2.0000
 writing wavefunctions

and the vasp_std log contains

Code: Select all

DAV:  13    -0.204926764686E+04   -0.60327E-04   -0.59639E-04  2720   0.389E-02    0.116E-01
DAV:  14    -0.204926779616E+04   -0.14931E-03   -0.31674E-04  2720   0.313E-02    0.503E-02
DAV:  15    -0.204926785050E+04   -0.54339E-04   -0.14515E-04  2624   0.192E-02    0.303E-02
DAV:  16    -0.204926785737E+04   -0.68661E-05   -0.14076E-05  2624   0.520E-03    0.200E-02
DAV:  17    -0.204926785787E+04   -0.50088E-06   -0.94431E-06  2688   0.558E-03
 isp1:  1  isp2:  1  IDO:      1492992
 isp1:  1  isp2:  2  IDO:      1486080
 isp1:  2  isp2:  1  IDO:      1486080
 isp1:  2  isp2:  2  IDO:      1479200
 Jij:     0.0020561     0.0009124    -0.0029684   971.9737531   971.9717022   971.9738489
 isp1:  1  isp2:  1  IDO:      1492992
 isp1:  1  isp2:  2  IDO:      1486080
 isp1:  2  isp2:  1  IDO:      1486080
 isp1:  2  isp2:  2  IDO:      1479200
-Kij:    -0.0000968    -0.0008064     0.0009033    58.0110791    58.0106129    58.0099581
 D1c:     0.0000477     0.0000225    -0.0000702    11.4499716    11.4500132    11.4499822
   1 F= -.20492679E+04 E0= -.20492679E+04  d E =-.619615E-14  mag=     2.0000
 writing wavefunctions

It seems that the (2,1) pair contribution is ignored in vasp_gam.

Looking back, this seems to be a fairly common issue. If I had included the specific version of the VASP program I used when describing the problem, the issue might have been resolved more quickly. Anyway, big thanks for your suggestions; I apologize for taking up so much of your time on this issue, but I deeply appreciate the effort you have dedicated to it.