Merge pull request #60 from SMTG-UCL/develop

Develop
SMTG-Bham · Sep 10, 2023 · 31c02ed · 31c02ed
2 parents c14326c + 47730c6
commit 31c02ed
Show file tree

Hide file tree

Showing 23 changed files with 291,596 additions and 708 deletions.
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -1,6 +1,12 @@
 Change Log
 ==========
 
+v3.2.2
+----------
+- Consolidate ``SnB``/``doped`` ``INCAR`` defaults and remove redundant settings.
+- Ensure backwards compatiblity in defect folder name handling.
+- Fix bug in ``get_site_magnetizations``.
+
 v3.2.1
 ----------
 - Update CLI config handling.

diff --git a/README.md b/README.md
@@ -46,7 +46,7 @@ If using `VASP`, in order for `ShakeNBreak` to automatically generate the pseudo
   PMG_VASP_PSP_DIR: <Path to VASP pseudopotential top directory>
 ```
    Within your `VASP` pseudopotential top directory, you should have a folder named `POT_GGA_PAW_PBE`
-   which contains the `POTCAR.X(.gz)` files (in this case for PBE `POTCAR`s). Please refer to the `doped Installation docs <https://doped.readthedocs.io/en/latest/Installation.html>`_ if you have
+   which contains the `POTCAR.X(.gz)` files (in this case for PBE `POTCAR`s). Please refer to the [`doped` Installation docs](https://doped.readthedocs.io/en/latest/Installation.html) if you have
    difficulty with this.
 
 The font Montserrat ([Open Font License](https://scripts.sil.org/cms/scripts/page.php?site_id=nrsi&id=OFL)) will be installed with the package, and will be used by default for plotting.

diff --git a/SnB_input_files/incar.yaml b/SnB_input_files/incar.yaml
@@ -1,27 +1,15 @@
-'# May want to change NCORE, KPAR, AEXX, ENCUT, NUPDOWN, ISPIN, POTIM': ''
-'# ShakeNBreak INCAR with coarse settings to maximise speed with sufficient accuracy for qualitative structure searching': ''
-'# KPAR': '# No KPAR, only one kpoint'
+# ShakeNBreak uses the default doped INCAR settings described here: https://doped.readthedocs.io/en/latest/doped.VASP_sets.html
+# as well as these modifications:
+# Changed from doped default:
+LCHARG: false  # reduce file sizes
+LREAL: Auto  # real-space force projection, runtime speedup with acceptable loss in force accuracy
+LWAVE: false  # reduce file sizes
+NELM: 40  # slightly reduced NELM to speedup with incomplete SCF convergence far from minimum
+NSW: 300  # increased NSW from doped default to allow relaxation from initially distorted structure
+
+# Same as doped but with added comment information in INCAR file:
 ALGO: "Normal  # Change to All if ZHEGV, FEXCP/F or ZBRENT errors encountered (done automatically by snb-run)"
-EDIFFG: -0.01
-ENCUT: 300
-HFSCREEN: 0.208  # correct HSE screening parameter; see https://aip.scitation.org/doi/10.1063/1.2404663
-# Note this HFSCREEN value differs from the Materials Project MPHSERelaxSet default of 0.2! This should
-# be consistent between all your defect/bulk/competing-phase calculations.
 IBRION: '2  # Typically more stable / reliable than "1" (RMM-DIIS), but change if ionic convergence is poor (done automatically by snb-run)'
-ISIF: 2
-ISMEAR: 0
-ISPIN: '2  # Spin polarisation likely for defects'
-ISYM: '0  # Symmetry breaking extremely likely for defects'
-LASPH: true
-LCHARG: false
-LHFCALC: true
-LORBIT: 11  # lm-decomposed orbital projections; useful for DOS analysis (e.g. sumo-dosplot)
-LREAL: Auto
-LWAVE: false
-NCORE: 16
-NEDOS: 2000  # dense energy mesh output in OUTCAR; useful for DOS analysis (e.g. sumo-dosplot)
-NELM: 40
-NSW: 300
-PREC: Accurate
-PRECFOCK: Fast
-SIGMA: 0.05
+
+# Comments:
+'# ShakeNBreak INCAR with coarse settings to maximise speed with sufficient accuracy for qualitative structure searching': ''
diff --git a/docs/conf.py b/docs/conf.py
@@ -25,7 +25,7 @@
 author = 'Irea Mosquera-Lois, Seán R. Kavanagh'
 
 # The full version, including alpha/beta/rc tags
-release = '3.2.1'
+release = '3.2.2'
 
 
 # -- General configuration ---------------------------------------------------

diff --git a/setup.py b/setup.py
@@ -130,7 +130,7 @@ def package_files(directory):
 
 setup(
     name="shakenbreak",
-    version="3.2.1",
+    version="3.2.2",
     description="Package to generate and analyse distorted defect structures, in order to "
     "identify ground-state and metastable defect configurations.",
     long_description=long_description,

diff --git a/shakenbreak/analysis.py b/shakenbreak/analysis.py
@@ -108,7 +108,7 @@ def _get_distortion_filename(distortion) -> str:
         distortion (:obj:`str`):
             distortion label used for file names.
     """
-    if isinstance(distortion, float) or isinstance(distortion, int):
+    if isinstance(distortion, (float, int)):
         if distortion != 0:
             distortion_label = f"Bond_Distortion_{round(distortion * 100, 1)+0}%"
             # as percentage with 1 decimal place (e.g. 50.0%)
@@ -118,10 +118,11 @@ def _get_distortion_filename(distortion) -> str:
         if "_from_" in distortion and "Rattled" not in distortion:
             distortion_label = f"Bond_Distortion_{distortion}"
             # runs from other charge states
-        elif "Rattled_from_" in distortion:
+        elif "Rattled_from_" in distortion or distortion in [
+            "Unperturbed",
+            "Rattled",
+        ]:
             distortion_label = distortion
-        elif distortion == "Unperturbed" or distortion == "Rattled":
-            distortion_label = distortion  # e.g. "Unperturbed"/"Rattled"
         else:
             try:  # try converting to float, in case user entered '0.5'
                 distortion = float(distortion)
@@ -281,15 +282,13 @@ def _sort_data(
     if verbose:
         if energy_diff and float(energy_diff) < -min_e_diff:
             print(
-                f"{defect_name}: Energy difference between minimum, found with "
-                f"{gs_distortion} bond distortion, and unperturbed: "
-                f"{energy_diff:+.2f} eV."
+                f"{defect_name}: Energy difference between minimum, found with {gs_distortion} bond "
+                f"distortion, and unperturbed: {energy_diff:+.2f} eV."
             )
         elif energy_diff is None:
             print(
-                f"{defect_name}: Unperturbed energy not found in {energies_file}. "
-                f"Lowest energy structure found with {gs_distortion} bond "
-                f"distortion."
+                f"{defect_name}: Unperturbed energy not found in {energies_file}. Lowest energy "
+                f"structure found with {gs_distortion} bond distortion."
             )
     return defect_energies_dict, energy_diff, gs_distortion
 
@@ -1122,15 +1121,18 @@ def get_site_magnetizations(
     if not os.path.exists(f"{output_path}/{defect_species}"):
         raise FileNotFoundError(f"{output_path}/{defect_species} does not exist!")
 
-    if not defect_site:  # look for defect site, in order to include the distance
+    defect_site_coords = None
+    if isinstance(defect_site, list) or isinstance(defect_site, np.ndarray):
+        defect_site_coords = defect_site
+    elif not defect_site:  # look for defect site, in order to include the distance
         # between sites with significant magnetization and the defect
         if os.path.exists(f"{output_path}/distortion_metadata.json"):
             with open(f"{output_path}/distortion_metadata.json", "r") as f:
                 try:
                     defect_species_without_charge = "_".join(
                         defect_species.split("_")[0:-1]
                     )  # remove charge state
-                    defect_site = json.load(f)["defects"][
+                    defect_site_coords = json.load(f)["defects"][
                         defect_species_without_charge
                     ]["unique_site"]
                 except KeyError:
@@ -1162,10 +1164,12 @@ def get_site_magnetizations(
                 "found. Skipping magnetisation analysis."
             )
             continue
-        if isinstance(defect_site, list) or isinstance(defect_site, np.ndarray):
+        if isinstance(defect_site_coords, list) or isinstance(
+            defect_site_coords, np.ndarray
+        ):
             # for vacancies, append fake atom
             structure.append(
-                species="V", coords=defect_site, coords_are_cartesian=False
+                species="V", coords=defect_site_coords, coords_are_cartesian=False
             )
             defect_site = -1  # index of the added fake atom
         if not os.path.exists(f"{output_path}/{defect_species}/{dist_label}/OUTCAR"):

diff --git a/shakenbreak/cli.py b/shakenbreak/cli.py
@@ -863,10 +863,10 @@ def parse(defect, all, path, code):
     Parsed energies are written to a `yaml` file in the corresponding defect directory.
     """
     if defect:
-        io.parse_energies(defect, path, code)
+        _ = io.parse_energies(defect, path, code)
     elif all:
         defect_dirs = _parse_defect_dirs(path)
-        [io.parse_energies(defect, path, code) for defect in defect_dirs]
+        _ = [io.parse_energies(defect, path, code) for defect in defect_dirs]
     else:
         # assume current directory is the defect folder
         try:
@@ -878,7 +878,7 @@ def parse(defect, all, path, code):
             cwd = os.getcwd()
             defect = cwd.split("/")[-1]
             path = cwd.rsplit("/", 1)[0]
-            io.parse_energies(defect, path, code)
+            _ = io.parse_energies(defect, path, code)
         except Exception:
             raise Exception(
                 f"Could not parse defect '{defect}' in directory '{path}'. Please either specify "
@@ -952,8 +952,15 @@ def analyse(defect, all, path, code, ref_struct, verbose):
 
     def analyse_single_defect(defect, path, code, ref_struct, verbose):
         if not os.path.exists(f"{path}/{defect}") or not os.path.exists(path):
-            raise FileNotFoundError(f"Could not find {defect} in the directory {path}.")
-        io.parse_energies(defect, path, code)
+            orig_defect_name = defect
+            defect = defect.replace("+", "")  # try with old name format
+
+            if not os.path.exists(f"{path}/{defect}") or not os.path.exists(path):
+                raise FileNotFoundError(
+                    f"Could not find {orig_defect_name} in the directory {path}."
+                )
+
+        _ = io.parse_energies(defect, path, code)
         defect_energies_dict = analysis.get_energies(
             defect_species=defect, output_path=path, verbose=verbose
         )
@@ -973,12 +980,13 @@ def analyse_single_defect(defect, path, code, ref_struct, verbose):
         for defect in defect_dirs:
             print(f"\nAnalysing {defect}...")
             analyse_single_defect(defect, path, code, ref_struct, verbose)
+
     elif defect is None:
         # assume current directory is the defect folder
         if path != ".":
             warnings.warn(
-                "`--path` option ignored when running from within defect folder ("
-                "i.e. when `--defect` is not specified."
+                "`--path` option ignored when running from within defect folder (i.e. when `--defect` is "
+                "not specified."
             )
         cwd = os.getcwd()
         defect = cwd.split("/")[-1]
@@ -1002,10 +1010,10 @@ def analyse_single_defect(defect, path, code, ref_struct, verbose):
             analyse_single_defect(defect, orig_path, code, ref_struct, verbose)
         except Exception:
             raise Exception(
-                f"Could not analyse defect '{defect}' in directory '{path}'. Please either "
-                f"specify a defect to analyse (with option --defect), run from within a single "
-                f"defect directory (without setting --defect) or use the --all flag to analyse "
-                f"all defects in the specified/current directory."
+                f"Could not analyse defect '{defect}' in directory '{path}'. Please either specify a "
+                f"defect to analyse (with option --defect), run from within a single defect directory ("
+                f"without setting --defect) or use the --all flag to analyse all defects in the "
+                f"specified/current directory."
             )
 
 
@@ -1148,7 +1156,7 @@ def plot(
         for defect in defect_dirs:
             if verbose:
                 print(f"Parsing {defect}...")
-            io.parse_energies(defect, path, code)
+            _ = io.parse_energies(defect, path, code)
         # Create defects_dict (matching defect name to charge states)
         defects_wout_charge = [defect.rsplit("_", 1)[0] for defect in defect_dirs]
         defects_dict = {
@@ -1192,14 +1200,17 @@ def plot(
     else:
         orig_path = None
     try:
-        io.parse_energies(defect, path, code)
+        energies_file = io.parse_energies(defect, path, code)
+        defect_species = energies_file.rsplit("/", 1)[-1].replace(
+            ".yaml", ""
+        )  # in case '+' removed
         defect_energies_dict = analysis.get_energies(
-            defect_species=defect,
+            defect_species=defect_species,
             output_path=path,
             verbose=verbose,
         )
         plotting.plot_defect(
-            defect_species=defect,
+            defect_species=defect_species,
             energies_dict=defect_energies_dict,
             output_path=path,
             add_colorbar=colorbar,
@@ -1212,14 +1223,17 @@ def plot(
         )
     except Exception:
         try:
-            io.parse_energies(defect, orig_path, code)
+            energies_file = io.parse_energies(defect, orig_path, code)
+            defect_species = energies_file.rsplit("/", 1)[-1].replace(
+                ".yaml", ""
+            )  # in case '+' removed
             defect_energies_dict = analysis.get_energies(
-                defect_species=defect,
+                defect_species=defect_species,
                 output_path=orig_path,
                 verbose=verbose,
             )
             plotting.plot_defect(
-                defect_species=defect,
+                defect_species=defect_species,
                 energies_dict=defect_energies_dict,
                 output_path=orig_path,
                 add_colorbar=colorbar,
@@ -1392,31 +1406,24 @@ def groundstate(
     """
     # determine if running from within a defect directory or from the top level directory
     if any(
-        [
-            dir
-            for dir in os.listdir()
-            if os.path.isdir(dir)
-            and any(
-                [
-                    substring in dir
-                    for substring in ["Bond_Distortion", "Rattled", "Unperturbed"]
-                ]
-            )
-        ]
+        dir
+        for dir in os.listdir()
+        if os.path.isdir(dir)
+        and any(
+            substring in dir
+            for substring in ["Bond_Distortion", "Rattled", "Unperturbed"]
+        )
     ):  # distortion subfolders in cwd
         # check if defect folders also in cwd
         for dir in [dir for dir in os.listdir() if os.path.isdir(dir)]:
             defect_name = None
             try:
                 defect_name = _format_defect_name(dir, include_site_info_in_name=False)
             except Exception:
-                try:
+                with contextlib.suppress(Exception):
                     defect_name = _format_defect_name(
                         f"{dir}_0", include_site_info_in_name=False
                     )
-                except Exception:
-                    pass
-
             if (
                 defect_name
             ):  # recognised defect folder found in cwd, warn user and proceed