diff --git a/src/lobsterpy/cohp/analyze.py b/src/lobsterpy/cohp/analyze.py index 7d649396..8c61bde5 100644 --- a/src/lobsterpy/cohp/analyze.py +++ b/src/lobsterpy/cohp/analyze.py @@ -15,7 +15,7 @@ from pymatgen.core.structure import Structure from pymatgen.electronic_structure.cohp import CompleteCohp from pymatgen.electronic_structure.core import Spin -from pymatgen.electronic_structure.dos import LobsterCompleteDos +from pymatgen.electronic_structure.dos import CompleteDos, LobsterCompleteDos from pymatgen.io.lobster import ( Bandoverlaps, Charge, @@ -1450,6 +1450,7 @@ def get_lobster_calc_quality_summary( bandoverlaps_obj: Bandoverlaps | None = None, lobster_completedos_obj: LobsterCompleteDos | None = None, vasprun_obj: Vasprun | None = None, + ref_dos_obj: CompleteDos | None = None, dos_comparison: bool = False, e_range: list = [-5, 0], n_bins: int | None = None, @@ -1474,6 +1475,7 @@ def get_lobster_calc_quality_summary( :param bandoverlaps_obj: pymatgen lobster.io.BandOverlaps object :param lobster_completedos_obj: pymatgen.electronic_structure.dos.LobsterCompleteDos object :param vasprun_obj: pymatgen vasp.io.Vasprun object + :param ref_dos_obj: pymatgen.electronic_structure.dos.CompleteDos object :param dos_comparison: will compare DOS from VASP and LOBSTER and return tanimoto index :param e_range: energy range for DOS comparisons :param n_bins: number of bins to discretize DOS for comparisons @@ -1635,7 +1637,7 @@ def get_lobster_calc_quality_summary( stacklevel=2, ) if dos_comparison: - if "LSO" not in str(path_to_doscar).split("."): + if "LSO" not in str(path_to_doscar).split(".") and lobster_completedos_obj is None: warnings.warn( "Consider using DOSCAR.LSO.lobster, as non LSO DOS from LOBSTER can have negative DOS values", stacklevel=2, @@ -1655,63 +1657,68 @@ def get_lobster_calc_quality_summary( "Dos comparison is requested, so please provide either path_to_doscar or lobster_completedos_obj" ) - if path_to_vasprun: - vasprun = Vasprun(path_to_vasprun, parse_potcar_file=False, parse_eigen=False) - elif vasprun_obj: - vasprun = vasprun_obj + if path_to_vasprun and not ref_dos_obj: + dos_vasp = Vasprun(path_to_vasprun, parse_potcar_file=False, parse_eigen=False).complete_dos + elif vasprun_obj and not ref_dos_obj: + dos_vasp = vasprun_obj.complete_dos + elif ref_dos_obj: + dos_vasp = ref_dos_obj else: raise ValueError( - "Dos comparison is requested, so please provide either path to vasprun.xml or vasprun_obj" + "Dos comparison is requested, so please provide either path to vasprun.xml or " + "vasprun_obj or ref_dos_obj" ) - dos_vasp = vasprun.complete_dos quality_dict["dos_comparisons"] = {} # type: ignore - for orb in dos_lobster.get_spd_dos(): - if e_range[0] >= min(dos_vasp.energies) and e_range[0] >= min(dos_lobster.energies): - min_e = e_range[0] - else: - warnings.warn( - "Minimum energy range requested for DOS comparisons is not available " - "in VASP or LOBSTER calculation. Thus, setting min_e to -5 eV", - stacklevel=2, - ) - min_e = -5 + min_e = int(max(e_range[0], min(dos_vasp.energies), min(dos_lobster.energies))) + max_e = int(min(e_range[-1], max(dos_vasp.energies), max(dos_lobster.energies))) - if e_range[-1] <= max(dos_vasp.energies) and e_range[-1] <= max(dos_lobster.energies): - max_e = e_range[-1] - else: - warnings.warn( - "Maximum energy range requested for DOS comparisons is not available " - "in VASP or LOBSTER calculation. Thus, setting max_e to 0 eV", - stacklevel=2, - ) - max_e = 0 - - if np.diff(dos_vasp.energies)[0] >= 0.1 or np.diff(dos_lobster.energies)[0] >= 0.1: - warnings.warn( - "Input DOS files have very few points in the energy interval and thus " - "comparisons will not be reliable. Please rerun the calculations with " - "higher number of DOS points. Set NEDOS and COHPSteps tags to >= 2000 in VASP and LOBSTER " - "calculations, respectively.", - stacklevel=2, - ) + if min_e > e_range[0]: + warnings.warn( + f"Minimum energy range requested for DOS comparisons is not available in " + "VASP or LOBSTER calculation. " + f"Thus, setting `min_e` to the minimum possible value of {min_e} eV", + stacklevel=2, + ) + if max_e < e_range[-1]: + warnings.warn( + f"Maximum energy range requested for DOS comparisons is not available in " + "VASP or LOBSTER calculation. " + f"Thus, setting `max_e` to the maximum possible value of {max_e} eV", + stacklevel=2, + ) + + minimum_n_bins = min( + len(dos_vasp.energies[(dos_vasp.energies >= min_e) & (dos_vasp.energies <= max_e)]), + len(dos_lobster.energies[(dos_lobster.energies >= min_e) & (dos_lobster.energies <= max_e)]), + ) + + n_bins = n_bins or minimum_n_bins + + if n_bins > minimum_n_bins: + warnings.warn( + f"Number of bins requested for DOS comparisons is larger than the " + "number of points in the energy interval. Thus, setting " + f"`n_bins` to {minimum_n_bins}.", + stacklevel=2, + ) + n_bins = minimum_n_bins - if not n_bins: - n_bins = 56 + dos_fp_kwargs = { + "min_e": min_e, + "max_e": max_e, + "n_bins": n_bins, + "normalize": True, + } + for orb in dos_lobster.get_spd_dos(): fp_lobster_orb = dos_lobster.get_dos_fp( - min_e=min_e, - max_e=max_e, - n_bins=n_bins, - normalize=True, + **dos_fp_kwargs, fp_type=orb.name, ) fp_vasp_orb = dos_vasp.get_dos_fp( - min_e=min_e, - max_e=max_e, - n_bins=n_bins, - normalize=True, + **dos_fp_kwargs, fp_type=orb.name, ) @@ -1722,17 +1729,11 @@ def get_lobster_calc_quality_summary( quality_dict["dos_comparisons"][f"tanimoto_orb_{orb.name}"] = tani_orb # type: ignore fp_lobster = dos_lobster.get_dos_fp( - min_e=min_e, - max_e=max_e, - n_bins=n_bins, - normalize=True, + **dos_fp_kwargs, fp_type="summed_pdos", ) fp_vasp = dos_vasp.get_dos_fp( - min_e=min_e, - max_e=max_e, - n_bins=n_bins, - normalize=True, + **dos_fp_kwargs, fp_type="summed_pdos", ) diff --git a/tests/cli/test_cli.py b/tests/cli/test_cli.py index 29c27be4..67789a5d 100644 --- a/tests/cli/test_cli.py +++ b/tests/cli/test_cli.py @@ -553,7 +553,7 @@ def test_calc_quality_summary_nacl(self, tmp_path): "The atomic charge signs from Mulliken population analysis agree with the bond valence analysis. " "The atomic charge signs from Loewdin population analysis agree with the bond valence analysis. " "The Tanimoto index from DOS comparisons in the energy range between -20, 0 eV for s, p, summed orbitals " - "are: 0.9935, 0.9983, 0.9822." + "are: 0.9795, 0.9864, 0.9628." ) assert calc_quality_text == ref_text @@ -595,7 +595,7 @@ def test_calc_quality_summary_k3sb(self, tmp_path): "The atomic charge signs from Mulliken population analysis agree with the bond valence analysis. " "The atomic charge signs from Loewdin population analysis agree with the bond valence analysis. " "The Tanimoto index from DOS comparisons in the energy range between -20, 0 eV for s, p, summed orbitals " - "are: 0.8367, 0.9565, 0.9357." + "are: 0.7953, 0.8994, 0.8763." ) assert calc_quality_text == ref_text diff --git a/tests/cohp/test_describe.py b/tests/cohp/test_describe.py index d1c29933..455630ab 100644 --- a/tests/cohp/test_describe.py +++ b/tests/cohp/test_describe.py @@ -362,28 +362,26 @@ def test_warnings(self): e_range=[-50, 60], dos_comparison=True, bva_comp=False, + n_bins=500, ) - messages = [] - for warning in w: - messages.append(str(warning.message)) - count0 = 0 - count1 = 0 - count2 = 0 - count3 = 0 - for msg in messages: - if "Consider using DOSCAR.LSO.lobster" in msg: - count0 += 1 - if "Minimum energy range requested" in msg: - count1 += 1 - if "Maximum energy range requested" in msg: - count2 += 1 - if "Input DOS files have very few points" in msg: - count3 += 1 + # messages = [] + actual_warnings = [str(warning.message) for warning in w] - assert count0 == 1 - assert count1 == 1 - assert count2 == 1 - assert count3 == 1 + expected_warnings = [ + "Consider using DOSCAR.LSO.lobster, as non LSO DOS from LOBSTER can have negative DOS values", + "Minimum energy range requested for DOS comparisons is not available in VASP or LOBSTER calculation. " + "Thus, setting `min_e` to the minimum possible value of -15 eV", + "Maximum energy range requested for DOS comparisons is not available in VASP or LOBSTER calculation. " + "Thus, setting `max_e` to the maximum possible value of 5 eV", + "Number of bins requested for DOS comparisons is larger than the number of points in the energy interval. " + "Thus, setting `n_bins` to 107.", + "Input DOS files have very few points in the energy interval and thus comparisons will not be reliable. " + "Please rerun the calculations with higher number of DOS points. " + "Set NEDOS and COHPSteps tags to >= 2000 in VASP and LOBSTER calculations, respectively.", + ] + + for actual, expected in zip(actual_warnings, expected_warnings): + assert actual == expected calc_des = Description.get_calc_quality_description(calc_quality_warnings) @@ -392,8 +390,8 @@ def test_warnings(self): "The absolute and total charge spilling for the calculation is 2.255 and 12.72 %, respectively.", "The projected wave function is completely orthonormalized as no bandOverlaps.lobster file is " "generated during the LOBSTER run.", - "The Tanimoto index from DOS comparisons in the energy range between -5, 0 eV for s, p, summed orbitals " - "are: 0.4057, 0.2831, 0.2762.", + "The Tanimoto index from DOS comparisons in the energy range between -15, 5 eV for s, p, summed orbitals " + "are: 0.6712, 0.8113, 0.8064.", ] with warnings.catch_warnings(record=True) as w2: