Marginal agreement with C_N Wif value reported in the tutorial

[Foresttrek]

I ran the tutorial for C_N, but the Wif value for the dominant valence band overlap (#189) is about 15-20% larger than seen in the tutorial (Wif=0.0978, vs Wif=0.0808 in the tutorial). I am using VASP 6.4, but the original work seems to have been done using VASP 5.4.4. Should VASP 6.4 reproduce the VASP 5.4 results for Wif to within a percent or less, or is the larger difference that I see just a result of using VASP 6 in place of VASP 5?

[Foresttrek]

I should add that your INCARs, POSCARs, KPOINTS, and POTCAR files posted at Nomad and GitHub were used for the WSWQ calculations. So, the different results that I find should not be because of any newer version POTCARs that one could use with VASP 6.4.

Thanks for any insights or comments.

[mturiansky]

If all the calculation details are exactly the same, you should get exactly the same number independent of the version number. That being said, there could be calculation details that are different that were overlooked. For example, VASP 6 has ACE implemented for the Fock operator, which is in principle a difference with respect to VASP 5. There may be other differences. That said, are you certain you have the same pseudopotentials? The VASP 5 pseudopotentials had several variants floating around (some without any documentation of the change). I don't believe the actual POTCAR was directly loaded anywhere since that's not allowed.

Ultimately, a 15-20% error in the matrix element is not a huge deal compared to other parameters. A 20% error in the dE, dQ, or phonon frequencies would be a big problem.

[Foresttrek]

Mark, regarding your pseudopotentials as seen at NOMAD for your work on NonRad, there are two versions posted there that are named POTCAR and POTCAR.stripped. I am assuming that the stripped version was created to let others see the just POTCAR header information and thus verify their POTCAR against yours, but without posting the full set of pseudo-potential data. If the intent was to then delete the unedited complete POSCAR, it has accidentally been left posted, and you may wish to delete it.

Based on your comments, I went back and recalculated WAVECAR and WSWQ files with ACE disabled in VASP-6.4, and that has no significant effect (of order 1% change in W_if) on the 20-ish% difference in my W_if results with yours. Because I do reproduce your groundstate total energies for q=-1 and q=0 substitutional C, and I also get your results for delta-Q between the ground states, it seems likely that the differing W_if that I get originates from the WSWQ files generated by my use of VASP 6.4.2.

While I agree that VASP-5.4.4 and VASP-6.4.2 "should" give the same numerical results for the same inputs (accounting for differing defaults), I was wondering if you have ever directly tested this to see if it is true for the WSWQ files and W_if (and not just the defect total energies)? Maybe the overlap routine has been modified since VASP-5.4.4, as that would explain why the energies fully reproduce, but W_if does so only approximately. Perhaps your collaborator on NonRad (Prof. Kresse) would know.

While the only somewhat differing W_if values are not themselves critical, I was hoping to understand the source of irreproducibility.

Thanks again.

[mturiansky]

Thanks for pointing out the issue with Nomad. I contacted their administrators, and they confirmed that there was a bug and will be removing the file.

There was a bug fixed going from VASP 5 to 6 (described here), but this would not explain a 15-20% error. That produced nonsense results. Actually, this is more like the difference between kpoints. Do you append the gamma point as a 0-weight kpoint before calculating the WSWQ's? I didn't mention that in the tutorial, but it's mentioned in the Comput Phys Commun paper.

[Foresttrek]

Regarding the KPOINTS files, I do append a 0-weight gamma point in addition to the unity-weight Baldereschi point. However, I append it for BOTH the construction of the WAVECARs for the interpolated structures and for the subsequent calculation of WSWQ files. Do you instead append the gamma point only during the WSWQ calculation? (Seems like maybe this should not matter, but it could be a k-point-related difference in what I have done.)

[mturiansky]

Did you specify kpoint=-1 or kpoint=2 so that it reads the last kpoint rather than the first when running get_Wif_from_WSWQ?

[Foresttrek]

Each k-points was tried with each spin prior to posting here, just to see what happens; kpoint=2 and the second spin channel (spin down) gave the best results (these are the ones that are 15-20% more than yours).

[mturiansky]

Can you please provide the plot from passing the fig kwarg? Maybe also the INCAR, KPOINTS, POSCAR, and vasprun.xml for a couple of the CCD and WSWQ calculations? I am inclined to believe that this is not related to something changing between VASP 5 and 6. Certainly early versions of VASP 6 reproduced the results. If I can't find an issue in the figure or those files, then maybe I can try to re-run the calculations myself.

[Foresttrek]

Here are the W_if and V*W_if^2 values that I get in my most recent attempts to reproduce the tutorial's W_if values:

Wifs
[(np.int64(189), np.float64(0.10074087128409234)), (np.int64(190), np.float64(0.005792509079581725)), (np.int64(191), np.float64(0.014800031327711622))]
np.mean([x[1]**2 for x in Wifs])*ground_struct.volume
np.float64(3.8217052254022974)

Here is the figure made by get_Wif_from_WSWQ, where the interpolation parameter step size is 0.05 (5%) and the middle of the plot corresponds to the neutral-defect ground state: (please let me know if this is not the figure you intend)

I will send the related VASP files separately.

[Foresttrek]

Here are the files requested, attached below

files_to_MET_250813.tar.gz

as a gzipped tarfile. There are two sets of files, one where the WAVECARs were made, and one where the WSWQ files were made. Within each set there are four examples, consisting of the two ground states and the two extremal Q values that were used to calculate W_if values. The filename/directory suffixes of the style "_pX.XX" indicate the underlying interpolation parameter, where p0.00 corresponds to the q=-1 ground state and p1.00 corresponds to the q=0 ground state. I have also included four actual WSWQ files. The WSWQ files for p = 0.75, 1.00, and 1.25 "should" give approximately W_if=0.097 to 0.107 for band 189 (instead of the tutorial's reported 0.0808), if you were to try them out in your NonRad install and the calculation worked out similarly to here. Maybe see if you differ from me for exactly the same WSWQ files?.

[Foresttrek]

To prevent confusion, I should mention that my POSCARs were made using a pre-existing PERL script that predates my use of NonRad, which I mainly became interested in because of the get_Wif_from_WSWQ capability. So the POSCARs will perhaps look slightly different than expected. The delta-Q that the PERL script reports for your two fully relaxed C_N ground-state POSCARs is dQ=1.6859, while your 2021 paper reports 1.69. So, I am assuming the interpolated POSCARs are likely to be good.

[mturiansky]

That is the figure that I intended, thanks.

One thing I notice is that I think you're plotting over a larger range of Q than did. I believe I only went up to 20% of dQ (0.338 amu^{1/2} Ang). Band 189 is a little non-linear, and I wonder if this contributes partially to the discrepancy.

I also notice that in the vasprun.xml of your make_WSWQ_* directories, the ALGO = none seems to be ignored. The IALGO corresponds to ALGO = Normal. I think you might need to write "None", on the wiki it says: "Tip: Except for 'None', 'Nothing', 'Exact' and 'Eigenval' (which must be spelled out), the first letter determines the applied algorithm.". Having ALGO doubly defined could also confuse it. In any case, I don't think this should be a major source of discrepancy because the wavefunctions should only be getting "more accurate". It is, however, a potential difference. When you set ALGO = None, you should also set NELM = 1, since each electronic step isn't doing anything.

The dQ value doesn't really tell me anything about the interpolation that the script is doing, since that is determined without doing any interpolation. I assume it's just linearly interpolating the positions of each site, right?

[Foresttrek]

Mark, The range of Q is a completely minor problem so long as both ranges (yours and mine) are centered on the minimum point of the parabola (as they are), in which case symmetrically changing the range of Q tends to move the fitted slope up or down a bit, but with almost no change in slope. Also, the ALGO issue is similarly immaterial, interestingly though, ALGO=none (or even "None") gives a nonsensical total energy [whereas ALGO=normal does not]. While the ALGO=none restart's total energy is wrong, the restored/restarted wavefunctions and WSWQ results remain completely fine and essential match the results given using by ALGO=normal at the make WSWQ step.

Since these suggestions did not solve the reproducibility problem, I spent a great deal of time revisiting and analyzing (1) all my files intended to mimic yours, (2) all your similar available-online files (Nomad and GitHub), (3) various NonRad output and plots from the tutorial, and importantly (4) the archived older versions of the NonRad code and its tutorials. I won't retrace my messy path through all of this effort, but I did figure out exactly why your original results (as reported in your paper and in the tutorial) cannot be reproduced using NonRad's latest version.

Basically, NonRad V1.2.0 and your paper-recommended k-point sampling can NEVER reproduce the original paper/tutorial results for two simple reasons:
(1) The paper and the tutorial were done using NonRad V1.0.0, which contains the significant and KNOWN-TO-YOU eigenvalue error problem, which is fixed in all later code versions.
(2) The paper and the tutorial effectively used only Baldereschi kpoint sampling for all steps, and unlike statements in the NonRad paper, the unweighted gamma-point WAS NOT used in conjunction with the unity-weighted Baldereschi point to obtain WSWQ files. This is seen in the tutorial text where get_Wif_from_WSWQ is invoked. At this step kpoint=2 needs to be invoked but you instead used: " Wifs = get_Wif_from_WSWQ(WSWQs, str(ground_files / 'vasprun.xml'), 192, [189, 190, 191], spin=1, fig=fig) ". Thus kpoint=2 (the zero-weight gamma point in KPOINTS), was ignored by NonRad, and the default (kpoint=1) was used instead such the WSWQ files were parsed by NonRad using the Baldereschi-based WSWQ values.

To fully test the above assertions, I restored the eigenvalue error to NonRad V1.2.0 and reran my implementations of your input files using kpoint=1 at the get_Wif_from_WSWQ step when using NonRad. With these two simple changes I immediately reproduce the tutorial results -- QED.

While I appreciate your efforts to make NonRad available to the physics community, you should make it completely clear to prospective users that the tutorial cannot be reproduced for the reasons noted above (and also because PAW versions may differ).

Alternatively, you could redo the tutorial to give correct results when using the current versions of NonRad, VASP, and VASP PAWs. Sure, the redone results will differ from your paper and will reveal to some few that you published an error, but as you said at the outset herein, the moderate differences are immaterial (to the paper's conclusions, but not to your Users' time).

Thanks for your volunteered time spent corresponding with the physics community, Steve.

[mturiansky]

I'm very sorry for the significant amount of time you spent on this. Thank you for getting to the bottom of it. You're referring to a bug that was fixed almost five years ago, and I haven't thought about for equally as long. It didn't come to mind during our discussion but clearly should have. I should have updated the tutorial back when I fixed the bug, but that didn't occur to me at the time (which is silly in hindsight). Again, I apologize for the time you sank into this. I will update the tutorial.