Skip to content

Basic GPAW support - #149

Open
sblisesivdin wants to merge 37 commits into
SMTG-Bham:developfrom
sblisesivdin:main
Open

Basic GPAW support#149
sblisesivdin wants to merge 37 commits into
SMTG-Bham:developfrom
sblisesivdin:main

Conversation

@sblisesivdin

Copy link
Copy Markdown
  • New doped/gpaw.py is implemented for PW and LCAO workflows.
  • Test script is located at tests/test_gpaw.py.
  • examples/Graphene_with_GPAW is added with 2 Python scripts and a CIF file. The first py script is examples/Graphene_with_GPAW/gpaw_workflow.py, which creates the necessary files, and the second is examples/Graphene_with_GPAW/run_and_analyze.py, which runs all files and prints the analysis results.
  • There is also a documentation for GPAW implementation at docs/GPAW_Support.rst
  • Also, there are some minor corrections at doped/core.py and doped/corrections.py.

@coderabbitai

coderabbitai Bot commented Jan 31, 2026

Copy link
Copy Markdown

Important

Review skipped

Auto reviews are disabled on base/target branches other than the default branch.

Please check the settings in the CodeRabbit UI or the .coderabbit.yaml file in this repository. To trigger a single review, invoke the @coderabbitai review command.

⚙️ Run configuration

Configuration used: Repository UI

Review profile: CHILL

Plan: Pro Plus

Run ID: 8a01aad6-e68f-4cdf-8a6e-7db8f587e6e8

You can disable this status message by setting the reviews.review_status to false in the CodeRabbit configuration file.

Use the checkbox below for a quick retry:

  • 🔍 Trigger review

Walkthrough

Introduces comprehensive GPAW integration to doped library with a new module providing defect relaxation setup (GPAWDefectRelaxSet), calculation parsing (GPAWParser and GPAWDefectsParser), and utility functions. Includes bug fixes for Kumagai correction empty-array handling and improved defect region radius calculation. Adds documentation, workflow examples, and test coverage.

Changes

Cohort / File(s) Summary
Documentation and API Reference
docs/GPAW_Support.rst, docs/doped.gpaw.rst, docs/doped.rst, docs/index.rst
New GPAW support documentation page with installation, workflow steps, and usage examples; API reference and toctree entries for doped.gpaw module.
GPAW Core Module
doped/gpaw.py
New comprehensive module introducing GPAWDefectRelaxSet for generating GPAW input with CIF and relax.py script, GPAWParser for loading and extracting properties from .gpw files, GPAWDefectsParser for batch-parsing multiple defect calculations with auto-detection of bulk folders, and utility functions for site potentials and planar-averaged potentials.
Corrections Logic
doped/corrections.py
Changes defect_region_radius calculation from pydefect default to derived reciprocal lattice norm when not specified; affects sampling region for eFNV correction.
Core Logic Bug Fix
doped/core.py
Adds safeguard in Kumagai correction to handle empty sampling arrays by conditionally computing correction_error or setting to NaN.
Example Workflow
examples/Graphene_with_GPAW/gpaw_workflow.py, examples/Graphene_with_GPAW/graphene4x4.cif, examples/Graphene_with_GPAW/run_and_analyze.py
New example demonstrating end-to-end GPAW defect workflow for graphene; includes structure file, defect generation, input generation, calculation execution, and thermodynamic analysis.
Tests
tests/test_gpaw.py
Unit tests for GPAWDefectRelaxSet with default and custom settings, LCAO mode, and GPAWParser with mocked calculations.

Sequence Diagram(s)

sequenceDiagram
    participant User
    participant GPAWDefectRelaxSet
    participant FileSystem
    participant Structure as Structure/DefectEntry

    User->>GPAWDefectRelaxSet: __init__(defect_entry, charge_state, gpaw_settings)
    GPAWDefectRelaxSet->>Structure: Extract structure info
    User->>GPAWDefectRelaxSet: write_input(output_path)
    GPAWDefectRelaxSet->>FileSystem: Write structure.cif
    GPAWDefectRelaxSet->>FileSystem: Generate and write relax.py script
    FileSystem-->>User: Input files ready
Loading
sequenceDiagram
    participant GPAWDefectsParser
    participant FileSystem as File System
    participant GPAWParser
    participant Bulk as Bulk Calc
    participant Defect as Defect Calcs
    participant Corrections as Corrections Engine
    participant DefectEntry

    GPAWDefectsParser->>FileSystem: Auto-detect bulk folder
    GPAWDefectsParser->>Bulk: Initialize GPAWParser(bulk_path)
    GPAWParser->>Bulk: Load .gpw file
    Bulk-->>GPAWParser: Computed structure, potentials
    GPAWDefectsParser->>FileSystem: Iterate defect folders
    loop For each defect
        GPAWDefectsParser->>Defect: Initialize GPAWParser(defect_path)
        GPAWParser->>Defect: Load .gpw file
        Defect-->>GPAWParser: Computed structure, potentials
        GPAWDefectsParser->>Corrections: Apply Kumagai correction
        Corrections-->>GPAWDefectsParser: Corrected energy
        GPAWDefectsParser->>DefectEntry: Construct DefectEntry
    end
    GPAWDefectsParser-->>GPAWDefectsParser: parse_all() returns Dict[DefectEntry]
Loading

Estimated code review effort

🎯 4 (Complex) | ⏱️ ~70 minutes

Poem

