Implement current step

examples/coulomb_blockade/pentacene/curr.py

@@ -0,0 +1,133 @@
+#!/usr/bin/env python
+
+import pickle
+from argparse import ArgumentParser
+from pathlib import Path
+
+import numpy as np
+from ase.units import _e, _hplanck, kB
+
+G0 = 2.0 * _e**2 / _hplanck
+
+
+def natural_sort(filepath: Path) -> float:
+    """docstring"""
+    numerical_part = filepath.stem.split("_")[-1]
+    return float(numerical_part)
+
+
+def fermidistribution(energy, kt):
+    # fermi level is fixed to zero
+    # energy can be a single number or a list
+    assert kt >= 0.0, "Negative temperature encountered!"
+
+    if kt == 0:
+        if isinstance(energy, float):
+            return int(energy / 2.0 <= 0)
+        else:
+            return (energy / 2.0 <= 0).astype(int)
+    else:
+        return 1.0 / (1.0 + np.exp(energy / kt))
+
+
+def get_current(bias, E, transmission, T=300, unit="uA"):
+    """Get the current in nA.
+
+    Fermi-Dirac distribution = 1 / (1 + e^((E - mu) / kT))
+    """
+
+    # TODO refactor redundant computations
+
+    if not isinstance(bias, (int, float)):
+        bias = bias[np.newaxis]
+        E = E[:, np.newaxis]
+        transmission = transmission[:, np.newaxis]
+
+    mu = bias / 2.0
+    kT = kB * T
+
+    fl = fermidistribution(E - mu, kT)
+    fr = fermidistribution(E + mu, kT)
+
+    return G0 * np.trapz((fl - fr) * transmission, x=E, axis=0) * 1e6  # uA
+
+
+def numerical_derivative(x, y):
+    """docstring"""
+    dy_dx = np.diff(y) / np.diff(x)
+    dy_dx = np.append(dy_dx, np.nan)
+
+    return dy_dx
+
+
+def compute_current(
+    energies: np.ndarray,
+    V_min=-2.5,
+    V_max=2.5,
+    dV=0.1,
+    temperature=9,
+    transmission_folder_name="transmission_folder",
+) -> None:
+    """docstring"""
+
+    files = []
+    for file in sorted(
+        Path(transmission_folder_name).iterdir(),
+        key=natural_sort,
+    ):
+        files.append(file)
+
+    bias = np.linspace(V_min, V_max, int((V_max - V_min) / dV) + 1)
+    transmissions = np.asarray([np.load(fn) for fn in files])
+
+    energies = energies.real
+
+    current = np.asarray(
+        [
+            get_current(bias, energies, transmission, temperature)
+            for transmission in transmissions
+        ]
+    )
+
+    derivative = np.asarray([numerical_derivative(bias, i) for i in current])
+
+    np.save("current.npy", current)
+    np.save("derivative.npy", derivative)
+
+
+if __name__ == "__main__":
+    """docstring"""
+
+    parser = ArgumentParser()
+
+    parser.add_argument(
+        "-pf",
+        "--parameters-filename",
+        help="name of parameters file",
+    )
+
+    parser.add_argument(
+        "-ef",
+        "--energies-filename",
+        help="name of energies file",
+    )
+
+    parser.add_argument(
+        "-tf",
+        "--transmission-folder-name",
+        help="name of folder containing transmission files",
+    )
+
+    args = parser.parse_args()
+
+    with open(args.parameters_filename, "rb") as file:
+        parameters = pickle.load(file)
+
+    with open(args.energies_filename, "rb") as file:
+        energies = np.load(file)
+
+    compute_current(
+        energies,
+        transmission_folder_name=args.transmission_folder_name,
+        **parameters,
+    )

pyproject.toml

@@ -34,6 +34,7 @@ requires-python = ">=3.9"
 "quantum_transport.hybridize" = "aiida_quantum_transport.calculations.hybridize:HybridizationCalculation"
 "quantum_transport.dmft" = "aiida_quantum_transport.calculations.dmft:DMFTCalculation"
 "quantum_transport.transmission" = "aiida_quantum_transport.calculations.transmission:TransmissionCalculation"
