materialsproject / pymatgen

Python Materials Genomics (pymatgen) is a robust materials analysis code that defines classes for structures and molecules with support for many electronic structure codes. It powers the Materials Project.
https://pymatgen.org
Other
1.53k stars 868 forks source link

Incorrect Gvec generation in Wavecar class of pymatgen.io.vasp.outputs module #2361

Closed RyanNKHall closed 2 years ago

RyanNKHall commented 2 years ago

Describe the bug Within the pymatgen.io.vasp.outputs module's Wavecar class is a function called _generate_G_points (line 5012) that generates the G vectors for the plane wave expansion of the wavefunction using the cutoff energy from the Wavecar file and the kpoints. This function fails to generate the correct number of Gpoints when tested against a range of VASP version 5 Wavecar files. I found the issue arises due to the _encut attribute, which represents the cutoff energy for the plane wave expansion being cast to an integer type in the Wavecar class init function (line 4833). I was able to fix the bug by leaving encut as a float.

To Reproduce Steps to reproduce the behavior:

  1. Use the example WAVECAR file provided. It was generated for Silicon (mp149 in materials project) with hybrid functionals and a Monkhorst-pack grid of 9 x 9 x 9. The resultant Wavecar file contains 35 kpoints, 96 bands and has an encut of 245.345 when not cast to an integer -> this is cast to 245 in the original code
  2. I have provided a test Python script called Wavefunction.py that creates an instance of the Wavecar class with the WAVECAR file to reproduce the error
  3. A value error "Incorrect value of vasp_type given..." is produced on the third kpoint (ink = 2) (the casting/rounding does not affect the first two iterations likely because the kpoint magnitude is small)
  4. Adding a print statement to outputs.py if statement that raises the value error provides some more verbose information on the issue. Specifically, adding: print("\n# Gpoints:", len(self.Gpoints[ink]), "# plane waves:", nplane, "# extra Gpts:", len(extra_gpoints), "\n") before the value error is raises will display to the terminal the number of Gpoints generated versus the required number nplane. In this example, 244 Gpoints are generated instead of the expected 253.
  5. If you modify line 4833 of outputs.py from:
    # extract kpoint, bands, energy, and lattice information
    self.nk, self.nb, self.encut = np.fromfile(f, dtype=np.float64, count=3).astype(np.int_)

    to:

    self.nk, self.nb = np.fromfile(f, dtype=np.float64, count=2).astype(np.int_)
    self.encut = np.fromfile(f, dtype=np.float64, count=1)

    The error/issue is resolved. I have not tested this with earlier VASP WAVECAR outputs

Provide any example files that are needed to reproduce the error, especially if the bug pertains to parsing of a file.

Expected behavior The _generate_G_points should generate the correct number of Gpoints for the specified VASP WAVECAR type e.g. std, gam, or ncl

Files required to reproduce error bug_report.zip

Desktop (please complete the following information):

Additional context Feel free to contact me via email: ryan.n.hall@unsw.edu.au

Andrew-S-Rosen commented 2 years ago

@RyanNKHall: Just a head's up, the line you want is probably

self.encut = np.fromfile(f, dtype=np.float64, count=1)[0]

Without the [0] indexing, self.encut is returned as a length-one array instead of a float.

RyanNKHall commented 2 years ago

Thanks @arosen93. I actually spotted this the other day too. It doesn't seem to cause any issues without it because I think most functions that use encut would be able to unpack the value from a size 1 array but I guess it could cause issues and it's more correct to reference the value.

I noticed there's a variable that's unpacked a few lines down that indexes the array in this way and it's obvious now that a numpy function would create an array by default.

Andrew-S-Rosen commented 2 years ago

This is now addressed via #2410. @RyanNKHall or @mkhorton, looks like you can close this one.

RyanNKHall commented 2 years ago

Addressed via #2410 Thanks @arosen93 for completing the PR and @mkhorton for the support.