🐰 A new GPAW garden we grow,
With relax.py scripts all aglow,
Defects parsed with perfect care,
Potentials flowing through the air,
Corrections balanced, graphene spreads fair!

🚥 Pre-merge checks | ✅ 2 | ❌ 1
❌ Failed checks (1 warning)
Check name Status Explanation Resolution
Docstring Coverage ⚠️ Warning Docstring coverage is 77.78% which is insufficient. The required threshold is 80.00%. Write docstrings for the functions missing them to satisfy the coverage threshold.
✅ Passed checks (2 passed)
Check name Status Explanation
Title check ✅ Passed The title 'Basic GPAW support' directly and clearly summarizes the main objective of adding GPAW integration to the doped package, matching the primary changes across multiple files.
Description check ✅ Passed The description comprehensively covers all major changes including the new gpaw.py module, tests, examples, documentation, and minor corrections, all directly related to the changeset.

✏️ Tip: You can configure your own custom pre-merge checks in the settings.

✨ Finishing Touches
🧪 Generate unit tests (beta)
  • Create PR with unit tests

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands.

@coderabbitai coderabbitai Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actionable comments posted: 1

🤖 Fix all issues with AI agents
In `@examples/Graphene_with_GPAW/run_and_analyze.py`:
- Around line 19-48: The print statements use unnecessary f-strings (Ruff F541)
— remove the f-prefix on the string literals in the prints inside the
run_and_analyze workflow: change the prints that reference "relaxed.gpw exists,
skipping calculation.", "gpaw_output.txt exists but relaxed.gpw does not.
Skipping likely non-conserving/failed calculation.", "relax.py not found in
{folder}, skipping." and the two messages inside the try/except ("Calculation
completed.", "Calculation failed: {e}" — keep f only when interpolating e) so
that only the messages that actually include variables remain formatted; update
the prints near the subprocess.run block and the FileNotFoundError handler
accordingly, referencing the existing variables folder and exception variable e
and ensure any literal messages are plain strings (no leading f).
🧹 Nitpick comments (10)
doped/gpaw.py (4)

156-170: Consider adding error handling for grid point lookup.

The gd.get_nearest_grid_point and subsequent indexing could fail if the atom position is outside expected bounds or if the grid descriptor is not properly initialized.

💡 Suggested improvement
 def _get_site_potentials_from_calc(calc) -> np.ndarray:
     """
     Helper to extract site potentials from a GPAW calculator.
     """
     atoms = calc.get_atoms()
     v_ext = calc.get_electrostatic_potential()
     gd = calc.hamiltonian.finegd
     
     site_potentials = []
     for atom in atoms:
-        indices = gd.get_nearest_grid_point(atom.position)
-        val = v_ext[tuple(indices % gd.N_c)]
-        site_potentials.append(val)
+        try:
+            indices = gd.get_nearest_grid_point(atom.position)
+            val = v_ext[tuple(indices % gd.N_c)]
+            site_potentials.append(val)
+        except (IndexError, AttributeError) as e:
+            raise RuntimeError(f"Failed to get potential at atom position {atom.position}: {e}")
         
     return np.array(site_potentials)

262-273: Edge case: empty eigenvalues list could cause issues.

If get_eigenvalues returns empty lists for all k-points and spins, the sorted(energies) will be empty, and the list comprehensions on lines 269-270 will produce empty results, leading to vbm = efermi and cbm = efermi with band_gap = 0. While the fallback is reasonable, a warning might be helpful for debugging.

💡 Suggested improvement
         energies = sorted(energies)
+        if not energies:
+            import warnings
+            warnings.warn("No eigenvalues found; band gap, VBM, and CBM set to Fermi level.")
+            return 0.0, efermi, efermi, efermi
         # Identify VBM and CBM based on efermi
         vbm = max([e for e in energies if e <= efermi]) if any(e <= efermi for e in energies) else efermi

397-405: Silent exception swallowing in charge state parsing.

The try-except-pass block silently ignores parsing errors when guessing charge state from folder names. This could hide issues where folder naming conventions are misunderstood.

🔧 Proposed fix to log parsing failures
             try:
                 if "_" in folder:
                     suffix = folder.rsplit("_", 1)[-1]
                     if suffix.startswith(("+", "-")):
                         charge_state = int(suffix)
-            except Exception:
-                pass
+            except (ValueError, IndexError):
+                pass  # Keep default charge_state=0 for non-standard naming

416-425: Consider more specific exception handling.

The broad Exception catches could mask unexpected errors. Consider catching more specific exceptions for better debuggability.

🔧 Proposed improvement
                 # Apply Kumagai correction if possible
                 if self.dielectric is not None and charge_state != 0:
                     try:
                         defect_entry.get_kumagai_correction()
-                    except Exception as e:
+                    except (ValueError, RuntimeError) as e:
                         print(f"Warning: Kumagai correction failed for {folder}: {e}")
                 
                 defect_dict[defect_entry.name] = defect_entry
-            except Exception as e:
+            except (FileNotFoundError, ValueError, RuntimeError) as e:
                 print(f"Failed to parse {folder}: {e}")
examples/Graphene_with_GPAW/gpaw_workflow.py (2)

58-60: Unused loop variable entry in the defect listing loop.

The loop on line 59 iterates over items but only uses name. Consider using _ for the unused variable or removing the .items() call.

💡 Proposed fix
     print("Generated defects:")
-    for name, entry in generator.defect_entries.items():
+    for name in generator.defect_entries:
         print(f" - {name}")

34-39: Consider adding validation before species replacement.

