GalvCalc 1.0.1

From atomic-scale descriptors to macroscopic corrosion polarization curves

GalvCalc models the micro-galvanic corrosion of alloys with coupled anodic dissolution and hydrogen-evolution kinetics. Starting from a bulk crystal structure, it builds the surface models, estimates the electrochemical descriptors (equilibrium potentials, work functions, surface energies and hydrogen adsorption free energies), and assembles multi-phase polarization curves — with built-in anode/cathode area-ratio optimization and alloy-content scans of the corrosion current and potential.

The package is developed alongside the manuscript "GalvCalc: A Framework for Modeling Micro-galvanic Corrosion of Alloys with Coupled Anodic Dissolution and Hydrogen Evolution Kinetics".

Features

Highlights

• Full pipeline: bulk crystal → slab generation → Wulff shape → facet-weighted descriptors → corrosion polarization curves.
• Nernst-equation equilibrium potentials for arbitrary ion combinations (Mg, Fe, Zn, Al, ...), backed by a bundled database of aqueous-ion formation free energies.
• Calibrated Butler–Volmer kinetics for both Mg- and Fe-based systems.
• DFT-grounded second-phase exchange currents for nine common Mg-alloy intermetallics (Mg₁₇Al₁₂, Mg₂Al₃, MgZn₂, LaMg₁₂, CaMg₂, Y₅Mg₂₄, NdMg₃, Mg₂Si, CeMg₁₂).
• CGCNN prediction of surface energies and work functions, and TabPFN-based prediction of hydrogen adsorption free energies.

Installation

pip install GalvCalc

Core dependencies are installed by default. The ml extra adds the pre-trained predictors (CGCNN / TabPFN), and the defects extra (Python ≥ 3.10) enables the substitution-based doping workflow in GalvCalc.anode.

Requirements: Python >= 3.9; numpy, scipy, pandas, matplotlib, pymatgen, PyYAML, sympy, tqdm, joblib, openpyxl. The ml extra adds torch, scikit-learn, tabpfn.

Quick start