+"quantum_transport.current" = "aiida_quantum_transport.calculations.current:CurrentCalculation"
 
 [project.entry-points."aiida.cmdline.data"]
 "quantum_transport" = "aiida_quantum_transport.cli.commands:data_cli"
@@ -48,6 +49,7 @@ requires-python = ">=3.9"
 "quantum_transport.hybridize" = "aiida_quantum_transport.parsers.hybridize:HybridizationParser"
 "quantum_transport.dmft" = "aiida_quantum_transport.parsers.dmft:DMFTParser"
 "quantum_transport.transmission" = "aiida_quantum_transport.parsers.transmission:TransmissionParser"
+"quantum_transport.current" = "aiida_quantum_transport.parsers.current:CurrentParser"
 
 [project.optional-dependencies]
 docs = ["sphinx-design~=0.4.1", "pydata-sphinx-theme==0.13.3"]

src/aiida_quantum_transport/calculations/__init__.py

@@ -1,3 +1,4 @@
+from .current import CurrentCalculation
 from .custom import CustomCalculation
 from .dft import DFTCalculation
 from .dmft import DMFTCalculation
@@ -14,4 +15,5 @@
     "HybridizationCalculation",
     "DMFTCalculation",
     "TransmissionCalculation",
+    "CurrentCalculation",
 ]

src/aiida_quantum_transport/calculations/current.py

@@ -0,0 +1,116 @@
+from __future__ import annotations
+
+import pickle
+from typing import TYPE_CHECKING
+
+from aiida import orm
+from aiida.common.datastructures import CalcInfo, CodeInfo
+from aiida.common.folders import Folder
+from aiida.engine import CalcJob
+
+if TYPE_CHECKING:
+    from aiida.engine.processes.calcjobs.calcjob import CalcJobProcessSpec
+
+
+class CurrentCalculation(CalcJob):
+    """docstring"""
+
+    _default_parser_name = "quantum_transport.current"
+
+    @classmethod
+    def define(cls, spec: CalcJobProcessSpec) -> None:
+        """docstring"""
+
+        super().define(spec)
+
+        spec.input(
+            "code",
+            valid_type=orm.AbstractCode,
+            help="The current script",
+        )
+
+        spec.input(
+            "parameters",
+            valid_type=orm.Dict,
+            required=False,
+            default=lambda: orm.Dict({}),
+            help="parameters used to compute current",
+        )
+
+        spec.input(
+            "hybridization.energies_file",
+            valid_type=orm.SinglefileData,
+            help="file containing computed energies",
+        )
+
+        spec.input(
+            "transmission.transmission_folder",
+            valid_type=orm.FolderData,
+            help="folder containing transmission files",
+        )
+
+        spec.input(
+            "metadata.options.parser_name",
+            valid_type=str,
+            default=cls._default_parser_name,
+        )
+
+        spec.output(
+            "current_file",
+            valid_type=orm.SinglefileData,
+            help="The current data file",
+        )
+
+        spec.output(
+            "derivative_file",
+            valid_type=orm.SinglefileData,
+            help="The current derivative data file",
+        )
+
+        spec.exit_code(
+            400,
+            "ERROR_ACCESSING_OUTPUT_FILE",
+            "an issue occurred while accessing an expected retrieved file",
+        )
+
+    def prepare_for_submission(self, folder: Folder) -> CalcInfo:
+        """docstring"""
+
+        pickled_parameters_filename = "parameters.pkl"
+        with folder.open(pickled_parameters_filename, "wb") as file:
+            parameters: orm.Dict = self.inputs.parameters
+            pickle.dump(parameters.get_dict(), file)
+
+        codeinfo = CodeInfo()
+        codeinfo.code_uuid = self.inputs.code.uuid
+        codeinfo.cmdline_params = [
+            "--parameters-filename",
+            pickled_parameters_filename,
+            "--energies-filename",
+            self.inputs.hybridization.energies_file.filename,
+            "--transmission-folder-name",
+            "transmission_folder",
+        ]
+
+        self.node.get_remote_workdir()
+
+        calcinfo = CalcInfo()
+        calcinfo.codes_info = [codeinfo]
+        calcinfo.local_copy_list = [
+            (
+                self.inputs.hybridization.energies_file.uuid,
+                self.inputs.hybridization.energies_file.filename,
+                self.inputs.hybridization.energies_file.filename,
+            ),
+            (
+                self.inputs.transmission.transmission_folder.uuid,
+                ".",
+                "transmission_folder",
+            ),
+        ]
+        calcinfo.retrieve_list = [
+            "current.npy",
+            "derivative.npy",
+        ]
+
+        return calcinfo