The replace_species({"N": "C"}) call assumes N atoms exist in the structure. If the CIF file doesn't contain N, this will silently do nothing, potentially leading to confusion.

💡 Suggested improvement
 bulk_structure = defect_structure.copy()
+if "N" not in [str(sp) for sp in bulk_structure.composition.elements]:
+    print("Warning: No N atoms found in structure to replace with C")
 bulk_structure.replace_species({"N": "C"})
tests/test_gpaw.py (1)

1-106: Consider adding tests for GPAWDefectsParser and get_gpaw_defect_entry.

The current tests cover GPAWDefectRelaxSet and GPAWParser well, but get_gpaw_defect_entry and GPAWDefectsParser are not tested. These functions have more complex logic that would benefit from test coverage.

Would you like me to help draft additional tests for get_gpaw_defect_entry and GPAWDefectsParser?

examples/Graphene_with_GPAW/run_and_analyze.py (3)

12-40: Make core count configurable instead of hard-coded.

Line 12 hard-codes CORES = 8, which is brittle across machines. Consider pulling from an env var (or CLI) with a safe default.

Proposed change
-CORES = 8 # Set number of cores for MPI
+CORES = int(os.environ.get("GPAW_CORES", "8"))  # Set number of cores for MPI

24-43: Ensure the gpaw_output.txt sentinel is actually written.

The skip logic depends on gpaw_output.txt, but stdout is currently discarded. Unless relax.py writes that file itself, the sentinel won’t exist. Redirecting stdout/stderr to the file makes the skip logic reliable and preserves logs for debugging.

Proposed change
-        subprocess.run(cmd, cwd=folder, check=True, stdout=subprocess.DEVNULL)
+        log_path = os.path.join(folder, "gpaw_output.txt")
+        with open(log_path, "w") as log:
+            subprocess.run(
+                cmd,
+                cwd=folder,
+                check=True,
+                stdout=log,
+                stderr=subprocess.STDOUT,
+            )

49-50: Generic Exception catches should be more specific or include re-raising for unexpected errors.

The blanket except Exception blocks (lines 49–50 and 97–100) mask unexpected failures. However, note that line 97–100 already includes traceback.print_exc(), which does log the full error detail.

For lines 49–50 (in the batch processing loop): Consider whether swallowing exceptions silently is intentional. If failures should be visible, either log more detail (like the traceback) or re-raise to fail fast. If best-effort processing is desired, add logging similar to the traceback approach used at lines 97–100.

For lines 97–100: Already includes traceback printing, which addresses the core concern. If proceeding silently is acceptable here, the current implementation is reasonable; otherwise, consider re-raising after the traceback.

Comment on lines +19 to +48
# Check if already run (relaxed.gpw exists)
if os.path.exists(os.path.join(folder, "relaxed.gpw")):
print(f" relaxed.gpw exists, skipping calculation.")
return

# Check if failed/non-conserving (gpaw_output.txt exists but relaxed.gpw doesn't)
if os.path.exists(os.path.join(folder, "gpaw_output.txt")):
print(f" gpaw_output.txt exists but relaxed.gpw does not. Skipping likely non-conserving/failed calculation.")
return

if not os.path.exists(relax_script):
print(f" relax.py not found in {folder}, skipping.")
return

try:
# Construct command
if CORES > 1:
# Use mpirun for parallel execution
# Note: gpaw python is often recommended but python usually works if gpaw installed correctly
cmd = ["mpirun", "-np", str(CORES), "gpaw", "python", "relax.py"]
else:
cmd = ["python", "relax.py"]

# Run in the folder to handle relative paths correctly
subprocess.run(cmd, cwd=folder, check=True, stdout=subprocess.DEVNULL)
print(f" Calculation completed.")
except subprocess.CalledProcessError as e:
print(f" Calculation failed: {e}")
except FileNotFoundError:
print(f" Execution failed. Check if mpirun/gpaw are in your PATH.")

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Remove empty f-strings (Ruff F541).

Lines 21, 26, 44, and 48 use f-strings without placeholders, which triggers lint errors.

Proposed change
-        print(f"  relaxed.gpw exists, skipping calculation.")
+        print("  relaxed.gpw exists, skipping calculation.")
@@
-        print(f"  gpaw_output.txt exists but relaxed.gpw does not. Skipping likely non-conserving/failed calculation.")
+        print("  gpaw_output.txt exists but relaxed.gpw does not. Skipping likely non-conserving/failed calculation.")
@@
-        print(f"  Calculation completed.")
+        print("  Calculation completed.")
@@
-        print(f"  Execution failed. Check if mpirun/gpaw are in your PATH.")
+        print("  Execution failed. Check if mpirun/gpaw are in your PATH.")
🧰 Tools
🪛 Ruff (0.14.14)

[error] 21-21: f-string without any placeholders

Remove extraneous f prefix

(F541)


[error] 26-26: f-string without any placeholders

Remove extraneous f prefix

(F541)


[error] 43-43: subprocess call: check for execution of untrusted input

(S603)


[error] 44-44: f-string without any placeholders

Remove extraneous f prefix

(F541)


[error] 48-48: f-string without any placeholders

Remove extraneous f prefix

(F541)

