Using doped with GPAW

This tutorial demonstrates how to use the doped Python API to generate defects, set up GPAW calculators, and parse the thermodynamic results completely within a Python environment.

Note: To keep this tutorial fast and accessible without requiring a High-Performance Computing (HPC) cluster, we will parse pre-calculated MgO (Magnesium Oxide) data included in the doped test suite.

import os
from pymatgen.core.structure import Structure
from pymatgen.core.lattice import Lattice
from doped.generation import DefectsGenerator
from doped.gpaw import GPAWDefectRelaxSet, GPAWDefectsParser

# 1. Generate the optimized bulk structure directly (MgO rock-salt)
lattice = Lattice.cubic(4.21)
bulk_structure = Structure.from_spacegroup("Fm-3m", lattice, ["Mg", "O"], [[0, 0, 0], [0.5, 0.5, 0.5]])

# 2. Generate Defects
generator = DefectsGenerator(bulk_structure)
print(f"Generated {len(generator.defect_entries)} defects.")

Writing GPAW Input Scripts

We can use GPAWDefectRelaxSet to automatically generate relax.py scripts for our bulk and defect supercells, passing in our desired GPAW calculator parameters.

Executing the cell below will create actual calculation folders (calc_bulk and calc_...) in your current directory so you can inspect the generated inputs.

# Define GPAW settings
gpaw_settings = {
    "mode": {"name": "pw", "ecut": 400},
    "xc": "PBE",
    "kpts": {"size": (2, 2, 2), "gamma": True},
}

# Write Bulk Input
bulk_dir = "calc_bulk"
bulk_set = GPAWDefectRelaxSet(generator.bulk_supercell, charge_state=0, gpaw_settings=gpaw_settings)
bulk_set.write_input(bulk_dir)
print("Written input for bulk supercell")

# Write Defect Inputs (Just the first 3 to keep the directory clean)
for name, entry in list(generator.defect_entries.items())[:3]:
    folder_name = f"calc_{name}"
    defect_set = GPAWDefectRelaxSet(entry, charge_state=entry.charge_state, gpaw_settings=gpaw_settings)
    defect_set.write_input(folder_name)
    print(f"Written input for {name}")

print("\nGPAW input generation complete! In a real workflow, you would now submit these folders to your HPC job scheduler.")

Parsing and Charge Corrections

For this tutorial, we will load the pre-calculated gpw output files from the doped test suite.

By default, doped calculates the Kumagai (eFNV) finite-size charge correction.

Note on Defect Radius: doped natively supports flexible defect region definitions. You can optionally pass radius_method="max" (default), "min", or "average" to the Kumagai function to adjust the sampling region.

Note on FNV: If you prefer the Freysoldt (FNV) correction, you can manually calculate it after parsing using defect_entry.get_freysoldt_correction().

# Locate the pre-calculated test data. Change the folder address if necessary.
data_dir = "../tests/data/gpaw_mgo_test"
bulk_dir = os.path.join(data_dir, "bulk")

# Initialize parser (using a dummy dielectric for MgO = 10.0)
parser = GPAWDefectsParser(
    output_path=data_dir, 
    bulk_path=bulk_dir,
    dielectric=10.0
)

# Parse all completed calculations
defect_dict = parser.parse_all()

# Print Thermodynamic Results and compare Kumagai vs Freysoldt
for name, entry in defect_dict.items():
    if entry.charge_state == 0:
        continue # Skip neutral defects
        
    # 1. Recalculate Kumagai with custom region sampling ('max', 'min', or 'average')
    entry.get_kumagai_correction(radius_method="max")
    
    # 2. Manually calculate Freysoldt (FNV) for comparison
    entry.get_freysoldt_correction()
    
    print(f"Defect: {name}")
    print(f"  Charge: {entry.charge_state}")
    print(f"  Supercell Energy: {entry.sc_entry.energy:.4f} eV")
    
    kumagai_corr = entry.corrections.get('kumagai_charge_correction', 0.0)
    fnv_corr = entry.corrections.get('freysoldt_charge_correction', 0.0)
    
    print(f"  Kumagai (eFNV) Correction (max): {kumagai_corr:.4f} eV")
    print(f"  Freysoldt (FNV) Correction:      {fnv_corr:.4f} eV\n")