Merge pull request #56 from SMTG-UCL/develop

`v3.1.0`; extending compatibility to `python=3.11`

.github/dependabot.yml

@@ -8,3 +8,7 @@ updates:
     python-packages:
       patterns:
         - "*"
+  open-pull-requests-limit: 2
+  allow:
+    - dependency-type: direct
+    - dependency-type: indirect

.github/workflows/build_and_test.yml

@@ -18,8 +18,7 @@ jobs:
       max-parallel: 5
 
       matrix:
-        python-version: ['3.8', '3.9', '3.10']
-        # hiphive requires numba -> currently incompatible with Python 3.11
+        python-version: ['3.8', '3.9', '3.10', '3.11']
 
     name: Python ${{ matrix.python-version }} Test Pop
 

.github/workflows/pip_install_test.yml

@@ -15,7 +15,7 @@ jobs:
       max-parallel: 5
 
       matrix:
-        python-version: ['3.8', '3.9', '3.10']
+        python-version: ['3.8', '3.9', '3.10', '3.11']
 
     name: Python ${{ matrix.python-version }} pip install
 

.github/workflows/release.yml

@@ -57,7 +57,7 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
         with:
-          tag_name: ${{ env.version }}
+          tag_name: v${{ env.version }}
           release_name: ${{ env.version }}
           body_path: release_info.txt
           draft: false

CHANGELOG.rst

@@ -1,6 +1,11 @@
 Change Log
 ==========
 
+v3.1.0
+----------
+- Update dependencies, as `hiphive=1.2` has been released, making `ShakeNBreak` compatible with
+  `python=3.11`🎉
+
 v3.0.0
 ----------
 - Switch to semantic versioning

README.md

@@ -126,6 +126,7 @@ Automatic testing is run on the master and develop branches using Github Actions
 
 ## Studies using `ShakeNBreak`
 We'll add papers that use `ShakeNBreak` to this list as they come out!