🤖 Prompt for AI Agents
In `@examples/Graphene_with_GPAW/run_and_analyze.py` around lines 19 - 48, The
print statements use unnecessary f-strings (Ruff F541) — remove the f-prefix on
the string literals in the prints inside the run_and_analyze workflow: change
the prints that reference "relaxed.gpw exists, skipping calculation.",
"gpaw_output.txt exists but relaxed.gpw does not. Skipping likely
non-conserving/failed calculation.", "relax.py not found in {folder}, skipping."
and the two messages inside the try/except ("Calculation completed.",
"Calculation failed: {e}" — keep f only when interpolating e) so that only the
messages that actually include variables remain formatted; update the prints
near the subprocess.run block and the FileNotFoundError handler accordingly,
referencing the existing variables folder and exception variable e and ensure
any literal messages are plain strings (no leading f).

@kavanase

Copy link
Copy Markdown
Member

Hi @sblisesivdin! This is fantastic, thank you very much for contributing to doped!

I will set aside some time to go through this properly, but some general thoughts for now:

  • Thank you very much for adding tests – very useful. I think it would really help robustness and confidence if we can also add some quick tests – primarily for the charge corrections – using the same structures/functional as those in the VASP tests. I imagine this would be a pain to do for any hybrid DFT tests we have, but having some like-for-like tests with GGA data (e.g. from examples/MgO/Defects/Pre_Calculated_Results/) would be great to ensure consistency in the charge correction implementations.
  • I see for the parsing tests you use mocking, and then a script for a full workflow test. I think mocking is useful for unit testing, but experience has shown us that testing with real data is best for identifying bugs and making maintenance easier. Given that you have the full workflow test script, would it be possible to add some of the calculation results from that script as test data, and add some parsing tests with that?
  • There is some draft code to implement Quantum Espresso support in Adding QE support to doped.  #133, which is currently being worked on offline. In that code, the defect parser classes are abstracted to make support for different DFT engines more modular and easier to maintain and expand in future. Not asking anything from you here, just noting that in the (hopefully near) future we will likely rearrange some of this code to work with an abstracted & modular structure similar to Adding QE support to doped.  #133.

Thank you for your work on this!

@kavanase
kavanase changed the base branch from main to develop January 31, 2026 18:19
Comment thread doped/corrections.py Outdated
defect_region_radius = 0.5 * np.min(
[1 / np.linalg.norm(lattice.reciprocal_lattice.matrix[i]) for i in range(3)]
)
defect_region_radius *= np.pi

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this radius multiplied by pi?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Looking again I realise this is because it's going to/from the reciprocal lattice)

@kavanase kavanase Feb 1, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This new code returns half the defect region radius as with the previous code (with the previous code matching the pydefect implementation) for the isotropic case, and then relatively smaller again for anisotropic cases, where the difference is changing from the pydefect code which uses the max distance between parallel cell planes, to instead be the shortest distance between parallel planes (-> inscribed sphere radius).

This change (with a 2pi factor) makes sense to me. I think it would be good to raise this as an issue with the pydefect code to (1) fix any issues with the current implementation and ensure consistency between doped and pydefect for the eFNV correction, and (2) to run it by the author of the eFNV correction (Yu Kumagai, who maintains pydefect). Would you mind doing this @sblisesivdin? Thank you!
https://github.com/kumagai-group/pydefect/blob/a9d96d3c996f2e173a1ba16945c26a85e432f5ef/pydefect/cli/vasp/make_efnv_correction.py#L96

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm... thinking about it again, the purpose of defect_region_radius is just to determine the sampling region (everything outside of this), which is then just used to get the cell potential at a point as far away from the defect (to determine the alignment-like correction, if I'm correct).
So for a highly-anisotropic system, I think you would still want it closer to the largest distance between parallel planes rather than shortest, I would think? I can imagine in some rare, anisotropic cases (e.g. graphene), taking the max could cause this to result in a sampling region with no atoms (which I guess is what you found?)... As an alternative approach, one could get the shortest distances to all other atoms in the cell, and set the cutoff as one which gives (at least) 10 atoms outside this radius? Or within ~3 Å of the furthest defect-to-atom distance, or some other tolerance like that?

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for being late. I opened an issue at pydefect as you suggested.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @sblisesivdin!
Just noting that this is the same implementation, just slightly cleaner (with numpy vectorisation) and including the 2pi rather than 1pi factor mentioned above:

defect_region_radius = 2*np.pi * np.min(1/np.linalg.norm(lattice.reciprocal_lattice.matrix, axis=1))
defect_region_radius *= 0.5  # half the shortest distance between parallel planes

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@sblisesivdin I think this factor-of-2 bug noted here is still not fixed?

Comment thread doped/gpaw.py Outdated

# Write structure to a file
structure_filename = "structure.cif"
self.defect_supercell.to(filename=os.path.join(output_path, structure_filename), fmt="cif")

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nitpicking, but from what I remember with pymatgen.to(... fmt='cif') it writes the unsymmetrized (P1) CIF (in all cases), but using the CifWriter(..., symprec=0.01) class one can write the symmetrized CIF, which I guess is preferred for most cases?

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Solved with d9a5a0a

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ups. Not solved with the commit d9a5a0a. Regarding the CifWriter(symprec=0.01) suggestion: I implemented this with d9a5a0a, but we actually have to omit the symprec argument for the relax sets! If symprec is passed, pymatgen automatically reduces perfect bulk supercells down to their primitive unit cells. This breaks the parsing step when it attempts to map the 31-atom defect cell back to the bulk. I kept it as the standard CifWriter(supercell) to ensure it preserves the user's supercell matrix, as in commit 74dff1d.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok! Can we add a comment to the code here to note this? (So it doesn't get changed in the future when one has forgotten about this issue) – a good example of the importance of tests :)

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comment added with dd442ba

