My Community

I think that I found a bug in the write out of matrix elements to at least the WAVEDERF file (and perhaps WAVEDER, but I don't know) when using different settings of NCORE/NPAR. I am currently using VASP 6.2.1 for these calculations compiled with Intel 18.0.1.163 and Intel MKL on NERSC's Cori supercomputer. The only modification to the VASP code that I made was to uncomment the line in linear_optics.F so that the WAVEDERF file would be written out. That way I can use a postprocessing code to recalculate the dielectric function or extract specific matrix elements for doing excitation rate and lifetime calculations.

As a test to play with dipole transition matrix elements, I was doing a relatively simple calculation of the primitive diamond unit cell with an 8x8x8 gamma centered grid (29 k-points) and NBANDS = 32 using the PBE functional and LOPTICS = .TRUE.

When I have NPAR = 1, I get sensible results for the dielectric function and the real part's limit is approximately equal (within convergence/numerics) to the value obtained through DFPT. Furthermore, running my python postprocessing script on the WAVEDERF file yields the same result.
The only issue that I have is that in the WAVEDERF, only information for matrix elements up to band 8 are outputted.

When I set NCORE to other values than the equivalent of NPAR = 1, for example NCORE = 2 or NCORE = 4, the density-density version of the dielectric function remains the same in the OUTCAR as before. However, the values of the matrix elements are no longer the same as those in the above file and I get different results with the postprocessing code.
In this case, all band pairs are included in the WAVEDERF though (despite the matrix elements in WAVEDERF being different than expected).

Below are the files I used for inputs. I can also provide any other information about the runs as needed.

As an aside, I have also run similar tests on a different supercomputer that I have access to and had similar issues as reported above with the values of the matrix elements in WAVEDERF changing depending upon the setting of NCORE/NPAR.

Thank you for any help and advice,
Andy


--------------
KPOINTS
--------------
Automatic mesh
0
Gamma
8 8 8
0 0 0


--------------
POSCAR
--------------
Diamond
3.574131
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.0000000000 0.5000000000
0.5000000000 0.5000000000 0.0000000000
C
2
Direct
0.000000000 0.000000000 0.000000000
0.250000000 0.250000000 0.250000000


--------------
INCAR
--------------
SYSTEM = Diamond

ISTART = 0
ICHARG = 2
INIWAV = 1

NBANDS = 32

ENCUT = 450
EDIFF = 1E-5
ALGO = Fast
NELMIN = 2
NELM = 80
LREAL = .FALSE.
ISPIN = 1
LOPTICS = .TRUE.

ISMEAR = 0
SIGMA = 0.05
NEDOS = 5000

NPAR = 1
!NCORE = 2
!NCORE = 32

After doing a little more testing, I wanted to provide an update. My problem is partially solved (i.e. re-calculating the dielectric function is solved, but I have some questions).

After working over my python script again and re-reading the linear_optics.F source, I realized there was a unit / scaling factor issue in my code. My code originally worked on the OPTIC file from optics.F along the lines of the Fortran code originally by Dr. Juergen Furthmueller. After realizing that we are now reading in -<u_m|r_j|u_n> rather than the nabla operator matrix element, I was able to properly update my code to normalize terms correctly and I do obtain within numerical accuracy the same that comes out of the OUTCAR file, regardless of NCORE/NPAR.

However, I am still curious about the dependence of the entries in WAVEDERF on NCORE/NPAR. Again, a little more checking through the different WAVEDERF files that I generated indicated that the modulus of the matrix elements is the same. That is |<u_m|r_j|u_n>|^2 is the same no matter which NCORE setting I use. So this means that there is only a phase factor between the different elements. While it is obvious to me that this won't affect diagonal elements of the dielectric tensor, I'm not sure about the off diagonal terms in systems with lower than orthorhombic symmetry.

I don't understand why there is such a phase factor and would appreciate any insight that can be given.

Further, I do still have an issue that for certain settings of NCORE rather than obtaining matrix elements between all 32 bands in the WAVEDERF file, I only obtain a subset even if NBANDS = 32 and the number of entries in EIGENVAL is also 32. For example, using 32 cores and 32 bands, NCORE = 1 and 2 yield a full WAVEDERF file. For NCORE =4, only matrix elements for bands 1 through 16 are written out. For NCORE = 8, only matrix elements for bands 1 through 8 are written out.

Thanks again for any help.

Regarding the phase factor of the matrix elements:
The eigenvectors belonging to a degenerate subspace are defined up to a phase factor so you might get slightly different results in different runs. This should however not affect your final results.
When you say that you observe that the matrix elements are the same up to a phase factor do you mean this for all the states or just for states in degenerate subspaces?

Regarding the incomplete WAVEDERF files:
If you want all the bands to be written you should set LVEL=.TRUE. in the INCAR file.

Let me know if this answers your question.