diff --git a/.gitignore b/.gitignore index 617fe48..ecaef68 100644 --- a/.gitignore +++ b/.gitignore @@ -34,6 +34,7 @@ examples/hBN/negf_output_k20/self_energy/self_energy_leadL.h5 examples/hBN/negf_output_k20/self_energy/self_energy_leadR.h5 examples/hBN/negf_output_k50/self_energy/* examples/hBN/negf_output_k70/self_energy/* +examples/hBN/conductance_* examples/CNT/negf_profiling examples/CNT/output* examples/CNT/*long* diff --git a/dpnegf/runner/NEGF.py b/dpnegf/runner/NEGF.py index b8eb8d0..2654117 100644 --- a/dpnegf/runner/NEGF.py +++ b/dpnegf/runner/NEGF.py @@ -19,6 +19,9 @@ from dpnegf.utils.constants import Boltzmann, eV2J from dpnegf.utils.make_kpoints import kmesh_sampling_negf from dpnegf.utils.band_edge import validate_fermi_in_band_gap +from dpnegf.utils.energy_grid import build_energy_grid, uniform_nodes +from dpnegf.utils.conductance import integrate_conductance +from dpnegf.utils.argcheck import validate_energy_options from dpnegf.negf.poisson_init import Grid,Interface3D,Dirichlet,Dielectric from dpnegf.negf.scf_method import PDIISMixer,DIISMixer,BroydenFirstMixer,BroydenSecondMixer,AndersonMixer @@ -42,12 +45,12 @@ def __init__(self, model: torch.nn.Module, structure: Union[AtomicData, ase.Atoms, str], ele_T: float, - emin: float, emax: float, espacing: float, - density_options: dict, - unit: str, - scf: bool, poisson_options: dict, - stru_options: dict,eta_lead: float,eta_device: float, - block_tridiagonal: bool, + emin: float=None, emax: float=None, espacing: float=None, + density_options: dict=None, + unit: str="eV", + scf: bool=False, poisson_options: dict=None, + stru_options: dict=None,eta_lead: float=1e-5,eta_device: float=0., + block_tridiagonal: bool=False, self_energy_options: Optional[dict] = None, rgf_options: Optional[dict] = None, hs_cache: Optional[dict] = None, @@ -55,9 +58,40 @@ def __init__(self, e_fermi: float=None, results_path: Optional[str]=None, plot_blocks: Optional[bool]=False, AtomicData_options: Optional[dict]=None, + energy_grid: Optional[dict]=None, + conductance_options: Optional[dict]=None, **kwargs): - + # Use an explicitly supplied energy grid; otherwise preserve + # the legacy uniform-grid configuration. + self.energy_grid_options = ( + dict(energy_grid) + if energy_grid is not None + else { + "method": "uniform", + "emin": emin, + "emax": emax, + "espacing": espacing, + } + ) + self.conductance_options = dict(conductance_options or {}) + self.out_conductance = (output_options or {}).get( + "conductance",False,) + self.spin_degeneracy = self.conductance_options.get( + "spin_degeneracy") + if self.spin_degeneracy is None: + self.spin_degeneracy = 1 if hasattr(model, "soc_param") else 2 + validate_energy_options( + { + "energy_grid": self.energy_grid_options, + "conductance_options": self.conductance_options, + "output_options": output_options or {}, + "ele_T": ele_T, + "unit": unit, + "scf": scf, + "stru_options": stru_options, + } + ) # self.model = model # No need to set model as property for memory saving self.results_path = results_path self.cdtype = torch.complex128 @@ -350,7 +384,7 @@ def __init__(self, # geting the output settings out_opts = dict(output_options or {}) - self.out_tc = out_opts.get("tc", False) + self.out_tc = out_opts.get("tc", False) or self.out_conductance self.out_dos = out_opts.get("dos", False) self.out_density = out_opts.get("density", False) self.out_potential = out_opts.get("potential", False) @@ -435,10 +469,39 @@ def generate_energy_grid(self): if self.out_current: cal_int_grid = True - if self.out_dos or self.out_tc or self.out_current_nscf or self.out_ldos: + if (self.out_dos or self.out_tc or self.out_current_nscf or self.out_ldos + or self.energy_grid_options.get("method") == "clenshaw_curtis"): # Energy gird is set relative to Fermi level - self.uni_grid = torch.linspace(start=self.emin, end=self.emax, steps=int((self.emax-self.emin)/self.espacing)) + if self.out_conductance or self.energy_grid_options.get("method") == "clenshaw_curtis": + if abs(self.chemiPot["lead_L"]-self.chemiPot["lead_R"]) > 5e-4: + if self.out_conductance: + raise ValueError("Conductance calculation is not supported for non-zero bias case in this version.") + if self.energy_grid_options.get("method") == "clenshaw_curtis": + raise ValueError("Clenshaw-Curtis energy grid is not supported for non-zero bias case in this version.") + lead = self.conductance_options.get("lead", "lead_L") + labels = self.conductance_options.get("mu", ["Ef"]) + reference = float(self.deviceprop.E_ref) + centers = [] + for label in labels: + if label == "Ef": + absolute = self.e_fermi[lead] + elif label == "Ev": + absolute = self.E_v[lead] + elif label == "Ec": + absolute = self.E_c[lead] + else: + # label is a number that relates to the reference energy + centers.append(float(label)) + continue + centers.append(absolute - reference) + self.quadrature_grid = build_energy_grid( + self.energy_grid_options, centers, self.ele_T, labels) + self.uni_grid = torch.from_numpy(self.quadrature_grid.energies.copy()) + else: + self.uni_grid = torch.from_numpy(uniform_nodes(self.energy_grid_options)) + self.energy_grid = self.uni_grid + if cal_pole and self.density_options["method"] == "Ozaki": self.poles, self.residues = ozaki_residues(M_cut=self.density_options["M_cut"]) self.poles = 1j* self.poles * self.kBT + self.deviceprop.lead_L.chemiPot - self.deviceprop.chemiPot @@ -718,7 +781,13 @@ def negf_compute(self,scf_require=False,Vbias=None): assert scf_require is not None, "scf_require should be set to True or False" self.out['k']=[];self.out['wk']=[] - if hasattr(self, "uni_grid"): self.out["uni_grid"] = self.uni_grid + if hasattr(self, "uni_grid"): + self.out["uni_grid"] = self.uni_grid # compatibility alias, possibly non-unifrom + self.out["energy_grid"] = self.uni_grid # compatibility alias, possibly non-unifrom + self.out["E_ref"] = float(self.deviceprop.E_ref) + self.out["energy_grid_options"] = dict(self.energy_grid_options) + if hasattr(self, "quadrature_grid"): + self.out["energy_grid_metadata"] = self.quadrature_grid.metadata() # Self-Energy Calculaiton or Loading self.prepare_self_energy(scf_require) @@ -1001,7 +1070,17 @@ def negf_compute(self,scf_require=False,Vbias=None): if scf_require==False: self.out["k"] = np.array(self.out["k"]) - self.out['T_avg'] = torch.tensor(self.out['wk']) @ torch.stack(list(self.out["T_k"].values())).cpu() + if "T_k" in self.out: + self.out['T_avg'] = torch.tensor(self.out['wk']) @ torch.stack(list(self.out["T_k"].values())).cpu() + if self.out_conductance: + records = integrate_conductance(self.out["T_avg"], self.quadrature_grid, self.spin_degeneracy) + for record in records: + record["E_ref_eV"] = float(self.deviceprop.E_ref) + record["mu_absolute_eV"] = record["mu_eV"] + record["E_ref_eV"] + self.out["conductance"] = records + with open(os.path.join(self.results_path, "conductance.json"), 'w') as f: + json.dump(records, f, indent=2) + log.info("Landauer conductance (G/G0): %s", [r["G_over_G0"] for r in records]) # TODO:check the following code for multiple k points calculation if self.out_current_nscf: self.out["BIAS_POTENTIAL_NSCF"], self.out["CURRENT_NSCF"] = self.compute_current_nscf(self.uni_grid, self.out["T_avg"]) diff --git a/dpnegf/tests/test_conductance.py b/dpnegf/tests/test_conductance.py new file mode 100644 index 0000000..e89a03c --- /dev/null +++ b/dpnegf/tests/test_conductance.py @@ -0,0 +1,365 @@ +import copy +import json +from pathlib import Path + +import numpy as np +import pytest +import torch +from scipy.integrate import quad +from scipy.special import expit + +from dpnegf.utils.argcheck import normalize_run +from dpnegf.utils.energy_grid import ( + build_energy_grid, + clenshaw_curtis_nodes, + thermal_energy, + uniform_nodes, +) +from dpnegf.utils.conductance import ( + conductance_clenshaw_curtis, + conductance_uniform, + integrate_conductance, + G0_S, +) + +ROOT = Path(__file__).resolve().parents[2] +DATA = ROOT / "dpnegf/tests/data/test_negf/test_negf_run" + + +@pytest.mark.parametrize("temperature", [10, 300, 1000]) +@pytest.mark.parametrize("n", [17, 33, 65]) +def test_cc_mass_symmetry_nesting(temperature, n): + offsets, weights = clenshaw_curtis_nodes(n, temperature) + + assert np.isfinite(offsets).all() + np.testing.assert_array_equal(offsets, -offsets[::-1]) + assert np.all(weights > 0) + + mass = np.tanh(0.4 / (2 * thermal_energy(temperature))) + assert weights.sum() == pytest.approx(mass, abs=1e-14) + + refined, _ = clenshaw_curtis_nodes(2 * n - 1, temperature) + np.testing.assert_allclose( + offsets, + refined[::2], + atol=1e-15, + rtol=0, + ) + + grid = build_energy_grid( + { + "method": "clenshaw_curtis", + "num_points": n, + }, + temperature=temperature, + ) + result = conductance_clenshaw_curtis(np.ones(n), grid) + + assert result["G_over_G0"] == pytest.approx(mass) + assert result["G_S"] == pytest.approx(G0_S * mass) + + +@pytest.mark.parametrize("shape", ["smooth", "narrow_peak", "step"]) +def test_physical_spectra_against_quad(shape): + kbt = thermal_energy(300) + + if shape == "smooth": + spectrum = lambda e: 1 + 0.5 * np.cos(e / 0.03) + elif shape == "narrow_peak": + spectrum = lambda e: ( + 0.002**2 / ((e - 0.02)**2 + 0.002**2) + ) + else: + spectrum = lambda e: np.asarray(e >= 0.017, dtype=float) + + kernel = lambda e: ( + expit(-abs(e / kbt)) + * (1 - expit(-abs(e / kbt))) + / kbt + ) + + reference = quad( + lambda e: spectrum(e) * kernel(e), + -0.4, + 0.4, + points=[0.017, 0.02], + epsabs=1e-11, + )[0] + + energies = np.linspace(-0.4, 0.4, 16001) + result = conductance_uniform(energies, spectrum(energies)) + + assert result["G_over_G0"] == pytest.approx( + reference, + rel=1e-3, + ) + + grid = build_energy_grid( + { + "method": "clenshaw_curtis", + "num_points": 1025, + } + ) + cc = conductance_clenshaw_curtis( + spectrum(grid.energies), + grid, + ) + + assert cc["G_over_G0"] == pytest.approx( + reference, + rel=0.01, + ) + + +def test_multi_target_and_k_weighting(): + options = {"method": "clenshaw_curtis"} + grid = build_energy_grid(options, [-0.1, 0.1, -0.1]) + + assert len(grid.energies) == 66 + + spectrum = np.stack( + [ + 1 + np.cos(grid.energies * 10), + 0.5 + grid.energies**2, + ] + ) + k_weights = np.array([0.3, 0.7]) + + results = integrate_conductance(spectrum, grid) + averaged = integrate_conductance(k_weights @ spectrum, grid) + + for result, avg in zip(results, averaged): + assert ( + k_weights @ result["G_over_G0"] + ) == pytest.approx(avg["G_over_G0"]) + + assert averaged[0] == averaged[2] + + shifted = build_energy_grid(options, [0.2, 0.4, 0.2]) + values = np.stack( + [ + 1 + np.cos((shifted.energies - 0.3) * 10), + 0.5 + (shifted.energies - 0.3)**2, + ] + ) + + for a, b in zip( + results, + integrate_conductance(values, shifted), + ): + np.testing.assert_allclose( + a["G_over_G0"], + b["G_over_G0"], + ) + + spin = integrate_conductance( + spectrum, + grid, + spin_degeneracy=1, + ) + np.testing.assert_allclose( + spin[0]["G_S"], + np.array(results[0]["G_S"]) / 2, + ) + + +def config(): + return json.loads( + (DATA / "negf_chain_new.json").read_text() + ) + + +def test_legacy_grid_unchanged_and_new_precedence(): + old = config() + task = old["task_options"] + + expected = torch.linspace( + task["emin"], + task["emax"], + int( + (task["emax"] - task["emin"]) + / task["espacing"] + ), + ) + + normalized = normalize_run(copy.deepcopy(old)) + actual = uniform_nodes( + normalized["task_options"]["energy_grid"] + ) + + np.testing.assert_array_equal(actual, expected.numpy()) + assert actual.dtype == expected.numpy().dtype + + old["task_options"]["energy_grid"] = { + "method": "clenshaw_curtis" + } + old["task_options"]["out_current_nscf"] = False + + task = normalize_run(old)["task_options"] + + assert task["energy_grid"] == { + "method": "clenshaw_curtis", + "num_points": 33, + "half_width": 0.4, + } + assert "emin" not in task + + +@pytest.mark.parametrize( + "change", + [ + {"ele_T": 0}, + {"unit": "Hartree"}, + {"scf": True}, + {"conductance_options": {"mu": []}}, + {"conductance_options": {"mu": ["bad"]}}, + { + "energy_grid": { + "method": "clenshaw_curtis", + "num_points": 32, + } + }, + {"output_options": {"density": True}}, + ], +) +def test_invalid_configuration(change): + raw = config() + raw["task_options"]["energy_grid"] = { + "method": "clenshaw_curtis" + } + raw["task_options"].update(change) + + with pytest.raises(ValueError): + normalize_run(raw) + + +def test_bad_quadrature_inputs(): + with pytest.raises(ValueError, match="cache precision"): + build_energy_grid( + {"method": "clenshaw_curtis"}, + [0.0, 1e-10], + ) + + with pytest.raises(ValueError): + conductance_uniform( + [0, 1, 1.1], + [1, 1, 1], + ) + + with pytest.raises(ValueError): + conductance_clenshaw_curtis(np.ones(32)) + + with pytest.raises(ValueError): + conductance_clenshaw_curtis(np.full(33, np.nan)) + + +@pytest.mark.parametrize( + "device", + [ + "cpu", + pytest.param( + "cuda", + marks=pytest.mark.skipif( + not torch.cuda.is_available(), + reason="CUDA unavailable", + ), + ), + ], +) +def test_chain_runner_cache_and_batching( + tmp_path, + monkeypatch, + device, +): + from dpnegf.entrypoints.run import run + import importlib + + runner = importlib.import_module("dpnegf.runner.NEGF") + raw = config() + task = raw["task_options"] + + task["energy_grid"] = {"method": "clenshaw_curtis"} + task["conductance_options"] = { + "mu": [0.0, 0.05, 0.0] + } + task["output_options"] = {"conductance": True} + + for key in list(task): + if key.startswith("out_"): + del task[key] + + task["rgf_options"] = { + "device": device, + "e_batch_size": 8, + } + + cache = tmp_path / "self_energy" + task["self_energy_options"] = { + "cache": { + "save_path": str(cache), + }, + "numba_jit": False, + "parallel": { + "n_workers": 1, + "blas_threads": 1, + }, + } + + def execute(name): + path = tmp_path / (name + ".json") + path.write_text(json.dumps(raw)) + out = tmp_path / name + + run( + str(path), + str(DATA / "nnsk_C_new.json"), + str(DATA / "chain.vasp"), + str(out), + None, + ) + + return torch.load( + out / "results/negf.out.pth", + weights_only=False, + ) + + first = execute("cold") + + assert len(first["energy_grid"]) == 66 + assert first["conductance"][0]["G_over_G0"] == pytest.approx( + 1, + rel=1e-4, + ) + assert first["conductance"][0] == first["conductance"][2] + + task["self_energy_options"]["cache"]["use_saved"] = True + task["rgf_options"]["e_batch_size"] = 1 + + def no_precompute(*args, **kwargs): + raise AssertionError( + "cached run must not precompute self energy" + ) + + monkeypatch.setattr( + runner, + "compute_all_self_energy", + no_precompute, + ) + + second = execute("warm") + + torch.testing.assert_close( + first["T_avg"], + second["T_avg"], + rtol=1e-9, + atol=1e-10, + ) + + grid = build_energy_grid( + {"method": "clenshaw_curtis"}, + [0, 0.05, 0], + ) + post = integrate_conductance(second["T_avg"], grid) + + for a, b in zip(post, second["conductance"]): + assert a["G_S"] == pytest.approx(b["G_S"]) \ No newline at end of file diff --git a/dpnegf/utils/argcheck.py b/dpnegf/utils/argcheck.py index a7ab431..609c327 100644 --- a/dpnegf/utils/argcheck.py +++ b/dpnegf/utils/argcheck.py @@ -88,6 +88,81 @@ def tbtrans_negf(): +def energy_grid_options(): + return Variant("method", [ + Argument("uniform", dict, [ + Argument("emin", [int, float], optional=False, doc="Lower energy relative to E_ref."), + Argument("emax", [int, float], optional=False, doc="Upper energy relative to E_ref."), + Argument("num_points", [int, None], optional=True, default=None, + doc="Number of nodes including endpoints; mutually exclusive with espacing."), + Argument("espacing", [int, float, None], optional=True, default=None, + doc="Legacy parameter: point count is int((emax-emin)/espacing), including endpoints."), + ]), + Argument("clenshaw_curtis", dict, [ + Argument("num_points", int, optional=True, default=33, + doc="Odd node count >=3 per target. 17/33/65 nodes are nested at fixed temperature and window."), + Argument("half_width", [int, float], optional=True, default=0.4, + doc="Half window in eV about each mu. Tails are not extrapolated or renormalized."), + ]), + ], optional=True, default_tag="uniform") + + +def conductance_options(): + return [ + Argument("mu", list, optional=True, default=["Ef"], + doc="Nonempty list of Ef/Ev/Ec or numeric eV relative to E_ref. Does not change Hamiltonian or occupancy."), + Argument("lead", str, optional=True, default="lead_L", + doc="lead_L or lead_R for Ef/Ev/Ec. Band edges are calculated automatically when requested."), + Argument("spin_degeneracy", [int, None], optional=True, default=None, + doc="1 for explicit spin, 2 for spin-degenerate transmission. None uses model SOC convention."), + ] + + +def validate_energy_options(task): + """Cross-field validation shared by normalized JSON and the Python API.""" + import numpy as np + from dpnegf.utils.energy_grid import uniform_nodes, clenshaw_curtis_nodes + + grid = task["energy_grid"] + method = grid.get("method", "uniform") + conductance = task.get("output_options", {}).get("conductance", False) + if method == "uniform": + uniform_nodes(grid) + elif method == "clenshaw_curtis": + clenshaw_curtis_nodes(grid.get("num_points", 33), task["ele_T"], grid.get("half_width", 0.4)) + else: + raise ValueError(f"Unknown energy-grid method: {method}") + if conductance or method == "clenshaw_curtis": + if task.get("unit", "Hartree") != "eV": + raise ValueError("Conductance and Clenshaw-Curtis currently require unit='eV'") + if task.get("scf", False): + raise ValueError("Conductance and Clenshaw-Curtis currently require scf=false") + if not np.isfinite(task["ele_T"]) or task["ele_T"] <= 0: + raise ValueError("Conductance requires ele_T > 0") + for lead in ("lead_L", "lead_R"): + voltage = task["stru_options"][lead]["voltage"] + if not np.isfinite(voltage) or abs(voltage) > 1e-12: + raise ValueError("Conductance and Clenshaw-Curtis require zero lead voltages") + opts = task.get("conductance_options") or {} + mu = opts.get("mu", ["Ef"]) + if not isinstance(mu, list) or not mu: + raise ValueError("conductance_options.mu must be a nonempty list") + for center in mu: + if isinstance(center, str): + if center not in ("Ef", "Ev", "Ec"): + raise ValueError("mu labels must be Ef, Ev or Ec") + elif isinstance(center, bool) or not isinstance(center, (int, float)) or not np.isfinite(center): + raise ValueError("numeric mu must be a finite energy in eV") + if opts.get("lead", "lead_L") not in ("lead_L", "lead_R"): + raise ValueError("conductance_options.lead must be lead_L or lead_R") + if opts.get("spin_degeneracy") not in (None, 1, 2): + raise ValueError("spin_degeneracy must be 1 or 2") + if any(center in ("Ev", "Ec") for center in mu): + task["stru_options"]["compute_band_edges"] = True + if method == "clenshaw_curtis": + outputs = task.get("output_options", {}) + if any(outputs.get(key, False) for key in ("density", "potential", "current", "current_nscf", "lcurrent")): + raise ValueError("Clenshaw-Curtis is a transmission grid, not a density or finite-bias current grid") def negf(): doc_scf = "Whether to run a self-consistent (Poisson-NEGF) loop. Default: False (non-self-consistent transport)." @@ -127,9 +202,10 @@ def negf(): Argument("stru_options", dict, optional=False, sub_fields=stru_options(), doc=doc_stru_options), Argument("poisson_options", dict, optional=True, default={}, sub_fields=[], sub_variants=[poisson_options()], doc=doc_poisson_options), Argument("self_energy_options", dict, optional=True, default={}, sub_fields=self_energy_options(), doc=doc_self_energy_options), - Argument("espacing", [int, float], optional=False, doc=doc_espacing), - Argument("emin", [int, float], optional=False, doc=doc_emin), - Argument("emax", [int, float], optional=False, doc=doc_emax), + Argument("energy_grid", dict, optional=True, default={}, sub_fields=[], sub_variants=[energy_grid_options()], + doc="Fixed transmission grad. Legacy emin/emax/espacing migrate to uniform. CC uses eV relative to E_ref."), + Argument("conductance_options", dict, optional=True, default={}, + sub_fields=conductance_options(), doc="Finite-temperature Landauer integration targets"), Argument("rgf_options", dict, optional=True, default={}, sub_fields=rgf_options_group(), doc=doc_rgf_options), Argument("e_fermi", [int, float], optional=True, default=None ,doc=doc_e_fermi), Argument("density_options", dict, optional=True, default={}, sub_fields=[], sub_variants=[density_options()], doc=doc_density_options), @@ -143,6 +219,8 @@ def output_options(): return [ Argument("dos", bool, optional=True, default=False, doc="Write density of states."), Argument("tc", bool, optional=True, default=False, doc="Write transmission coefficient."), + Argument("conductance", bool, optional=True, default=False, + doc="Write zero-bias non-SCF conductance; requires unit=eV and ele_T>0."), Argument("density", bool, optional=True, default=False, doc="Write electron density."), Argument("potential", bool, optional=True, default=False, doc="Write electrostatic potential."), Argument("current", bool, optional=True, default=False, doc="Write self-consistent current."), @@ -490,12 +568,18 @@ def normalize_run(data): data = run_op.normalize_value(data) run_op.check_value(data, strict=True) + if data.get("task_options", {}).get("task") == "negf": + validate_energy_options(data["task_options"]) + return data # (flat legacy key) -> tuple describing where it now lives inside task_options. # Tuples of length 1 = top-level under task_options; length 2 = one sub-dict; etc. _NEGF_LEGACY_KEY_MAP = { + "emin": ("energy_grid", "emin"), + "emax": ("energy_grid", "emax"), + "espacing": ("energy_grid", "espacing"), "sgf_solver": ("self_energy_options", "solver"), "se_numba_jit": ("self_energy_options", "numba_jit"), "se_info_display": ("self_energy_options", "info_display"), @@ -540,6 +624,7 @@ def _migrate_legacy_negf_task_options(data): if task_options.get("task") != "negf": return + explicit_grid = "energy_grid" in task_options for legacy_key, path in _NEGF_LEGACY_KEY_MAP.items(): if legacy_key not in task_options: continue @@ -551,6 +636,13 @@ def _migrate_legacy_negf_task_options(data): f"Migrating value automatically for this run." ) + if explicit_grid and legacy_key in ("emin", "emax", "espacing"): + log.warning( + f"NEGF task_options: both legacy '{legacy_key}' and " + f"'energy_grid' are set; keeping the value from 'energy_grid'." + ) + continue + cursor = task_options for parent in path[:-1]: existing = cursor.get(parent) diff --git a/dpnegf/utils/conductance.py b/dpnegf/utils/conductance.py new file mode 100644 index 0000000..76d8694 --- /dev/null +++ b/dpnegf/utils/conductance.py @@ -0,0 +1,133 @@ +""" +Finite-temperature, zero-bias conductance calculation. +""" + +import numpy as np +from scipy.constants import elementary_charge, h + +from dpnegf.utils.energy_grid import ( + EnergyGrid, build_energy_grid, fermi_window, thermal_energy, trapezoid_weights, +) +from scipy.special import expit + +G0_S = 2 * elementary_charge ** 2 / h + +def _array(value): + """Convert input data to a finite, real-valued float64 NumPy array.""" + if hasattr(value, "detach"): + value = value.detach().cpu().numpy() + value = np.asarray(value) + if np.iscomplexobj(value) or not np.isfinite(value).all(): + raise ValueError("transmission must be real and finite") + return value.astype(np.float64, copy=False) + + +def _result(transmission, weights, spin_degeneracy, metadata): + if spin_degeneracy not in (1, 2): + raise ValueError("spin_degeneracy must be 1 or 2") + values = _array(transmission) + if values.ndim == 0 or values.shape[-1] != len(weights): + raise ValueError("last transmission dimension must match the energy nodes") + integral = values @ np.asarray(weights) + reduced = integral * (spin_degeneracy / 2) + return dict(metadata, G_over_G0=reduced.tolist(), G_S=(G0_S * reduced).tolist(), + G0_S=G0_S, spin_degeneracy=spin_degeneracy) + + +def conductance_uniform(energies, transmission, mu=0.0, temperature=300.0, + spin_degeneracy=2): + """Trapezoidal integral of T(E)*(-df/dE) on a uniform energy grid. + + The integral is restricted to the supplied window; tails are not + extrapolated or renormalized. Refinement is required to assess quadrature + error. The omitted kernel mass alone is not a conductance error bound. + """ + raw_energies = ( + energies.detach().cpu().numpy() if hasattr(energies, "detach") else np.asarray(energies) + ) + precision = ( np.finfo(raw_energies.dtype).eps if np.issubdtype(raw_energies.dtype, np.floating) else 0.0) + + energies = _array(raw_energies) + weights = trapezoid_weights(energies) + delta = np.diff(energies) + + spacing_atol = max(1e-10, 4 * precision * np.max(np.abs(energies)),) + + if not np.allclose( delta, delta.mean(), rtol=1e-4, atol=spacing_atol): + raise ValueError("conductance_uniform requires uniformly spaced energies") + + if not np.isfinite(mu): + raise ValueError("mu must be finite") + + kbt = thermal_energy(temperature) + tail = float( + expit((energies[0] - mu) / kbt) + + expit((mu - energies[-1]) / kbt) + ) + + return _result( transmission, weights * fermi_window(energies, mu, temperature), spin_degeneracy, + { + "method": "uniform", + "mu_eV": float(mu), + "temperature_K": float(temperature), + "energy_unit": "eV", + "window_eV": [ + float(energies[0]), + float(energies[-1]), + ], + "kernel_mass": 1 - tail, + "omitted_kernel_mass": tail, + }, + ) + +def conductance_clenshaw_curtis( + transmission, grid=None, *, mu=0.0, + temperature=300.0, half_width=0.4, + num_points=33,spin_degeneracy=2): + """Integrate transmission evaluated at CC nodes, never interpolated. + + With no grid, transmission must follow the single-target fixed rule + specified by the keyword arguments. Returns one record for a single target, + or a list for a multi-target grid. + """ + if grid is None: + grid = build_energy_grid( + { + "method": "clenshaw_curtis", + "num_points": num_points, + "half_width": half_width, + }, [mu], temperature, + ) + + if grid.method != "clenshaw_curtis": + raise ValueError("conductance_clenshaw_curtis requires a CC energy grid") + + results = integrate_conductance(transmission,grid,spin_degeneracy) + + return results[0] if len(results) == 1 else results + + + +def integrate_conductance(transmission,grid,spin_degeneracy=2): + """Integrate all targets on an EnergyGrid without resampling transmission.""" + if not isinstance(grid, EnergyGrid): + raise TypeError("grid must be an EnergyGrid returned by build_energy_grid") + values = _array(transmission) + + if values.ndim == 0 or values.shape[-1] != len(grid.energies): + raise ValueError("transmission must match the combined energy grid") + + return [ + _result( + values[..., target["indices"]], + target["weights"], + spin_degeneracy, + dict( + target, + method=grid.method, + temperature_K=grid.temperature, + energy_unit="eV", + ), + ) + for target in grid.targets + ] \ No newline at end of file diff --git a/dpnegf/utils/energy_grid.py b/dpnegf/utils/energy_grid.py new file mode 100644 index 0000000..15106ac --- /dev/null +++ b/dpnegf/utils/energy_grid.py @@ -0,0 +1,255 @@ +""" +Fixed transmission grids and Landauer quadrature weights (energies in eV) +""" + +from dataclasses import dataclass + +import numpy as np +import torch +from scipy.special import expit + +from dpnegf.utils.constants import Boltzmann, eV2J + + + +def thermal_energy(temperature): + if not np.isfinite(temperature) or temperature <= 0: + raise ValueError( "temperature must be finite and greater than zero Kelvin") + return Boltzmann * float(temperature) / eV2J + + +def fermi_window(energies, mu, temperature): + kbt = thermal_energy(temperature) + # expit(-abs(x)) avoids cancellation on the positive-energy tail. + p = expit(-np.abs((np.asarray(energies) - mu) / kbt)) + + return p * (1 - p) / kbt + + +def trapezoid_weights(energies): + energies = np.asarray(energies, dtype=float) + + if ( + energies.ndim != 1 + or len(energies) < 2 + or not np.isfinite(energies).all() + or np.any(np.diff(energies) <= 0) + ): + raise ValueError("energies must be finite, strictly increasing, and contain at least two points") + + delta = np.diff(energies) + return np.r_[ + delta[0], + delta[:-1] + delta[1:], + delta[-1], + ] / 2 + + +def uniform_nodes(options): + lo, hi = options["emin"], options["emax"] + + if not np.isfinite([lo, hi]).all() or hi <= lo: + raise ValueError("uniform energy grid requires finite emin < emax") + + n, spacing = options.get("num_points"), options.get("espacing") + + if (n is None) == (spacing is None): + raise ValueError( "uniform grid requires exactly one of num_points and espacing" ) + + if spacing is not None: + if not np.isfinite(spacing) or spacing <= 0: + raise ValueError("espacing must be finite and positive") + + n = int((hi - lo) / spacing) + + if n < 2: + raise ValueError("espacing must produce at least two energy points") + + # Preserve both the node count and torch default dtype of legacy inputs. + return torch.linspace(lo, hi, steps=n).numpy() + + if (isinstance(n, bool) or not isinstance(n, (int, np.integer)) or n < 2): + raise ValueError("num_points must be an integer at least two") + + return np.linspace(lo, hi, n, dtype=np.float64) + + +def clenshaw_curtis_nodes( + num_points=33, + temperature=300.0, + half_width=0.4, +): + """Return finite energy offsets and probability-space integration weights. + + The finite interval is [-half_width, half_width]. Weights sum to its + Fermi-window mass, not one. No transmission extrapolation is performed. + Odd orders include zero; orders 2**m + 1 are nested. + """ + n = num_points + + if ( + isinstance(n, bool) + or not isinstance(n, (int, np.integer)) + or n < 3 + or n % 2 == 0 + ): + raise ValueError( + "Clenshaw-Curtis requires an odd num_points >= 3" + ) + + kbt = thermal_energy(temperature) + + if not np.isfinite(half_width) or half_width <= 0: + raise ValueError("half_width must be finite and positive") + + order = n - 1 + theta = np.pi * np.arange(n) / order + + weights = np.empty(n) + weights[0] = weights[-1] = 1 / (order * order - 1) + + interior = np.ones(order - 1) + + for k in range(1, order // 2): + interior -= (2 * np.cos(2 * k * theta[1:-1])/ (4 * k * k - 1)) + + interior -= (np.cos(order * theta[1:-1])/ (order * order - 1)) + weights[1:-1] = 2 * interior / order + p0 = expit(-half_width / kbt) + mass = np.tanh(half_width / (2 * kbt)) + + # Compute the negative half and reflect to avoid logit(1) at low T. + p = (p0 + mass * np.sin(theta[: n // 2 + 1] / 2)**2) + + offsets = np.empty(n) + offsets[0] = -half_width + offsets[1 : n // 2] = kbt * (np.log(p[1:-1]) - np.log1p(-p[1:-1])) + offsets[n // 2] = 0.0 + offsets[n // 2 + 1 :] = -offsets[: n // 2][::-1] + + return offsets, weights * (mass / 2) + + +@dataclass +class EnergyGrid: + """Sorted unique nodes and per-target indices/weights into that array.""" + + energies: np.ndarray + method: str + temperature: float + targets: list + + def metadata(self): + return { + "method": self.method, + "temperature_K": self.temperature, + "energy_unit": "eV", + "targets": self.targets, + } + + +def _uniform(options, mu, temperature): + energies = uniform_nodes(options) + + return ( + energies, + trapezoid_weights(energies) + * fermi_window(energies, mu, temperature), + ) + + +def _clenshaw_curtis(options, mu, temperature): + offsets, weights = clenshaw_curtis_nodes( + options.get("num_points", 33), + temperature, + options.get("half_width", 0.4), + ) + + return mu + offsets, weights + + +GRID_METHODS = { + "uniform": _uniform, + "clenshaw_curtis": _clenshaw_curtis, +} + + +def build_energy_grid( + options, + mu=(0.0,), + temperature=300.0, + labels=None, +): + """Build fixed nodes for all chemical potentials, relative to one E_ref. + + mu is a sequence of numeric energies in eV. Band-edge labels must be + resolved by the caller. Repeated targets are retained but evaluated once. + """ + method = options.get("method", "uniform") + + if method not in GRID_METHODS: + raise ValueError(f"Unknown energy-grid method: {method}") + + thermal_energy(temperature) + mu = np.asarray(mu, dtype=float) + + if ( + mu.ndim != 1 + or len(mu) == 0 + or not np.isfinite(mu).all() + ): + raise ValueError( + "mu must be a nonempty sequence of finite energies" + ) + + labels = list(mu) if labels is None else list(labels) + + if len(labels) != len(mu): + raise ValueError("labels and mu must have equal lengths") + + rules = [ + GRID_METHODS[method](options, center, temperature) + for center in mu + ] + + energies = np.unique( + np.concatenate([e for e, w in rules]) + ) + + # Existing SE cache uses E_{energy:.8f}; do not silently merge collisions. + if len({f"{e:.8f}" for e in energies}) != len(energies): + raise ValueError( + "Distinct energy nodes collide at the self-energy " + "cache precision (8 decimals)" + ) + + targets = [] + kbt = thermal_energy(temperature) + + for center, label, (nodes, weights) in zip(mu, labels, rules): + tail = float( + expit((nodes[0] - center) / kbt) + + expit((center - nodes[-1]) / kbt) + ) + + targets.append( + { + "label": label, + "mu_eV": float(center), + "indices": np.searchsorted(energies, nodes).tolist(), + "weights": weights.tolist(), + "window_eV": [ + float(nodes[0]), + float(nodes[-1]), + ], + "kernel_mass": 1 - tail, + "omitted_kernel_mass": tail, + } + ) + + return EnergyGrid( + energies, + method, + float(temperature), + targets, + ) \ No newline at end of file diff --git a/examples/graphene/negf.json b/examples/graphene/negf.json index fe53d5d..a58d8d6 100644 --- a/examples/graphene/negf.json +++ b/examples/graphene/negf.json @@ -1,14 +1,11 @@ { - "structure": "./struct_negf.xyz", + "structure": "./stru_negf.xyz", "task_options": { "task": "negf", "scf": false, "block_tridiagonal": true, "ele_T": 300, "unit": "eV", - "espacing": 0.05, - "emin": -10, - "emax": 10, "eta_lead": 1e-05, "eta_device": 0.0, "stru_options": { @@ -61,13 +58,27 @@ "err": 1e-05 }, "self_energy_options": { - "solver": "Sancho-Rubio" + "solver": "Sancho-Rubio", + "numba_jit": true }, "output_options": { - "dos": true, "tc": true, - "ldos": true, - "current_nscf": false + "conductance": true + }, + "energy_grid": { + "method": "clenshaw_curtis", + "num_points": 33, + "half_width": 0.4 + }, + "conductance_options": { + "mu": [ + "Ef" + ], + "lead": "lead_L" + }, + "rgf_options": { + "device": "cpu", + "e_batch_size": 8 } } } diff --git a/examples/hBN/negf.json b/examples/hBN/negf.json index d567863..c76e71f 100644 --- a/examples/hBN/negf.json +++ b/examples/hBN/negf.json @@ -69,7 +69,8 @@ "dos": true, "tc": true, "ldos": true, - "current_nscf": false + "current_nscf": false, + "conductance": true } } } diff --git a/examples/hBN/negf_conductance_cc.json b/examples/hBN/negf_conductance_cc.json new file mode 100644 index 0000000..ebf76da --- /dev/null +++ b/examples/hBN/negf_conductance_cc.json @@ -0,0 +1,95 @@ +{ + "structure": "./stru_negf.xyz", + "task_options": { + "task": "negf", + "scf": false, + "block_tridiagonal": true, + "ele_T": 300, + "unit": "eV", + "eta_lead": 1e-05, + "eta_device": 0.0, + "stru_options": { + "gamma_center": false, + "time_reversal_symmetry": true, + "nel_atom": { + "B": 3, + "N": 5 + }, + "kmesh": [ + 1, + 1, + 1 + ], + "pbc": [ + false, + true, + false + ], + "device": { + "id": "32-64", + "sort": true + }, + "compute_band_edges": true, + "lead_L": { + "id": "0-32", + "voltage": 0.0, + "kmesh_lead_Ef": [ + 1, + 50, + 20 + ], + "useBloch": false + }, + "lead_R": { + "id": "64-96", + "voltage": 0.0, + "kmesh_lead_Ef": [ + 1, + 50, + 20 + ], + "useBloch": false + } + }, + "density_options": { + "method": "Fiori", + "integrate_way": "direct" + }, + "poisson_options": { + "solver": "fmm", + "err": 1e-05 + }, + "self_energy_options": { + "solver": "Sancho-Rubio", + "numba_jit": true, + "parallel": { + "cpu_budget": 32, + "n_workers": 1, + "blas_threads": 1 + } + }, + "output_options": { + "dos": true, + "tc": true, + "ldos": true, + "current_nscf": false, + "conductance": true + }, + "energy_grid":{ + "method": "clenshaw_curtis", + "num_points": 33, + "half_width": 0.4 + }, + "conductance_options": { + "mu": [ + "Ev", + "Ec" + ], + "lead": "lead_L" + }, + "rgf_options": { + "device": "cpu", + "e_batch_size": 8 + } + } +} diff --git a/examples/hBN/negf_conductance_uniform.json b/examples/hBN/negf_conductance_uniform.json new file mode 100644 index 0000000..1e2cb10 --- /dev/null +++ b/examples/hBN/negf_conductance_uniform.json @@ -0,0 +1,96 @@ +{ + "structure": "./stru_negf.xyz", + "task_options": { + "task": "negf", + "scf": false, + "block_tridiagonal": true, + "ele_T": 300, + "unit": "eV", + "eta_lead": 1e-05, + "eta_device": 0.0, + "stru_options": { + "gamma_center": false, + "time_reversal_symmetry": true, + "nel_atom": { + "B": 3, + "N": 5 + }, + "kmesh": [ + 1, + 1, + 1 + ], + "pbc": [ + false, + true, + false + ], + "device": { + "id": "32-64", + "sort": true + }, + "compute_band_edges": true, + "lead_L": { + "id": "0-32", + "voltage": 0.0, + "kmesh_lead_Ef": [ + 1, + 50, + 20 + ], + "useBloch": false + }, + "lead_R": { + "id": "64-96", + "voltage": 0.0, + "kmesh_lead_Ef": [ + 1, + 50, + 20 + ], + "useBloch": false + } + }, + "density_options": { + "method": "Fiori", + "integrate_way": "direct" + }, + "poisson_options": { + "solver": "fmm", + "err": 1e-05 + }, + "self_energy_options": { + "solver": "Sancho-Rubio", + "numba_jit": true, + "parallel": { + "cpu_budget": 32, + "n_workers": 1, + "blas_threads": 1 + } + }, + "output_options": { + "dos": true, + "tc": true, + "ldos": true, + "current_nscf": false, + "conductance": true + }, + "energy_grid":{ + "method": "uniform", + "emin": -10, + "emax": 10, + "num_points": 4001 + }, + "conductance_options": { + "mu": [ + "Ev", + "Ec" + ], + "lead": "lead_L" + }, + "rgf_options": { + "device": "cpu", + "e_batch_size": 8 + } + } +}