Comment thread doped/gpaw.py Outdated
Comment thread doped/gpaw.py Outdated
Comment thread docs/GPAW_Support.rst Outdated
-----------------------

GPAW calculations of charged defects require finite-size corrections to account for periodic
image interactions. ``doped`` supports the **Kumagai (eFNV)** correction for GPAW.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can also still support the FNV correction too right? As we have the planar averaged potentials. Most people will and should use eFNV, but it would be good to keep FNV as an option for backwards compatibility/comparability etc

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, the FNV correction is still supported. Documentation is updated as can be seen at 61e6892

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great. With the GPAW parser the Kumagai correction is currently hard-coded though.
In practice, this should always be preferred, and the user can manually apply the FNV correction after parsing if preferred, with defect_entry.get_freysoldt_correction() I believe? Can we note this in the parser docstring (that eFNV correction is used by default, but FNV can be used with defect_entry.get_freysoldt_correction()), and could we also add tests (using the MgO example data?) of the FNV correction with GPAW outputs?
Thanks!

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK. I have added a docstring in doped/gpaw.py and I have added a new test function in tests/test_gpaw.py for MgO with FNV correction. Everything is in 0187141

@kavanase

Copy link
Copy Markdown
Member
  • You show some nice MWEs in the documentation page. Would you mind adding a short Python notebook, which shows how to use doped with GPAW all through Python if possible? Which we could then include in the docs like the other tutorials. I think this would be very nice to illustrate how one can use doped with GPAW

Comment thread doped/gpaw.py Outdated
sblisesivdin added a commit to sblisesivdin/doped that referenced this pull request Apr 2, 2026
Replaces the single-point grid lookup in `_get_site_potentials_from_calc` with a volumetric average over the atomic core sphere.

Previously, the potential was extracted at the exact nearest grid point, making the calculation highly sensitive to the GPAW grid resolution (`h` parameter). This update introduces a volumetric mask that averages the electrostatic potential (`v_ext`) over a specified `core_radius` (default 1.0 Angs) around each atom, safely handling periodic boundary conditions via the minimum image convention.

This aligns the GPAW parsing behavior with VASP and `sxdefectalign`, ensuring numerically stable and grid-independent site potentials.

Addressing review feedback on PR SMTG-Bham#149.
@sblisesivdin

Copy link
Copy Markdown
Author

Hi @kavanase. Thanks for the comments and I am sorry for the late reply. You and @Ath088 are completely right—taking the potential at the exact grid point makes it far too sensitive to the GPAW grid resolution (h parameter) and ignores the spatial distribution of the core.

Since PR #133 is still in draft mode and those utility functions might still be evolving, I wanted to avoid creating a dependency on an unmerged branch for now.

Instead, I've pushed a commit that updates _get_site_potentials_from_calc() to use a self-contained volumetric averaging method. It now calculates the mean electrostatic potential over a spherical core region (defaulting to 1.0 Angs) around each atom and safely handles periodic boundaries via the minimum image convention. This successfully mimics the VASP and sxdefectalign behavior to provide grid-independent site potentials.

Once PR #133 is merged and the codebase is further modularized, I agree it would be great to refactor this to use the centralized core sphere averaging utility!

Let me know if the new implementation looks good to you.

@kavanase

kavanase commented Apr 2, 2026

Copy link
Copy Markdown
Member

Hi @sblisesivdin,
Great, thanks for your work on this!

Yes that makes sense about avoiding dependence on a draft PR for now. The updated implementation definitely looks like a good improvement to me, though I would guess that a weighted average around the atomic sites (like the Gaussian averaging used in sxdefectalign) would be better than a flat mean over some radius around the atomic sites. This has been implemented by @Ath088 here: https://github.com/Walser52/doped/blob/e8012e515c3a0a7d1641f513b05a66126952706b/doped/utils/parsing.py#L2444-L2543

So it could make sense to just copy over this code for now?

You can see at https://github.com/Ath088/doped_QE/blob/main/examples/qe_function_tests.ipynb that @Ath088 has added some like-for-like tests between VASP and Quantum Espresso for the MgO defect examples. As mentioned above, I think having some like-for-like comparisons like this between GPAW and VASP (and Quantum Espresso) would be very useful for testing the charge correction functionality here, with the MgO GGA defect examples likely being the best choice for comparisons to start. (We may also want to have similar tests for an anisotropic system in the future, to ensure the proper handling of anisotropic dielectrics and the defect_region_radius updates you added – we already have a few anisotropic test systems for the VASP tests).

@sblisesivdin

Copy link
Copy Markdown
Author

Hi @kavanase,

Thanks for the feedback and the links!

Regarding the site potential extraction, you make a great point about using a Gaussian-weighted average rather than a flat mean. I will review @Ath088's implementation from the doped_QE branch and integrate that Gaussian averaging logic into _get_site_potentials_from_calc.

Regarding the testing: I completely agree. Setting up a like-for-like comparison between GPAW and VASP using the MgO GGA defect examples is the best way to validate the charge correction functionality. I will start working on generating the GPAW outputs for those MgO structures to build the test cases. I hope I can do.

Also, an update on the pydefect side:
Yu Kumagai has replied to our issue regarding defect_region_radius. He agreed that evaluating alternative definitions for anisotropic cells makes sense, noting the trade-off between having enough sampling points vs. staying far enough from the defect.

