Merge pull request #33 from SMTG-UCL/develop

Develop Former-commit-id: b83105a
SMTG-Bham · Sep 4, 2023 · 312b026 · 312b026
2 parents 246960f + 83dd63d
commit 312b026
Show file tree

Hide file tree

Showing 56 changed files with 1,425 additions and 452 deletions.
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -1,6 +1,15 @@
 Change Log
 ==========
 
+v.2.0.4
+----------
+- Add supercell re-ordering tests for parsing
+- Ensure final _relaxed_ defect site (for interstitials and substitutions) is used for finite-size
+  charge corrections
+- Consolidate functions and input sets with `ShakeNBreak`
+- Update defect generation tests
+- Use more efficient Wyckoff determination code
+
 v.2.0.3
 ----------
 - Sort defect entries in ``DefectPhaseDiagram`` for deterministic behaviour (particularly for plotting).

diff --git a/README.md b/README.md
@@ -63,19 +63,17 @@ and morphed over time as more and more functionality was added. After breaking c
 entirely refactored and rewritten, to work with the new
 `pymatgen-analysis-defects` package.
 
-## Studies using `doped`
-
-We'll add papers that use `doped` to this list as they come out!
-
-- A. T. J. Nicolson et al. [Journal of Materials Chemistry A](https://doi.org/10.1039/D3TA02429F) 2023
-- Y. W. Woo, Z. Li, Y-K. Jung, J-S. Park, A. Walsh [ACS Energy Letters](https://doi.org/10.1021/acsenergylett.2c02306) 2023
-- P. A. Hyde et al. [Inorganic Chemistry](https://doi.org/10.1021/acs.inorgchem.3c01510) 2023
-- J. Willis, K. B. Spooner, D. O. Scanlon. [ChemRxiv](https://chemrxiv.org/engage/chemrxiv/article-details/64c29140ce23211b20a787bb) 2023
-- X. Wang et al. arXiv 2023
-- J. Cen et al. [Journal of Materials Chemistry A](https://doi.org/10.1039/D3TA00532A) 2023
-- J. Willis & R. Claes et al. [ChemRxiv](https://doi.org/10.26434/chemrxiv-2023-lttnf) 2023
-- I. Mosquera-Lois & S. R. Kavanagh, A. Walsh, D. O. Scanlon [npj Computational Materials](https://www.nature.com/articles/s41524-023-00973-1) 2023
-- Y. T. Huang & S. R. Kavanagh et al. [Nature Communications](https://www.nature.com/articles/s41467-022-32669-3) 2022
-- S. R. Kavanagh, D. O. Scanlon, A. Walsh, C. Freysoldt [Faraday Discussions](https://doi.org/10.1039/D2FD00043A) 2022
-- S. R. Kavanagh, D. O. Scanlon, A. Walsh [ACS Energy Letters](https://pubs.acs.org/doi/full/10.1021/acsenergylett.1c00380) 2021
-- C. J. Krajewska et al. [Chemical Science](https://doi.org/10.1039/D1SC03775G) 2021
+## Studies using `doped` (so far)
+
+- A. T. J. Nicolson et al. [_Journal of Materials Chemistry A_](https://doi.org/10.1039/D3TA02429F) 2023
+- Y. W. Woo, Z. Li, Y-K. Jung, J-S. Park, A. Walsh [_ACS Energy Letters_](https://doi.org/10.1021/acsenergylett.2c02306) 2023
+- P. A. Hyde et al. [_Inorganic Chemistry_](https://doi.org/10.1021/acs.inorgchem.3c01510) 2023
+- J. Willis, K. B. Spooner, D. O. Scanlon. [_ChemRxiv_](https://chemrxiv.org/engage/chemrxiv/article-details/64c29140ce23211b20a787bb) 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
+- J. Willis & R. Claes et al. [_ChemRxiv_](https://doi.org/10.26434/chemrxiv-2023-lttnf) 2023
+- I. Mosquera-Lois & S. R. Kavanagh, A. Walsh, D. O. Scanlon [_npj Computational Materials_](https://www.nature.com/articles/s41524-023-00973-1) 2023
+- Y. T. Huang & S. R. Kavanagh et al. [_Nature Communications_](https://www.nature.com/articles/s41467-022-32669-3) 2022
+- S. R. Kavanagh, D. O. Scanlon, A. Walsh, C. Freysoldt [_Faraday Discussions_](https://doi.org/10.1039/D2FD00043A) 2022
+- S. R. Kavanagh, D. O. Scanlon, A. Walsh [_ACS Energy Letters_](https://pubs.acs.org/doi/full/10.1021/acsenergylett.1c00380) 2021
+- C. J. Krajewska et al. [_Chemical Science_](https://doi.org/10.1039/D1SC03775G) 2021
diff --git a/ToDo.md b/ToDo.md
@@ -101,6 +101,7 @@
   - Note that if you edit the entries in a DefectPhaseDiagram after creating it, you need to `dpd.find_stable_charges()` to update the transition level map etc.
 - Should tag parsed defects with `is_shallow` (or similar), and then omit these from plotting/analysis
   (and note this behaviour in examples/docs)
+- Ideally our defect parsing would be able to get the final _relaxed_ position of vacancies / antisites that move significantly (or the centroid if a defect cluster), to then use for the charge correction. Not a big deal for larger supercells, but a slight mismatch in defect site prediction for smaller supercells can have a semi-significant effect on the predicted charge correction. `Int_Te_3_unperturbed_1` is a good example of this tricky case.
 - Change formation energy plotting and tabulation to DefectPhaseDiagram methods rather than standalone
   functions – with `pymatgen` update what's the new architecture?
 - Better automatic defect formation energy plot colour handling (auto-change colormap based on number of defects, set similar colours for similar defects (types and inequivalent sites)) – and more customisable?
@@ -143,12 +144,7 @@
 ## Housekeeping
 - Clean `README` with bullet-point summary of key features, and sidebar like `SnB`.
 - `ShakeNBreak` related updates:
-  - Update to use the `ShakeNBreak` voronoi node-finding functions, as this has been made to be more
-    efficient than the `doped` version (which is already far more efficient than the original...) and
-    isn't available in current `pymatgen`.
-  - Use doped naming conventions and functions, site-matching/symmetry functions, defect entry generation
-    functions, `DefectSet` (and anything else?) in `ShakeNBreak`. Streamline SnB notebook with these!!
-  - Should have quick scan through both codes to ensure no redundant/duplicate functions.
+  - Use doped naming conventions and functions and defect entry generation functions in `ShakeNBreak`.
 - Code tidy up:
   - Notebooks in `tests`; update or delete.
   - Test coverage?
@@ -183,6 +179,7 @@
     and `z`, and checking that everything is zero (not net magnetisation, as could have opposing spin
     bipolaron). This is automatically handled in `SnB_replace_mag.py` (to be added to ShakeNBreak) and
     will be added to `doped` VASP calc scripts.
+  - Setting `LREAL = Auto` can sometimes be worth doing if you have a very large supercell for speed up, _but_ it's important to do a final calculation with `LREAL = False` for accurate energies/forces, so only do if you're a power user and have a very large supercell.
   - Show usage of `get_conv_cell_site` in notebooks/docs.
   - Note in docs that `spglib` convention used for Wyckoff labels and conventional structure definition.
     Primitive structure can change, as can supercell / supercell matrix (depending on input structure,

diff --git a/docs/conf.py b/docs/conf.py
@@ -25,7 +25,7 @@
 author = 'Seán R. Kavanagh'
 
 # The full version, including alpha/beta/rc tags
-release = '2.0.3'
+release = '2.0.4'
 
 
 # -- General configuration ---------------------------------------------------

diff --git a/doped/VASP_sets/HSESet.yaml b/doped/VASP_sets/HSESet.yaml
@@ -1,8 +1,8 @@
 INCAR:  # These settings are only added to the INCAR settings if LHFCALC is not set to False in user_incar_settings
   AEXX: 0.25  # HSE06 assumed by default. Set AEXX in INCAR config settings to tune alpha xc parameter.
   HFSCREEN: 0.208  # correct HSE screening parameter; see https://aip.scitation.org/doi/10.1063/1.2404663
-  # Note this 👆 differs from the Materials Project MPHSERelaxSet default of 0.2! (Will cause systematic
-  # energy shifts in HSE supercell calculations). Change to 0 for PBE0.
+  # 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. Change to 0 for PBE0.
   LHFCALC: true
   PRECFOCK: Fast
   GGA: PE  # underlying GGA functional for HSE/PBE0 is PBE
diff --git a/doped/VASP_sets/RelaxSet.yaml b/doped/VASP_sets/RelaxSet.yaml
@@ -1,5 +1,5 @@
 INCAR:
-  '# May want to change NCORE, KPAR, AEXX, ENCUT, NUPDOWN, ISPIN, POTIM': 'typical variable parameters'
+  '# May want to change NCORE, KPAR, AEXX, ENCUT, IBRION, LREAL, NUPDOWN, ISPIN': 'Typical variable parameters'
   ALGO: "Normal  # Change to All if ZHEGV, FEXCP/F or ZBRENT errors encountered"
   EDIFF_PER_ATOM: 2e-07  # capped at a max EDIFF of 1e-4 for large supercells (N(atoms) > 500)
   EDIFFG: -0.01

diff --git a/doped/analysis.py b/doped/analysis.py
@@ -15,11 +15,11 @@
 import pandas as pd
 from monty.json import MontyDecoder
 from monty.serialization import dumpfn, loadfn
-from pymatgen.analysis.defects.utils import TopographyAnalyzer
 from pymatgen.analysis.structure_matcher import ElementComparator, StructureMatcher
 from pymatgen.ext.matproj import MPRester
 from pymatgen.io.vasp.inputs import Poscar
 from pymatgen.util.string import unicodeify
+from shakenbreak.input import _get_voronoi_nodes
 from tabulate import tabulate
 
 from doped import _ignore_pmg_warnings
@@ -35,7 +35,7 @@
     get_locpot,
     get_outcar,
     get_vasprun,
-    reorder_unrelaxed_structure,
+    reorder_s1_like_s2,
 )
 from doped.vasp import DefectDictSet
 
@@ -168,8 +168,9 @@ def defect_entry_from_paths(
             from the defect calculation outputs (requires `POTCAR`s to be set up
             with `pymatgen`).
         initial_defect_structure (str):
-            Path to the unrelaxed defect structure, if structure matching with the
-            relaxed defect structure(s) fails (rarely required). Default is None.
+            Path to the initial/unrelaxed defect structure. Only recommended for use
+            if structure matching with the relaxed defect structure(s) fails (rare).
+            Default is None.
         skip_corrections (bool):
             Whether to skip the calculation and application of finite-size charge
             corrections to the defect energy (not recommended in most cases).
@@ -272,26 +273,23 @@ def defect_entry_from_paths(
             "Please manually specify defect charge using the `charge_state` argument."
         )
 
-    # Can specify initial defect structure (to help PyCDT find the defect site if
-    # multiple relaxations were required, else use from defect relaxation OUTCAR:
+    # Add defect structure to calculation_metadata, so it can be pulled later on (eg. for Kumagai loader)
+    defect_structure = defect_vr.final_structure.copy()
+    calculation_metadata["defect_structure"] = defect_structure
+
+    # identify defect site, structural information, and create defect object:
+    # Can specify initial defect structure (to help find the defect site if
+    # multiple relaxations were required, else use from defect relaxation OUTCAR):
     if initial_defect_structure:
-        initial_defect_structure = Poscar.from_file(initial_defect_structure).structure.copy()
+        defect_structure_for_ID = Poscar.from_file(initial_defect_structure).structure.copy()
     else:
-        initial_defect_structure = defect_vr.initial_structure.copy()
-
-    # Add initial defect structure to calculation_metadata, so it can be pulled later on
-    # (eg. for Kumagai loader)
-    calculation_metadata["initial_defect_structure"] = initial_defect_structure
-
-    # identify defect site, structural information, and create defect object
+        defect_structure_for_ID = defect_structure.copy()
     try:
-        def_type, comp_diff = get_defect_type_and_composition_diff(
-            bulk_supercell, initial_defect_structure
-        )
+        def_type, comp_diff = get_defect_type_and_composition_diff(bulk_supercell, defect_structure_for_ID)
     except RuntimeError as exc:
         raise ValueError(
             "Could not identify defect type from number of sites in structure: "
-            f"{len(bulk_supercell)} in bulk vs. {len(initial_defect_structure)} in defect?"
+            f"{len(bulk_supercell)} in bulk vs. {len(defect_structure_for_ID)} in defect?"
         ) from exc
 
     # Try automatic defect site detection - this gives us the "unrelaxed" defect structure
@@ -301,7 +299,7 @@ def defect_entry_from_paths(
             defect_site_idx,
             unrelaxed_defect_structure,
         ) = get_defect_site_idxs_and_unrelaxed_structure(
-            bulk_supercell, initial_defect_structure, def_type, comp_diff
+            bulk_supercell, defect_structure_for_ID, def_type, comp_diff
         )
 
     except RuntimeError as exc:
@@ -311,17 +309,15 @@ def defect_entry_from_paths(
         ) from exc
 
     if def_type == "vacancy":
-        defect_site = bulk_supercell[bulk_site_idx]
+        defect_site = guessed_initial_defect_site = bulk_supercell[bulk_site_idx]
     else:
-        if unrelaxed_defect_structure:
-            defect_site = unrelaxed_defect_structure[defect_site_idx]
-        else:
-            defect_site = initial_defect_structure[defect_site_idx]
+        defect_site = defect_structure_for_ID[defect_site_idx]
+        guessed_initial_defect_site = unrelaxed_defect_structure[defect_site_idx]
 
     if unrelaxed_defect_structure:
         if def_type == "interstitial":
             # get closest Voronoi site in bulk supercell to final interstitial site as this is
-            # likely to be the initial interstitial site
+            # likely to be the _initial_ interstitial site
             try:
                 struc_and_node_dict = loadfn("./bulk_voronoi_nodes.json")
                 if not StructureMatcher(
@@ -340,13 +336,7 @@ def defect_entry_from_paths(
                 voronoi_frac_coords = struc_and_node_dict["Voronoi nodes"]
 
             except FileNotFoundError:  # first time parsing
-                topography = TopographyAnalyzer(
-                    bulk_supercell,
-                    bulk_supercell.symbol_set,
-                    [],
-                    check_volume=False,
-                )
-                voronoi_frac_coords = [voronoi_tuple[0] for voronoi_tuple in topography.labeled_sites]
+                voronoi_frac_coords = [site.frac_coords for site in _get_voronoi_nodes(bulk_supercell)]
                 struc_and_node_dict = {
                     "bulk_supercell": bulk_supercell,
                     "Voronoi nodes": voronoi_frac_coords,
@@ -362,22 +352,27 @@ def defect_entry_from_paths(
                 voronoi_frac_coords,
                 key=lambda node: defect_site.distance_and_image_from_frac_coords(node)[0],
             )
-            int_site = unrelaxed_defect_structure[defect_site_idx]
-            unrelaxed_defect_structure.remove_sites([defect_site_idx])
-            unrelaxed_defect_structure.insert(
+            guessed_initial_defect_structure = unrelaxed_defect_structure.copy()
+            int_site = guessed_initial_defect_structure[defect_site_idx]
+            guessed_initial_defect_structure.remove_sites([defect_site_idx])
+            guessed_initial_defect_structure.insert(
                 defect_site_idx,  # Place defect at same position as in DFT calculation
                 int_site.species_string,
                 closest_node_frac_coords,
                 coords_are_cartesian=False,
                 validate_proximity=True,
             )
-            defect_site = unrelaxed_defect_structure[defect_site_idx]
+            guessed_initial_defect_site = guessed_initial_defect_structure[defect_site_idx]
+
+        else:
+            guessed_initial_defect_structure = unrelaxed_defect_structure.copy()
 
-        # Use the unrelaxed_defect_structure to fix the initial defect structure
-        initial_defect_structure = reorder_unrelaxed_structure(
-            unrelaxed_defect_structure, initial_defect_structure
+        # ensure unrelaxed_defect_structure ordered to match defect_structure, for appropriate charge
+        # correction mapping
+        unrelaxed_defect_structure = reorder_s1_like_s2(
+            unrelaxed_defect_structure, defect_structure_for_ID
         )
-        calculation_metadata["initial_defect_structure"] = initial_defect_structure
+        calculation_metadata["guessed_initial_defect_structure"] = guessed_initial_defect_structure
         calculation_metadata["unrelaxed_defect_structure"] = unrelaxed_defect_structure
     else:
         warnings.warn(
@@ -389,7 +384,7 @@ def defect_entry_from_paths(
         "@module": "doped.core",
         "@class": def_type.capitalize(),
         "structure": bulk_supercell,
-        "site": defect_site,
+        "site": guessed_initial_defect_site,
     }  # note that we now define the Defect in the bulk supercell, rather than the primitive structure
     # as done during generation. Future work could try mapping the relaxed defect site back to the
     # primitive cell, however interstitials will be tricky for this...
@@ -403,12 +398,12 @@ def defect_entry_from_paths(
         if not StructureMatcher(
             stol=0.05,
             comparator=ElementComparator(),
-        ).fit(test_defect_structure, unrelaxed_defect_structure):
+        ).fit(test_defect_structure, guessed_initial_defect_structure):
             warnings.warn(
                 f"Possible error in defect object matching. Determined defect: {defect.name} for defect "
                 f"at {defect_path} in bulk at {bulk_path} but unrelaxed structure (1st below) does not "
                 f"match defect.get_supercell_structure() (2nd below):"
-                f"\n{unrelaxed_defect_structure}\n{test_defect_structure}"
+                f"\n{guessed_initial_defect_structure}\n{test_defect_structure}"
             )
 
     defect_entry = DefectEntry(
@@ -569,6 +564,21 @@ def _convert_anisotropic_dielectric_to_isotropic_harmonic_mean(
             # Check compatibility of defect corrections with loaded metadata, and apply
             dp.run_compatibility()
 
+            # check that charge corrections are not negative
+            summed_corrections = sum(
+                val
+                for key, val in dp.defect_entry.corrections.items()
+                if any(i in key.lower() for i in ["freysoldt", "kumagai", "fnv", "charge"])
+            )
+            if summed_corrections < 0:
+                warnings.warn(
+                    f"The calculated finite-size charge corrections for defect at {defect_path} and bulk "
+                    f"at {bulk_path} sum to a negative value of {summed_corrections:.3f}. This is likely "
+                    f"due to some error or mismatch in the defect and bulk calculations, as the defect "
+                    f"charge correction energy should (almost always) be positive. Please double-check "
+                    f"your calculations and parsed results!"
+                )
+
     return dp.defect_entry
 
 
@@ -1020,8 +1030,8 @@ def from_paths(
             from the defect calculation outputs (requires POTCARs to be set up with
             `pymatgen`).
         initial_defect_structure (str):
-            Path to the unrelaxed defect structure, if structure matching with the
-            relaxed defect structure(s) fails.
+            Path to the initial/unrelaxed defect structure. Only recommended for use
+            if structure matching with the relaxed defect structure(s) fails.
         skip_corrections (bool):
             Whether to skip the calculation of finite-size charge corrections for the
             DefectEntry.
@@ -1109,14 +1119,6 @@ def freysoldt_loader(self, bulk_locpot=None):
                 "defect_frac_sc_coords": self.defect_entry.sc_defect_frac_coords,
             }
         )
-        if "unrelaxed_defect_structure" in self.defect_entry.calculation_metadata:
-            self.defect_entry.calculation_metadata.update(
-                {
-                    "initial_defect_structure": self.defect_entry.calculation_metadata[
-                        "unrelaxed_defect_structure"
-                    ],
-                }
-            )
 
     def kumagai_loader(self, bulk_outcar=None):
         """
@@ -1183,11 +1185,11 @@ def _raise_incomplete_outcar_error(outcar_path, dir_type="bulk"):
         bulk_structure = self.defect_entry.bulk_entry.structure
         bulksites = [site.frac_coords for site in bulk_structure]
 
-        defect_structure = self.defect_entry.calculation_metadata["initial_defect_structure"]
+        defect_structure = self.defect_entry.calculation_metadata["unrelaxed_defect_structure"]
         initsites = [site.frac_coords for site in defect_structure]
 
         distmatrix = bulk_structure.lattice.get_all_distances(
-            bulksites, initsites
+            bulksites, initsites  # TODO: Should be able to take this from the defect ID functions?
         )  # first index of this list is bulk index
         min_dist_with_index = [
             [