All structure files used below (Mg2Si.vasp, Mg.vasp, surface/*.vasp) and the demo DFT adsorption dataset (adsorption_analysis.csv) ship with the package under GalvCalc/examples/, so the snippets run fully offline from that folder (as in the bundled notebook).

1. Equilibrium potential

from GalvCalc.core.structures import Bulk

# Load a bulk structure from a POSCAR / CIF / .vasp file
bulk = Bulk.from_file("Mg2Si.vasp")

# Single-ion reaction: Mg -> Mg[2+] + 2 e-
Ee = Bulk.get_equilibrium_potential(
    ions="Mg[2+]",
    ion_numbers=[1],
    energy_formation=0.0,
)
print(f"Mg E_eq = {Ee:.3f} V vs. SHE")

# Multi-ion compound dissolution: MgZn2 -> Mg[2+] + 2 Zn[2+] + 6 e-
Ee2 = Bulk.get_equilibrium_potential(
    ions=["Mg[2+]", "Zn[2+]"],
    ion_numbers=[1, 2],
    energy_formation=-0.24,
)
print(f"MgZn2 E_eq = {Ee2:.3f} V vs. SHE")

2. Surfaces, ML-predicted facet properties and Wulff construction

from GalvCalc.core.structures import Bulk
from GalvCalc.cathode.surfaces import SurfaceProperties

bulk = Bulk.from_file("Mg2Si.vasp")

# Batch-generate all surface terminations up to max_index = 1
surface_props = SurfaceProperties.from_bulk_structure(
    bulk,
    max_index=1,
    min_slab_size=15.0,
    min_vacuum_size=15.0,
    center_slab=True,
    symmetrize=True,
    max_normal_search=1,
)

# CGCNN-predicted facet surface energies / work functions (needs the `ml` extra)
df_pred = surface_props.get_predicted_surface_properties()

# Wulff shape from the predicted facet properties
wulff = surface_props.wulff_construct(
    surface_energies=surface_props.surface_energies_dict,
    work_functions=surface_props.work_functions_dict,
    output_dir="demo_output",
    save_plot=True,
    save_csv=True,
)

3. Hydrogen adsorption and exchange-current estimation

The demo DFT adsorption dataset is bundled as adsorption_analysis.csv; the facet-weighted hydrogen-evolution exchange current follows from the adsorption free energy through the BEP-type relation of the manuscript:

import pandas as pd
from GalvCalc.cathode import ic0_mg

df_ads = pd.read_csv("adsorption_analysis.csv")
mg2si = df_ads[df_ads["formula"] == "Mg2Si"]
print(mg2si[["miller_index", "termination", "ads_position", "Eads", "workfunction"]].head())

# Mg2Si (110) H6 site: G_ads = 0.028 eV, work function = 3.7378 eV
i0 = ic0_mg(delta_G=0.028, wf=3.7378)
print(f"ic0_mg = {i0:.3e} A/cm^2")

from GalvCalc.cathode.hydrogen import AdsorptionManager

adsorption_manager = AdsorptionManager(surface_props.surfaces)
results = adsorption_manager.H_adsorption_analysis(
    adsorbate="H",
    site_indices="all",
    output_dir="H_adsorption_analysis",
    include_summary=True,
    include_dataframe=True,
    include_visualization=True,
)

formatted_data = adsorption_manager.get_formatted_data()
adsorption_manager.save_formatted_data("adsorption_data.json")

AdsorptionManager.H_adsorption_analysis locates the adsorption sites on every surface, builds the H-adsorbed slabs, predicts the hydrogen adsorption free energies with the bundled TabPFN model and exports POSCAR files plus CSV/JSON summaries. This path needs the ml extra.

4. Substitutional doping of anode surfaces

from GalvCalc.core.structures import Bulk, Surface
from GalvCalc.anode import SurfaceDopingManager

mg_bulk = Bulk.from_file("Mg.vasp")
mg_surface = Surface.from_bulk(bulk_structure=mg_bulk, miller_index=(0, 0, 1))
manager = SurfaceDopingManager(mg_surface, "Mg_001")

# Substitute surface Mg by Zn, Al and Y (second layer by default)
manager.batch_dope(["Zn", "Al", "Y"], save_to_file=False)
for name, info in manager.doping_info.items():
    print(f"{name}: host {info.host_element} -> {info.dopant_element}, "
          f"site {info.site_index}, layer {info.layer}, "
          f"multiplicity {info.multiplicity}, depth {info.depth:.2f} A")

manager.set_property("Mg_001", "work_function", 3.72)
manager.set_property("Mg_001", "vacancy_energy", 0.84)
manager.set_property("Mg_001_Zn", "work_function", 3.75)
manager.set_property("Mg_001_Zn", "vacancy_energy", 0.79)
manager.set_property("Mg_001_Al", "work_function", 3.74)
manager.set_property("Mg_001_Al", "vacancy_energy", 0.84)
manager.set_property("Mg_001_Y", "work_function", 3.29)
manager.set_property("Mg_001_Y", "vacancy_energy", 0.78)

df = manager.calculate_electrochemical_properties(E00=-2.37, ia00=1e-5)
print(df[["surface_name", "E0_calculated", "ia0_calculated", "dopant_element"]])

Set save_to_file=True to export the doped structures as POSCAR files into doped_surfaces/; the full doping record (host/dopant, site, layer, symmetry multiplicity) is available from manager.doping_info. This workflow builds on pymatgen's defect framework and needs the defects extra (Python ≥ 3.10).

5. Polarization curves and area-ratio optimization

from GalvCalc.polarization import (
    ElectrodeParameters, Composition, plot_single_polarization,
    mg_second_phase_example, plot_mg_second_phases,
)
from GalvCalc.polarization.area_ratio import (
    AreaRatioAnalyzer, create_example_parameters,
)

# Fe-based system: n = 1, exponent (alpha_a + 1) * n, prefactor 2
anode = ElectrodeParameters(4.1e-8, -0.44, 0.5, name="Fe/Fe2+", kinetic_form="fe")
cathode = ElectrodeParameters(7.9e-8, -0.059, 0.5, name="H+/H2", kinetic_form="fe")
comp = Composition("Fe", anode=anode, cathodes=[cathode], area_ratios=[1, 1])
fig = plot_single_polarization(comp, reference_electrode="SHE")

# Mg-based system: nine intermetallic second-phase cathodes on one plot
fig_mg = plot_mg_second_phases(mg_second_phase_example())

# Anode/cathode area-ratio optimization
params = create_example_parameters()
fig, optimal_ratio, i_max = AreaRatioAnalyzer(reference_electrode="SCE").plot_area_ratio_analysis(params)
print(f"optimal anode area ratio = {optimal_ratio:.2%}")

The second-phase cathodic exchange currents in mg_second_phase_example() are DFT-calibrated for nine Mg intermetallics; supply your own values to study other phases.

6. Alloy-content scan

import numpy as np
from GalvCalc.polarization import (
    wt_to_vol_fraction, two_phase_kinetics, scan_corrosion_vs_content,
)

# Any intermetallic second phase works: swap the phase formula, densities
# and kinetic parameters. Mg-Nd / Mg3Nd is the example used here.
rho_matrix, rho_phase = 1.8, 3.58  # Mg matrix, Mg3Nd second phase
nd_wt = np.linspace(0, 6, 13)      # Nd content (wt%)

def kinetics_at_content(w):
    vol_fraction = wt_to_vol_fraction(
        w, "Mg3Nd", rho_matrix=rho_matrix, rho_phase=rho_phase
    )
    return two_phase_kinetics(
        vol_fraction=vol_fraction,
        anode_i0=6e-23 if w <= 0 else 3e-23,
        anode_equilibrium_potential=-2.37,
        cathode_i0s=[10 ** -8.1, 10 ** -8.732],
        cathode_equilibrium_potential=-0.61,
    )

df_scan = scan_corrosion_vs_content(
    nd_wt,
    params_fn=kinetics_at_content,
    domain=(0.0, 1.0),
    n_ratios=51,
    n_potentials=1500,
)
df_scan[["content", "max_log10_i_corr", "E_corr_at_max"]].round(4)

Example notebooks

Runnable notebooks are shipped under GalvCalc/examples/:

• galvcalc_demo.ipynb — a step-by-step walkthrough covering equilibrium potentials, slab generation, adsorption-site schematics, Wulff construction, facet-weighted exchange currents, anode doping, polarization curves, area-ratio optimization and alloy-content scans.

Example structures (*.vasp, surface/*.vasp) and the demo DFT adsorption dataset (adsorption_analysis.csv) are bundled in the same folder.

License

Distributed under the MIT License.

Contact

Gaoning Shi — gaoning_shi@sjtu.edu.cn