His proposed solution is to update pydefect to allow the user to choose the radius definition (min, max, or average), and he mentioned he is leaning towards making the average distance the new default.

Given this, how would you prefer to handle the defect_region_radius calculation in this PR? Should I update our doped logic to use the average distance instead of the minimum to anticipate his changes, or should I go ahead and implement a radius_method parameter in doped so we are fully synced with the upcoming pydefect architecture?

@sblisesivdin

Copy link
Copy Markdown
Author

Hi again @kavanase

We have an update. Kumagai just replied again. He offered to implement the changes in pydefect himself within the next week! He plans to update the defect_region_radius parameter to accept strings (min, max, average) alongside the existing float input.

Since he is handling it upstream, we aren't blocked. Because pydefect currently accepts a float, doped can calculate the radius that we want and just pass the float to pydefect in the meantime.

For this PR, would you prefer:

  1. Leave the defect_region_radius calculation as it is (using the minimum/inscribed sphere).
  2. Change the internal calculation to use the average distance (which Kumagai is leaning toward as his new default).

Once his update is merged, we can refactor doped to just pass the string arguments directly to pydefect.

BTW, I am working on MgO calculations with GPAW :) I am hopeful.

@kavanase

kavanase commented Apr 3, 2026

Copy link
Copy Markdown
Member

Great! Thanks for following up with this @sblisesivdin.

For the defect region radius discussion, for now I think I would stick with your usage of the minimum inscribed sphere in this PR, while the discussion is still going on in pydefect on how exactly that will be updated, and once that is finalised and integrated we can decide the best way to integrate in doped.

Great to hear about the GPAW calculations with MgO! 😃

@sblisesivdin

Copy link
Copy Markdown
Author

Hi @kavanase,

As you know, Yu Kumagai has released the update in pydefect v0.9.12! (He closed the issue).

At this point, your #149 (comment) comment may have expired due to pydefect's new release.

Therefore, what to do next? In corrections.py, there is a doped _make_efnv_correction() function, which is a modified version of pydefect's make_efnv_correction() function. Do we still need this? Maybe just use:

from pydefect.cli.vasp.make_efnv_correction import (
     HalfMaxFaceDistanceDefectRegion,
     HalfMinFaceDistanceDefectRegion)

and

def _get_defect_region_object(radius_method: str, sample_radius_ratio: float = 0.9):
    """Translates user string into the corresponding pydefect region object."""
    method = radius_method.lower()
    if method == "max":
        return HalfMaxFaceDistanceDefectRegion(sample_radius_ratio=sample_radius_ratio)
    elif method == "min":
        return HalfMinFaceDistanceDefectRegion(sample_radius_ratio=sample_radius_ratio)
    # Add other mappings as necessary...
    else:
        raise ValueError(f"Unknown radius method: {radius_method}")

I hope I am using this new pydefect usage correctly.

@kavanase

Copy link
Copy Markdown
Member

Hi @sblisesivdin,
Thanks for checking! While removing hacked functions is definitely preferable in most cases, the latest version of this function in doped develop has some other alterations to make it significantly more efficient (using doped site mapping functions), so it's probably best to keep it in doped.
This means we can still have the fully flexible defect region radius choice, but also using that class from pydefect as default, if we think it's the best default choice.

For next steps, I think prioritizing like-for-like tests with the MgO defects (for which we now have correction tests for VASP and Quantum Espresso) would be best. Then can focus on streamlining the DefectsParser for GPAW and getting this PR merged. 😃

@sblisesivdin

Copy link
Copy Markdown
Author