+- J. Willis, K. B. Spooner, D. O. Scanlon. [_ChemRxiv_](https://chemrxiv.org/engage/chemrxiv/article-details/64c29140ce23211b20a787bb) 2023
 - A. T. J. Nicolson et al. [_Journal of Materials Chemistry A_](https://doi.org/10.1039/D3TA02429F) 2023
 - X. Wang et al. [_arXiv_](https://arxiv.org/abs/2302.04901) 2023
 - J. Cen et al. [_Journal of Materials Chemistry A_](https://doi.org/10.1039/D3TA00532A) 2023
@@ -160,7 +161,7 @@ You may also find this Preview paper useful, which discusses the general problem
 `BibTeX` entries for these papers are provided in the [`CITATIONS.md`](CITATIONS.md) file.
 
 ## Requirements
-`ShakeNBreak` is compatible with Python 3.8, 3.9 & 3.10 and requires the following open-source python packages:
+`ShakeNBreak` is compatible with Python 3.8 - 3.11 and requires the following open-source python packages:
 * [Pymatgen](https://pymatgen.org/)
 * [Ase](https://wiki.fysik.dtu.dk/ase/)
 * [Hiphive](https://hiphive.materialsmodeling.org/)

docs/Analysis.rst

@@ -265,4 +265,4 @@ In particular, symmetry-breaking as a result of structural reconstruction from t
 the defect, which should be accounted for when later computing concentrations and Fermi level position.
 These considerations, as well as the importance of metastability and temperature effects for the free
 energies (and thus concentrations) for certain defects/systems are discussed in this Tutorial Review
-paper: `Imperfections are not 0 K: free energy of point defects in crystals, 10.48550/arXiv.2307.10451 <https://doi.org/10.48550/arXiv.2307.10451>`_.
+paper: `Imperfections are not 0 K: free energy of point defects in crystals (Chem Soc Rev 2023) <https://doi.org/10.1039/D3CS00432E>`_.

docs/index.rst

@@ -213,6 +213,7 @@ Studies using ``ShakeNBreak``
 
 We'll add papers that use `ShakeNBreak` to this list as they come out!
 
+- J\. Willis, K. B. Spooner, D. O. Scanlon. `ChemRxiv <https://chemrxiv.org/engage/chemrxiv/article-details/64c29140ce23211b20a787bb>`__ 2023
 - A\. T. J. Nicolson et al. `Journal of Materials Chemistry A <https://doi.org/10.1039/D3TA02429F>`__ 2023
 - X\. Wang et al. `arXiv`_ 2023
 - J\. Cen et al. `Journal of Materials Chemistry A`_ 2023
@@ -230,7 +231,6 @@ We'll add papers that use `ShakeNBreak` to this list as they come out!
 .. Kumagai collab paper
 .. Lavan LiNiO2
 .. Sykes Magnetic oxide polarons
-.. Joe CuI preprint out soon
 .. Kat YTOS
 .. Squires (and mention benchmark test against AIRSS? See Slack message)
 
@@ -268,7 +268,7 @@ You may also find this Preview paper useful, which discusses the general problem
 Requirements
 ========================
 
-``ShakeNBreak`` is compatible with Python 3.8, 3.9 & 3.10 and requires the following open-source python packages:
+``ShakeNBreak`` is compatible with Python 3.8 - 3.11 and requires the following open-source python packages:
 
 * `Pymatgen <https://pymatgen.org/>`_
 * `Ase <https://wiki.fysik.dtu.dk/ase/>`_

setup.py

@@ -131,7 +131,7 @@ def package_files(directory):
 
 setup(
     name="shakenbreak",
-    version="3.0.0",
+    version="3.1.0",
     description="Package to generate and analyse distorted defect structures, in order to "
     "identify ground-state and metastable defect configurations.",
     long_description=long_description,
@@ -153,21 +153,22 @@ def package_files(directory):
         "Programming Language :: Python :: 3.8",
         "Programming Language :: Python :: 3.9",
         "Programming Language :: Python :: 3.10",
+        "Programming Language :: Python :: 3.11",
         "Topic :: Scientific/Engineering :: Chemistry",
         "Topic :: Scientific/Engineering :: Physics",
     ],
     keywords="chemistry pymatgen dft defects structure-searching distortions symmetry-breaking",
     packages=find_packages(),
-    python_requires=">=3.8, <3.11",
+    python_requires=">=3.8",
     install_requires=[
-        "numpy<1.24,>=1.21.2",
+        "numpy", #>=1.21.2" needed for numpy.typing.NDArray?
         "pymatgen>=2022.10.22",
         "pymatgen-analysis-defects>=2022.10.28",
         "matplotlib",
         "ase",
         "pandas>=1.1.0",
         "seaborn",
-        "hiphive",
+        "hiphive>=1.0",  # nbr_cutoff not defined in previous versions of mc_rattle
         "monty",
         "click>8.0",
         "importlib_metadata",

shakenbreak/cli.py

@@ -1511,19 +1511,19 @@ def mag(outcar, threshold, verbose):
         outcar_obj = Outcar(outcar)
         abs_mag_values = [abs(m["tot"]) for m in outcar_obj.magnetization]
 
+        if (
+            max(abs_mag_values) < threshold  # no one atomic moment greater than threshold
+            and sum(abs_mag_values) < threshold * 10  # total moment less than 10x threshold
+        ):
+            if verbose:
+                print(f"Magnetisation is below threshold (<{threshold} μB/atom)")
+            sys.exit(0)
+        else:
+            if verbose:
+                print(f"Magnetisation is above threshold (>{threshold} μB/atom)")
+            sys.exit(1)
+
     except Exception:
         if verbose:
             print(f"Could not read magnetisation from OUTCAR file at {outcar}")
         sys.exit(1)
-
-    if (
-        max(abs_mag_values) < threshold  # no one atomic moment greater than threshold
-        and sum(abs_mag_values) < threshold * 10  # total moment less than 10x threshold
-    ):
-        if verbose:
-            print(f"Magnetisation is below threshold (<{threshold} μB/atom)")
-        sys.exit(0)
-    else:
-        if verbose:
-            print(f"Magnetisation is above threshold (>{threshold} μB/atom)")
-        sys.exit(1)

shakenbreak/distortions.py

@@ -270,7 +270,7 @@ def rattle(
     try:
         rattled_ase_struct = generate_mc_rattled_structures(
             ase_struct,
-            n_configs=1,
+            1,  # n_configs in hiphive <= 1.1, n_structures in hiphive >= 1.2
             rattle_std=stdev,
             d_min=d_min,
             n_iter=n_iter,
@@ -289,7 +289,7 @@ def rattle(
                 try:
                     rattled_ase_struct = generate_mc_rattled_structures(
                         ase_struct,
-                        n_configs=1,
+                        1,  # n_configs in hiphive <= 1.1, n_structures in hiphive >= 1.2
                         rattle_std=stdev,
                         d_min=reduced_d_min,
                         n_iter=n_iter,

shakenbreak/io.py

@@ -112,7 +112,9 @@ def parse_vasp_energy(defect_dir, dist, energy, outcar):
             outcar = os.path.join(defect_dir, dist, "OUTCAR")
         if outcar:  # regrep faster than using Outcar/vasprun class
             try:
-                energy = _match(outcar, r"energy\(sigma->0\)\s+=\s+([\d\-\.]+)")[0][0][
+                energy = _match(
+                    outcar, r"entropy=.*energy\(sigma->0\)\s+=\s+([\d\-\.]+)"
+                )[0][0][
                     0
                 ]  # Energy of first match
                 converged = _match(

shakenbreak/plotting.py

@@ -218,7 +218,6 @@ def _format_defect_name(
     recognised_post_interstitial_strings = sorted(
         [
             "_i",  # but not '_I' as could be iodine
-            "i",  # but not 'I' as could be iodine
             "_int",
             "_Int",
             "int",
@@ -496,12 +495,22 @@ def _defect_name_from_matching_elements(
                 if dummy_h.is_valid_symbol(character)
             ]
             if len(possible_one_character_elements) > 0:
+                # in this case, we don't know the order of the 1-character vs 2-character elements in
+                # the name, so we try both orderings:
                 defect_name = _defect_name_from_matching_elements(
-                    possible_one_character_elements + possible_two_character_elements,
+                    possible_two_character_elements + possible_one_character_elements,
                     pre_charge_name,  # trimmed_pre_charge_name name for finding elements,
                     # pre_charge_name for matching defect format
                     include_site_num_in_name,
                 )
+                if defect_name is None:
+                    defect_name = _defect_name_from_matching_elements(
+                        possible_one_character_elements
+                        + possible_two_character_elements,
+                        pre_charge_name,  # trimmed_pre_charge_name name for finding elements,
+                        # pre_charge_name for matching defect format
+                        include_site_num_in_name,
+                    )
 
     if defect_name is None:
         # try single-character element match
@@ -713,7 +722,9 @@ def _remove_high_energy_points(
             > max_energy_above_unperturbed
         ):
             energies_dict["distortions"].pop(key)
-            if disp_dict:  # only exists if user selected `add_colorbar=True`
+            if (
+                disp_dict and key in disp_dict
+            ):  # only exists if user selected `add_colorbar=True`
                 disp_dict.pop(key)
     return energies_dict, disp_dict
 
@@ -1335,7 +1346,7 @@ def plot_all_defects(
                         FileNotFoundError
                     ):  # distortion_metadata.json not found in this folder
                         pass
-                if distortion_metadata and type(distortion_metadata) == dict:
+                if distortion_metadata and isinstance(distortion_metadata, dict):
                     num_nearest_neighbours, neighbour_atom = _parse_distortion_metadata(
                         distortion_metadata, defect, charge
                     )
@@ -1786,7 +1797,7 @@ def plot_colorbar(
             # Datapoints from other charge states
             if imported_indices:
                 other_charges = len(
-                    set(
+                    list(
                         list(energies_dict["distortions"].keys())[i].split("_")[-1]
                         for i in imported_indices.keys()
                     )
@@ -1803,7 +1814,11 @@ def plot_colorbar(
                         edgecolors="k",
                         ls="-",
                         s=50,
-                        marker=["s", "v", "<", ">", "^", "p", "X"][j],
+                        marker=(
+                            ["s", "v", "<", ">", "^", "p", "X"] * 3
+                        )[  # repeat markers in case many imported charge states
+                            j
+                        ],  # different markers for different charge states
                         zorder=10,  # make sure it's on top of the other points
                         cmap=colormap,
                         norm=norm,
@@ -1891,8 +1906,8 @@ def plot_colorbar(
 
 def plot_datasets(
     datasets: list,
-    dataset_labels: list,
-    defect_species: str,
+    dataset_labels: Optional[list] = None,
+    defect_species: str = "defect",
     include_site_num_in_name: Optional[bool] = False,
     title: Optional[str] = None,
     neighbour_atom: Optional[str] = None,
@@ -1918,10 +1933,11 @@ def plot_datasets(
             matching distortion to final energy (eV), as produced by
             `analysis._organize_data()` or `analysis.get_energies()`)
         dataset_labels (:obj:`list`):
-            Labels for each dataset plot legend.
+            Labels for each dataset plot legend. If None, defaults to
+            ["Distortions"]*len(datasets).
         defect_species (:obj:`str`):
             Specific defect name that will appear in plot labels (in LaTeX form)
-             and file names (e.g 'vac_1_Cd_0')
+             and file names (e.g 'vac_1_Cd_0'). Defaults to 'defect'.
         include_site_num_in_name (:obj:`bool`):
             Whether to include the site number (as generated by doped) in the
             defect name. Useful for materials with many symmetry-inequivalent
@@ -1976,12 +1992,15 @@ def plot_datasets(
             as a mpl.figure.Figure object
     """
     # Validate input
-    if len(datasets) != len(dataset_labels):
-        raise ValueError(
-            f"Number of datasets and labels must match! "
-            f"You gave me {len(datasets)} datasets and"
-            f" {len(dataset_labels)} labels."
-        )
+    if dataset_labels is not None:
+        if len(datasets) != len(dataset_labels):
+            raise ValueError(
+                f"Number of datasets and labels must match! "
+                f"You gave me {len(datasets)} datasets and"
+                f" {len(dataset_labels)} labels."
+            )
+    else:
+        dataset_labels = ["Distortions"] * len(datasets)
 
     _install_custom_font()
     # Set up figure
@@ -2111,7 +2130,7 @@ def plot_datasets(
 
             if imported_indices:
                 other_charges = len(
-                    set(
+                    list(
                         list(dataset["distortions"].keys())[i].split("_")[-1]
                         for i in imported_indices
                     )
@@ -2128,7 +2147,9 @@ def plot_datasets(
                         ls="-",
                         s=50,
                         zorder=10,  # make sure it's on top of the other lines
-                        marker=["s", "v", "<", ">", "^", "p", "X"][
+                        marker=(
+                            ["s", "v", "<", ">", "^", "p", "X"] * 3
+                        )[  # repeat markers in case many imported charge states
                             j
                         ],  # different markers for different charge states
                         alpha=1,

shakenbreak/scripts/SnB_run.sh

@@ -121,10 +121,13 @@ SnB_run_loop() {
         # if electronic convergence is not being reached, change ALGO to All
         if grep -q "aborting loop EDIFF was not reached" OUTCAR && ! grep -q "aborting loop because EDIFF is reached" OUTCAR; then
           # unconverged electronic loops and no converged electronic loops
-          if [ "$verbose" = true ]; then
-            echo "${i%?} is showing poor electronic convergence, changing ALGO to All."
+          ialgo=$(grep IALGO  OUTCAR | awk '{print $3}')
+          if [ "$ialgo" != "58" ]; then
+            if [ "$verbose" = true ]; then
+              echo "${i%?} is showing poor electronic convergence, changing ALGO to All."
+            fi
+            sed -i 's/ALGO.*/ALGO = All/g' INCAR
           fi
-          sed -i 's/ALGO.*/ALGO = All/g' INCAR
         fi
 
         # check if multiple <=single-step OUTCARs present, and CONTCAR empty/less than 9 lines or same as POSCAR

tests/test_cli.py

@@ -2098,6 +2098,27 @@ def _test_OUTCAR_error(error_string):
         self.assertEqual(len([i for i in saved_files if "INCAR" in i]), 1)
         self.assertEqual(len([i for i in saved_files if "CONTCAR" in i]), 1)
         self.assertEqual(len([i for i in saved_files if "OUTCAR" in i]), 1)
+        for i in saved_files:
+            os.remove(f"Bond_Distortion_10.0%/{i}")
+        if_present_rm("Bond_Distortion_10.0%/job_file")
+
+        # test changing no message with poor electronic convergence when ALGO already = All
+        with open("Bond_Distortion_10.0%/OUTCAR", "w") as fp:
+            fp.write(poor_electronic_convergence_outcar_string + "\nIALGO = 58  # i.e. ALGO = All")
+
+        proc = subprocess.Popen(
+            ["snb-run", "-v", "-s echo", "-n this", "-j job_file"],
+            stdout=subprocess.PIPE,
+            stderr=subprocess.PIPE,
+        )  # setting 'job command' to 'echo' to
+        out = str(proc.communicate()[0])
+        self.assertNotIn(
+            "Bond_Distortion_10.0% is showing poor electronic convergence, changing ALGO to All.",
+            out,
+        )
+        self.assertIn("ALGO = All", open("Bond_Distortion_10.0%/INCAR").read())
+        files = os.listdir("Bond_Distortion_10.0%")
+        saved_files = [file for file in files if "on" in file and "CAR_" in file]
         for i in saved_files:
             os.remove(f"Bond_Distortion_10.0%/{i}")
         os.remove("Bond_Distortion_10.0%/OUTCAR")
@@ -2521,7 +2542,7 @@ def test_parse(self):
         )  # Bond_Distortion_-20.0%_not_converged now included
         not_converged_energies = copy.deepcopy(test_energies)
         not_converged_energies["distortions"].update(
-            {"Bond_Distortion_-20.0%_not_converged": -1151.8383839}
+            {"Bond_Distortion_-20.0%_not_converged": -1110.37833497}
         )
         self.assertEqual(
             not_converged_energies, energies