src/aiida_quantum_transport/parsers/__init__.py

@@ -1,3 +1,4 @@
+from .current import CurrentParser
 from .custom import CustomParser
 from .dft import DFTParser
 from .dmft import DMFTParser
@@ -12,4 +13,5 @@
     "HybridizationParser",
     "DMFTParser",
     "TransmissionParser",
+    "CurrentParser",
 ]

src/aiida_quantum_transport/parsers/current.py

@@ -0,0 +1,29 @@
+from __future__ import annotations
+
+from pathlib import Path
+
+from aiida import orm
+from aiida.engine import ExitCode
+from aiida.parsers import Parser
+
+
+class CurrentParser(Parser):
+    """docstring"""
+
+    _OUTPUT_FILE_MAPPING = {
+        "current": "current.npy",
+        "derivative": "derivative.npy",
+    }
+
+    def parse(self, **kwargs) -> ExitCode | None:
+        """docstring"""
+
+        try:
+            with self.retrieved.as_path() as retrieved_path:
+                for label, filename in self._OUTPUT_FILE_MAPPING.items():
+                    path = Path(retrieved_path) / filename
+                    self.out(f"{label}_file", orm.SinglefileData(path))
+        except OSError:
+            return self.exit_codes.ERROR_ACCESSING_OUTPUT_FILE
+
+        return None

src/aiida_quantum_transport/workchains/coulomb_diamonds.py

@@ -6,6 +6,7 @@
 from aiida.engine import ToContext, WorkChain
 
 from aiida_quantum_transport.calculations import (
+    CurrentCalculation,
     DFTCalculation,
     DMFTCalculation,
     HybridizationCalculation,
@@ -110,6 +111,12 @@ def define(cls, spec: WorkChainSpec) -> None:
             include=["code", "metadata"],
         )
 
+        spec.expose_inputs(
+            CurrentCalculation,
+            namespace="current",
+            include=["code", "parameters", "metadata"],
+        )
+
         spec.expose_outputs(
             DFTCalculation,
             namespace="dft.leads",
@@ -145,6 +152,11 @@ def define(cls, spec: WorkChainSpec) -> None:
             namespace="transmission",
         )
 
+        spec.expose_outputs(
+            CurrentCalculation,
+            namespace="current",
+        )
+
         spec.outline(
             cls.run_dft,
             cls.define_scattering_region,
@@ -153,7 +165,7 @@ def define(cls, spec: WorkChainSpec) -> None:
             cls.run_dmft_converge_mu,
             cls.run_dmft_sweep_mu,
             cls.compute_transmission,
-            # cls.compute_current,
+            cls.compute_current,
             cls.gather_results,
         )
 
@@ -333,6 +345,25 @@ def compute_transmission(self):
 
     def compute_current(self):
         """docstring"""
+        current_inputs = {
+            "parameters": self.ctx.hybridization.inputs.parameters,
+            "hybridization": {
+                "energies_file": self.ctx.hybridization.outputs.energies_file,
+            },
+            "transmission": {
+                "transmission_folder": self.ctx.transmission.outputs.transmission_folder,
+            },
+            **self.exposed_inputs(
+                CurrentCalculation,
+                namespace="current",
+            ),
+        }
+        return ToContext(
+            current=self.submit(
+                CurrentCalculation,
+                **current_inputs,
+            )
+        )
 
     def gather_results(self):
         """docstring"""
@@ -392,3 +423,11 @@ def gather_results(self):
                 namespace="transmission",
             )
         )
+
+        self.out_many(
+            self.exposed_outputs(
+                self.ctx.current,
+                CurrentCalculation,
+                namespace="current",
+            )
+        )