Ok. I have added the MgO example. I haven't control it with the other solvers' results yet. Sorry for any inconvenience. I will do that and add my comments in the future. But, for now, I can share the latest version of files with you. MgO example's output is like:

    --- Step 1: Running Bulk Calculation ---
    --- Step 2: Running Defect Calculations ---
    --- Step 3: Parsing Results with Doped ---
    Parsing v_Mg_-1...
    Calculated Kumagai (eFNV) correction is 0.096 eV
    Parsing v_Mg_+1...
    Calculated Kumagai (eFNV) correction is 0.300 eV
    Parsing v_Mg_0...
    Parsing v_Mg_-2...
    Calculated Kumagai (eFNV) correction is 0.589 eV
    Estimated error in the Kumagai (eFNV) charge correction for defect v_Mg_Oh_O2.10_-2 is 0.080 eV (i.e. which is greater than the `error_tolerance`: 0.050 eV). You may want to check the accuracy of the correction by plotting the site potential differences (using `defect_entry.get_kumagai_correction()` with `plot=True`). Large errors are often due to unstable or shallow defect charge states (which can't be accurately modelled with the supercell approach; see https://doped.readthedocs.io/en/latest/Tips.html#perturbed-host-states-shallow-defects). If this error is not acceptable, you may need to use a larger supercell for more accurate energies.
    
    Parsed 4 defects.
    
    --- Step 4: Thermodynamic Analysis ---
    Chemical potentials not present for elements: ['Mg']. Assuming zero chemical potentials for these elements! (Absolute formation energies will likely be very inaccurate)
    
    Defect: v_Mg_Oh_O2.10_-1
      Charge: -1
      Supercell Energy: -484.0953 eV
      Corrections: {'kumagai_charge_correction': np.float64(0.09605757784246306)}
      Formation Energy (at VBM, no chempots): 4.3256 eV
    Defect: v_Mg_Oh_O2.10_+1
      Charge: 1
      Supercell Energy: -484.0953 eV
      Corrections: {'kumagai_charge_correction': np.float64(0.30000927595554494)}
      Formation Energy (at VBM, no chempots): 14.7101 eV
    Defect: v_Mg_Oh_O2.10_0
      Charge: 0
      Supercell Energy: -484.0953 eV
      Corrections: {}
      Formation Energy (at VBM, no chempots): 9.3198 eV
    Defect: v_Mg_Oh_O2.10_-2
      Charge: -2
      Supercell Energy: -484.0953 eV
      Corrections: {'kumagai_charge_correction': np.float64(0.5887542543871135)}
      Formation Energy (at VBM, no chempots): -0.2719 eV
    
    Full Automation Test Complete.

@kavanase

kavanase commented May 2, 2026

Copy link
Copy Markdown
Member

Hi @sblisesivdin, this looks promising! For v_Mg^+1, is this the relaxed or unrelaxed structure? If the former, then this value (~0.3 eV) seems to match well with that for doped-VASP/QE (https://github.com/Ath088/doped_QE/blob/main/examples/qe_function_tests.ipynb).

For the GPAW tests, it seems you are running these as self-contained scripts. Is it possible to save the output files and then use these for the tests? (To avoid the same GPAW calculations needing to be re-run every time this part of doped is tested).
This is what we do with the VASP/QE tests, and would also allow us to avoid using mocking in test_gpaw.py here.

@sblisesivdin

Copy link
Copy Markdown
Author

Hi @kavanase,

Thank you! Yes, the v_Mg defects (including the +1 charge state) were generated from the fully relaxed bulk MgO structure. It is really exciting to see that the ~0.3 eV correction aligns so well with the VASP and QE baselines!

Regarding the tests: I completely agree! In my latest commit, I updated the testing approach to do exactly what you suggested. I generated a coarse set of output files locally and saved them to tests/data/gpaw_mgo_test/ (specifically for bulk and v_Mg_-2).

The test_gpaw.py suite now simply parses these static saved files to verify the site potential extraction and the Kumagai correction math. No live GPAW calculations are run during the CI tests, and all the old mocking has been completely removed.

Comment thread doped/gpaw.py Outdated
@kavanase

Copy link
Copy Markdown
Member

Hi @sblisesivdin, great, thank you!
Would be great to polish off this PR so we can merge to main.

There are some comments on the code I added before, and an additional one now above, would you be able to address these please?

For the tests, they look good. It would be great if we could expand the charge correction tests just to have a few more datapoints beyond v_Mg^+2, to include the other charge states, and ideally the case of v_Mg^+1 unrelaxed defect structure and relaxed, to test that difference. Would you mind adding these? I guess this is relatively cheap to run, as you had most of these test cases in the previous version which ran the calculations all within the script?

Would it also be possible to include a test using the graphene example you give, especially as it previously caused the failure mode you noted regarding the defect region radius handling?

I think with these, we will then be very close to being able to merge this branch to main.

Just to note; for the site potentials handling, we will likely want to merge the core interpolation function (which now is common to QE and GPAW parsing) in the future, but don't need to worry about it now. There were some slight robustness updates to the site potentials code in the QE development branch (https://github.com/Ath088/doped_QE/blob/main/doped/utils/parsing.py#L2462), if you could copy over the new code for this function to ensure the correct handling – thanks!

@kavanase

Copy link
Copy Markdown
Member

Something I missed from your previous comment:

    --- Step 4: Thermodynamic Analysis ---
    Chemical potentials not present for elements: ['Mg']. Assuming zero chemical potentials for these elements! (Absolute formation energies will likely be very inaccurate)
    
    Defect: v_Mg_Oh_O2.10_-1
      Charge: -1
      Supercell Energy: -484.0953 eV
      Corrections: {'kumagai_charge_correction': np.float64(0.09605757784246306)}
      Formation Energy (at VBM, no chempots): 4.3256 eV
    Defect: v_Mg_Oh_O2.10_+1
      Charge: 1
      Supercell Energy: -484.0953 eV
      Corrections: {'kumagai_charge_correction': np.float64(0.30000927595554494)}
      Formation Energy (at VBM, no chempots): 14.7101 eV
    Defect: v_Mg_Oh_O2.10_0
      Charge: 0
      Supercell Energy: -484.0953 eV
      Corrections: {}
      Formation Energy (at VBM, no chempots): 9.3198 eV
    Defect: v_Mg_Oh_O2.10_-2
      Charge: -2
      Supercell Energy: -484.0953 eV
      Corrections: {'kumagai_charge_correction': np.float64(0.5887542543871135)}
      Formation Energy (at VBM, no chempots): -0.2719 eV
    
    Full Automation Test Complete.

Here all the v_Mg charge states give the same supercell energy, implying that the charge state input is not being respected somehow? It would be good to test both the total value (as printed here) and the alignment-like term (also stored in the eFNV correction output), as here I suspect it's given the same alignment-like term (due to being duplicate calculations) with only the PC-like term changing. Including correction plot tests like we have done for VASP & QE would catch this, and would also be great to include if possible (otherwise I can add these later, if the issue is resolved).

@sblisesivdin

Copy link
Copy Markdown
Author

Hi @kavanase,

Wow, incredibly sharp eye on the supercell energies! You are exactly right—because the charge parameter wasn't properly getting passed to the GPAW() calculator initialization in the generated relax.py scripts, every calculation defaulted to neutral, resulting in the exact same energy and PC-like terms. Give me some time, I am super busy in these days. I want to work on it (and other problems) with a clear mind.

Sorry for the late reply.

@kavanase

Copy link
Copy Markdown
Member

No worries at all, thanks @sblisesivdin!

@sblisesivdin

Copy link
Copy Markdown
Author

Hi @kavanase,

You were absolutely right about the supercell energies. The charge parameter wasn't properly getting passed to the GPAW() calculator initialization, so it was silently defaulting to a neutral state for every calculation.

I've just pushed a few commits that fix all of the issues:

  • Fixed GPAWDefectRelaxSet to properly inject the charge state into the GPAW() object. The supercell energies are now different across charge states.
  • Updated GPAWParser to pull calc.parameters.get('charge') directly from the .gpw file, falling back to the folder name only if it's missing.
  • Replaced the interpolation logic with the robust implementation from the doped_QE branch.
  • I regenerated the coarse MgO .gpw test files with the actual charge states and expanded test_gpaw.py to test all charge states dynamically. I also added the unrelaxed +1 state as you requested.
  • Also, some cosmetic changes in the output are provided.

As I saw, everything works fine now. But, you are catching problems very well :) If anything seems wrong, I can work on it.

@kavanase

Copy link
Copy Markdown
Member

Hi @sblisesivdin, thanks!

I think some of the main things to tackle would be:

There are some remaining review comments on the code, would you be able to address these please?

Would it also be possible to include a test using the graphene example you give, especially as it previously caused the failure mode you noted regarding the defect region radius handling?

You show some nice MWEs in the documentation page. Would you mind adding a short Python notebook, which shows how to use doped with GPAW all through Python if possible? Which we could then include in the docs like the other tutorials. I think this would be very nice to illustrate how one can use doped with GPAW

Then I can have an in-depth look through all the code and aim to merge. :)

@sblisesivdin

Copy link
Copy Markdown
Author

The Jupyter tutorial was a new thing for me. Normally, I have not used Jupyter for my studies. I can not make gpaw run on multicore with Jupyter. Therefore, I included a workflow of preparing input files etc. However, parsing is done on ready gpw files in tests/data/gpaw-mgo-test folder. Done in f94d47f and 98e5a0f

Comment thread doped/gpaw.py Outdated
) -> DefectEntry:
"""
Convenience function to create a DefectEntry from GPAW directories.
Assumes 'relaxed.gpw' exists in both directories.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should ideally auto-detect the output .gpw file in the directory, rather than hard-coding a specific name, with a preference order if multiple .gpw files present – see the _find_calc_outputs usage:

def _get_calculation_folders_for_parsing(

Comment thread doped/gpaw.py Outdated
else:
self.bulk_path = bulk_path

def parse_all(self) -> Dict[str, DefectEntry]:

@kavanase kavanase Jun 23, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ideally the GPAWDefectsParser would have the same structure/workflow as the VASP DefectsParser, meaning no parse_all method (it parses all defects with the init call: GPAWDefectsParser(...)

Comment thread doped/gpaw.py
Comment on lines +477 to +483
defect_entry = get_gpaw_defect_entry(
defect_dir,
self.bulk_path,
dielectric=self.dielectric,
charge_state=charge_state,
bulk_parser=bulk_parser,
)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This re-parses the bulk site potentials every time, which I guess can be time-consuming?
In the VASP DefectsParser this is why we parse the bulk site potentials once and then re-use with each defect (see doped.analysis). Can we do that here for efficiency?

@kavanase

Copy link
Copy Markdown
Member

It would be great if we could show an example of the defect formation energy diagram for MgO, showing that most of the full defect workflow can now be performed with GPAW. Like this:
image

You already have the defect supercells for this, so all that would be missing is the Mg and O chemical potentials. Would you be able to prepare these for this example (using the workflow for how it is done via CompetingPhases in the MgO tutorial)? Thank you for your work on this @sblisesivdin!!

@kavanase

kavanase commented Jun 23, 2026

Copy link
Copy Markdown
Member

Noting that some review comments are appearing as 'outdated' after I rebased to develop and pushed formatting updates, but are not resolved.

The rendered docs can be viewed here: https://doped--149.org.readthedocs.build/en/149/

@sblisesivdin

Copy link
Copy Markdown
Author

After your pre-commit formating and clean up, I managed to lose a huge chunk of my new changes due to a rebase and aborting it... Totally my bad but it is annoying you guess. I think i will take a break and continue later.

@kavanase

Copy link
Copy Markdown
Member

Ah sorry to hear that @sblisesivdin, that's a pain! If you happen to be using an IDE like PyCharm/VSCode/Cursor (most of which have good education discounts/free pro accounts) they usually save local history, so you can revert and not lose any progress. Or ofc if you use Time Machine with Mac or similar regular file snapshots.
I'm possibly misunderstanding, but aborting a rebase/merge should put you back to where you were without losing any local changes (unless it was a hard reset?). git reflog can be useful for tracking what has saved by git at any point.

@sblisesivdin

Copy link
Copy Markdown
Author

I am using GitKraken and gnome text editor on Linux, and normally they are working fine for me. The thing is my unstaged files are gone. I did not push them. And now they are gone. Before, I did that rebase/merge thing, I pull and few hundred changes came, and I tried to unstage few changes (don't ask me why, I think it is the wrong thing), but when I understand it is the wrong thing, I abort the rebase. But, all my unstaged files are gone. Luckily, I have a hard copy of few days ago, i will continue from that point in time, but let me cool down firstly :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants