diff --git a/qmp/algorithms/chop_imag.py b/qmp/algorithms/chop_imag.py
index 5059e56a..4c7eb5a4 100644
--- a/qmp/algorithms/chop_imag.py
+++ b/qmp/algorithms/chop_imag.py
@@ -6,11 +6,20 @@
 import typing
 import dataclasses
 import omegaconf
+import torch
 import torch.utils.tensorboard
 from ..utility.context import RuntimeContext
 from ..utility.action_dict import action_dict
 
 
+def _get_devices_from_context(context: RuntimeContext) -> list[torch.device]:
+    """
+    Get the devices list from the runtime context.
+    """
+    assert context.devices is not None
+    return context.devices
+
+
 @dataclasses.dataclass
 class ChopImagConfig:
     """
@@ -35,6 +44,9 @@ def main(
         model = context.create_model(runtime_config.model)
         data = checkpoint_data
 
+        # Get devices for multi-GPU computation
+        devices = _get_devices_from_context(context)
+
         logging.info(
             "Arguments Summary: Chop Size: %d, Second Order Magnitude: %.10f",
             self.chop_size,
@@ -59,7 +71,7 @@ def main(
             logging.info("The number of configurations: %d", num_configs)
             writer.add_scalar("chop_imag/num_configs", num_configs, i)  # type: ignore[no-untyped-call]
             psi = psi / psi.norm()
-            hamiltonian_psi = model.apply_within(configs, psi, configs)
+            hamiltonian_psi = model.apply_within(configs, psi, configs, devices)
             psi_hamiltonian_psi = (psi.conj() @ hamiltonian_psi).real
             energy = psi_hamiltonian_psi
             logging.info(
diff --git a/qmp/algorithms/guide.py b/qmp/algorithms/guide.py
index df3e75f6..651774b6 100644
--- a/qmp/algorithms/guide.py
+++ b/qmp/algorithms/guide.py
@@ -12,6 +12,14 @@
 from ..utility.action_dict import action_dict
 
 
+def _get_devices_from_context(context: RuntimeContext) -> list[torch.device]:
+    """
+    Get the devices list from the runtime context.
+    """
+    assert context.devices is not None
+    return context.devices
+
+
 @dataclasses.dataclass
 class GuideConfig:
     """
@@ -50,6 +58,9 @@ def main(
 
         data = checkpoint_data
 
+        # Get devices for multi-GPU computation
+        devices = _get_devices_from_context(context)
+
         # Create Optimizers
         # We use the same optimizer configuration for both network and sampling
         optimizer_network = context.create_optimizer(
@@ -96,7 +107,7 @@ def main(
                     [
                         configs_src_network,
                         model.find_relative(
-                            configs_src_network, psi_src_network, self.relative_count - len(configs_src_network)
+                            configs_src_network, psi_src_network, self.relative_count - len(configs_src_network), devices=devices
                         ),
                     ]
                 )
@@ -107,7 +118,7 @@ def main(
                     [
                         configs_src_sampling,
                         model.find_relative(
-                            configs_src_sampling, psi_src_sampling, self.relative_count - len(configs_src_sampling)
+                            configs_src_sampling, psi_src_sampling, self.relative_count - len(configs_src_sampling), devices=devices
                         ),
                     ]
                 )
@@ -118,7 +129,7 @@ def energy_sampling() -> torch.Tensor:
                 psi_src = network(configs_src)
                 with torch.no_grad():
                     psi_dst = network(configs_dst)
-                    hamiltonian_psi_dst = model.apply_within(configs_dst, psi_dst, configs_src)
+                    hamiltonian_psi_dst = model.apply_within(configs_dst, psi_dst, configs_src, devices)
                 psi_src_reweight = psi_src.conj() * reweight_sampling
                 num = psi_src_reweight @ hamiltonian_psi_dst
                 den = psi_src_reweight @ psi_src.detach()
@@ -132,7 +143,7 @@ def energy_network() -> torch.Tensor:
                     configs_dst = configs_dst_network
                     psi_src = network(configs_src)
                     psi_dst = network(configs_dst)
-                    hamiltonian_psi_dst = model.apply_within(configs_dst, psi_dst, configs_src)
+                    hamiltonian_psi_dst = model.apply_within(configs_dst, psi_dst, configs_src, devices)
                     psi_src_reweight = psi_src.conj() * reweight_network
                     num = psi_src_reweight @ hamiltonian_psi_dst
                     den = psi_src_reweight @ psi_src.detach()
@@ -145,7 +156,7 @@ def distribution() -> torch.Tensor:
                 with torch.no_grad():
                     psi_src = network(configs_src)
                     psi_dst = network(configs_dst)
-                    hamiltonian_psi_dst = model.apply_within(configs_dst, psi_dst, configs_src)
+                    hamiltonian_psi_dst = model.apply_within(configs_dst, psi_dst, configs_src, devices)
                     psi_src_reweight = psi_src.conj() * reweight_sampling
                     num = psi_src_reweight @ hamiltonian_psi_dst
                     den = psi_src_reweight @ psi_src
diff --git a/qmp/algorithms/haar.py b/qmp/algorithms/haar.py
index 31e78339..f8168723 100644
--- a/qmp/algorithms/haar.py
+++ b/qmp/algorithms/haar.py
@@ -18,6 +18,14 @@
 from ..utility.optimizer import scale_learning_rate
 
 
+def _get_devices_from_context(context: RuntimeContext) -> list[torch.device]:
+    """
+    Get the devices list from the runtime context.
+    """
+    assert context.devices is not None
+    return context.devices
+
+
 @dataclasses.dataclass
 class _DynamicLanczos:
     """
@@ -35,6 +43,7 @@ class _DynamicLanczos:
     first_extend: bool
     eigen_count: int = 1
     period: int = 256
+    devices: list[torch.device] | None = None
 
     def _extend(self, psi: torch.Tensor, basic_configs: torch.Tensor | None = None) -> None:
         if basic_configs is None:
@@ -45,7 +54,7 @@ def _extend(self, psi: torch.Tensor, basic_configs: torch.Tensor | None = None)
         logging.info("Number of core configurations: %d", count_core)
 
         self.configs = torch.cat(
-            [self.configs, self.model.find_relative(basic_configs, psi, self.count_extend, self.configs)]
+            [self.configs, self.model.find_relative(basic_configs, psi, self.count_extend, self.configs, self.devices)]
         )
         count_selected = len(self.configs)
         self.psi = torch.nn.functional.pad(self.psi, (0, count_selected - count_core))
@@ -80,7 +89,7 @@ def package(
                 selected = (psi.conj() * psi).real.argsort(descending=True)[: self.count_extend]
                 configs = self.configs
                 self._extend(psi[selected], self.configs[selected])
-                psi = self.model.apply_within(configs, psi, self.configs)
+                psi = self._apply_within(configs, psi, self.configs)
             for _, [alpha, beta, v] in zip(range(1 + self.step), self._run()):
                 yield package(self._eigh_tridiagonal(alpha, beta, v))
         elif self.single_extend:
@@ -181,6 +190,7 @@ def _apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j:
                 configs_i[start_index:end_index],
                 psi_i[start_index:end_index],
                 configs_j,
+                self.devices,
             )
             if result is None:
                 result = local_result
@@ -429,6 +439,7 @@ def main(
             target_energy: torch.Tensor
             lanczos_results: list[tuple[torch.Tensor, torch.Tensor, torch.Tensor]]
 
+            devices = _get_devices_from_context(context)
             for lanczos_results in _DynamicLanczos(
                 model=model,
                 configs=configs,
@@ -441,6 +452,7 @@ def main(
                 first_extend=self.krylov_extend_first,
                 eigen_count=self.krylov_eigen_count,
                 period=self.krylov_period,
+                devices=devices,
             ).run():
                 target_energy, configs, original_psi = lanczos_results[0]
                 logging.info(
@@ -546,7 +558,7 @@ def closure() -> torch.Tensor:
 
             loss = typing.cast(torch.Tensor, torch.enable_grad(closure)())  # type: ignore[no-untyped-call,call-arg]
             psi: torch.Tensor = loss.psi  # type: ignore[attr-defined]
-            final_energy = ((psi.conj() @ model.apply_within(configs, psi, configs)) / (psi.conj() @ psi)).real
+            final_energy = ((psi.conj() @ model.apply_within(configs, psi, configs, devices)) / (psi.conj() @ psi)).real
             logging.info(
                 "Loss during local optimization: %.10f, Final energy: %.10f, Target energy: %.10f, Reference energy: %.10f, Final error: %.10f",
                 loss.item(),
diff --git a/qmp/algorithms/haar_with_orbit.py b/qmp/algorithms/haar_with_orbit.py
index 7206745f..44746d6d 100644
--- a/qmp/algorithms/haar_with_orbit.py
+++ b/qmp/algorithms/haar_with_orbit.py
@@ -23,6 +23,7 @@
     _DynamicLanczos,
     _sampling_from_last_iteration,
     _merge_pool_from_neural_network_and_pool_from_last_iteration,
+    _get_devices_from_context,
 )
 from .orbit import NaturalOrbitCalculator, _read_fcidump_tensors
 from ..models.optimized_basis import Model as OptimizedModel, ModelConfig as OptimizedModelConfig
@@ -281,6 +282,8 @@ def main(
             target_energy: torch.Tensor
             lanczos_results: list[tuple[torch.Tensor, torch.Tensor, torch.Tensor]]
 
+            # Get devices for multi-GPU computation
+            devices = _get_devices_from_context(context)
             for lanczos_results in _DynamicLanczos(
                 model=model,
                 configs=configs,
@@ -293,6 +296,7 @@ def main(
                 first_extend=self.krylov_extend_first,
                 eigen_count=self.krylov_eigen_count,
                 period=self.krylov_period,
+                devices=devices,
             ).run():
                 target_energy, configs, original_psi = lanczos_results[0]
                 logging.info("Current energy: %.10f, samples: %d", target_energy.item(), len(configs))
@@ -393,7 +397,7 @@ def closure() -> torch.Tensor:
 
             loss = typing.cast(torch.Tensor, torch.enable_grad(closure)())  # type: ignore[no-untyped-call,call-arg]
             psi: torch.Tensor = loss.psi  # type: ignore[attr-defined]
-            final_energy = ((psi.conj() @ model.apply_within(configs, psi, configs)) / (psi.conj() @ psi)).real
+            final_energy = ((psi.conj() @ model.apply_within(configs, psi, configs, devices)) / (psi.conj() @ psi)).real
             logging.info(
                 "Loss: %.10f, Final energy: %.10f, Target energy: %.10f, Reference energy: %.10f, Final error: %.10f",
                 loss.item(),
diff --git a/qmp/algorithms/pert.py b/qmp/algorithms/pert.py
index a48c04a0..8bd54e85 100644
--- a/qmp/algorithms/pert.py
+++ b/qmp/algorithms/pert.py
@@ -6,10 +6,19 @@
 import typing
 import dataclasses
 import omegaconf
+import torch
 from ..utility.context import RuntimeContext
 from ..utility.action_dict import action_dict
 
 
+def _get_devices_from_context(context: RuntimeContext) -> list[torch.device]:
+    """
+    Get the devices list from the runtime context.
+    """
+    assert context.devices is not None
+    return context.devices
+
+
 @dataclasses.dataclass
 class PerturbationConfig:
     """
@@ -29,13 +38,16 @@ def main(
         model = context.create_model(runtime_config.model)
         data = checkpoint_data
 
+        # Get devices for multi-GPU computation
+        devices = _get_devices_from_context(context)
+
         if "haar" not in data and "imag" in data:
             data["haar"] = data.pop("imag")
         configs, psi = data["haar"]["pool"]
         configs = configs.to(context.device)
         psi = psi.to(context.device)
 
-        energy0_num = psi.conj() @ model.apply_within(configs, psi, configs)
+        energy0_num = psi.conj() @ model.apply_within(configs, psi, configs, devices)
         energy0_den = psi.conj() @ psi
         energy0 = (energy0_num / energy0_den).real.item()
         logging.info("Current energy is %.8f", energy0)
@@ -46,7 +58,7 @@ def main(
         current_target_number = number
         logging.info("Starting finding relative configurations with %d.", number)
         while True:
-            other_configs = model.find_relative(configs, psi, current_target_number, configs)
+            other_configs = model.find_relative(configs, psi, current_target_number, configs, devices)
             current_result_number = other_configs.size(0)
             logging.info("Found %d relative configurations.", current_result_number)
             if current_result_number == last_result_number:
@@ -56,9 +68,9 @@ def main(
             logging.info("Doubling target number to %d.", current_target_number)
             break
 
-        hamiltonian_psi = model.apply_within(configs, psi, other_configs)
+        hamiltonian_psi = model.apply_within(configs, psi, other_configs, devices)
         energy2_num = (hamiltonian_psi.conj() * hamiltonian_psi).real / (psi.conj() @ psi).real
-        energy2_den = energy0 - model.diagonal_term(other_configs).real
+        energy2_den = energy0 - model.diagonal_term(other_configs, devices).real
         energy2 = (energy2_num / energy2_den).sum().item()
         logging.info("Correct energy is %.8f", energy2)
         logging.info("Error is reduced from %.8f to %.8f", energy0 - model.ref_energy, energy2 - model.ref_energy)
diff --git a/qmp/algorithms/vmc.py b/qmp/algorithms/vmc.py
index 67e84a5d..b68fd98d 100644
--- a/qmp/algorithms/vmc.py
+++ b/qmp/algorithms/vmc.py
@@ -12,6 +12,14 @@
 from ..utility.action_dict import action_dict
 
 
+def _get_devices_from_context(context: RuntimeContext) -> list[torch.device]:
+    """
+    Get the devices list from the runtime context.
+    """
+    assert context.devices is not None
+    return context.devices
+
+
 @dataclasses.dataclass
 class VmcConfig:
     """
@@ -47,6 +55,9 @@ def main(
             runtime_config.optimizer, network.parameters(), checkpoint_data.get("optimizer")
         )
 
+        # Get devices for multi-GPU computation
+        devices = _get_devices_from_context(context)
+
         logging.info(
             "Arguments Summary: Sampling Count: %d, Relative Count: %d, Local Steps: %d, Unique: %s",
             self.sampling_count,
@@ -79,7 +90,7 @@ def main(
             else:
                 configs_src = configs_i
                 configs_dst = torch.cat(
-                    [configs_i, model.find_relative(configs_i, psi_i, self.relative_count - len(configs_i))]
+                    [configs_i, model.find_relative(configs_i, psi_i, self.relative_count - len(configs_i), devices=devices)]
                 )
             logging.info("Relative configurations calculated, count: %d", len(configs_dst))
 
@@ -89,7 +100,7 @@ def closure() -> torch.Tensor:
                 psi_src = network(configs_src)
                 with torch.no_grad():
                     psi_dst = network(configs_dst)
-                    hamiltonian_psi_dst = model.apply_within(configs_dst, psi_dst, configs_src)
+                    hamiltonian_psi_dst = model.apply_within(configs_dst, psi_dst, configs_src, devices)
                 psi_src_reweight = psi_src.conj() * reweight
                 num = psi_src_reweight @ hamiltonian_psi_dst
                 den = psi_src_reweight @ psi_src.detach()
diff --git a/qmp/hamiltonian/hamiltonian.py b/qmp/hamiltonian/hamiltonian.py
index d56b9a3a..102a3f73 100644
--- a/qmp/hamiltonian/hamiltonian.py
+++ b/qmp/hamiltonian/hamiltonian.py
@@ -3,9 +3,28 @@
 """
 
 import os
+import logging
 import platformdirs
 import torch
 import torch.utils.cpp_extension
+from concurrent.futures import ThreadPoolExecutor
+
+
+# Cache for CUDA streams per device
+_stream_cache: dict[str, torch.cuda.Stream] = {}
+
+
+def _get_stream(device: torch.device) -> torch.cuda.Stream | None:
+    """
+    Get or create a CUDA stream for the given device.
+    Returns None for CPU devices.
+    """
+    if device.type != "cuda":
+        return None
+    device_key = str(device)
+    if device_key not in _stream_cache:
+        _stream_cache[device_key] = torch.cuda.Stream(device=device)
+    return _stream_cache[device_key]
 
 
 class Hamiltonian:
@@ -111,6 +130,12 @@ def __init__(
             case _:
                 raise ValueError(f"Unknown kind: {kind}")
 
+        # Multi-device support: cache tensors on each device
+        # These dictionaries store site, kind, coef tensors for each device
+        self._site_dict: dict[str, torch.Tensor] = {}
+        self._kind_dict: dict[str, torch.Tensor] = {}
+        self._coef_dict: dict[str, torch.Tensor] = {}
+
     def _sort_site_kind_coef(self) -> None:
         """
         Reorder the site, kind, and coefficient tensors in descending order of the norm of the coefficients.
@@ -120,9 +145,22 @@ def _sort_site_kind_coef(self) -> None:
         self.kind = self.kind[order]
         self.coef = self.coef[order]
 
+    def _prepare_data_for_device(self, device: torch.device) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
+        """
+        Prepare the site, kind, and coefficient tensors for computation on a specific device.
+        Uses caching to avoid repeated transfers for the same device.
+        """
+        device_key = str(device)
+        if device_key not in self._site_dict:
+            self._site_dict[device_key] = self.site.to(device=device).contiguous()
+            self._kind_dict[device_key] = self.kind.to(device=device).contiguous()
+            self._coef_dict[device_key] = self.coef.to(device=device).contiguous()
+        return self._site_dict[device_key], self._kind_dict[device_key], self._coef_dict[device_key]
+
     def _prepare_data(self, device: torch.device) -> None:
         """
         Prepare the site, kind, and coefficient tensors for computation on the given device.
+        Deprecated: use _prepare_data_for_device instead for multi-device support.
         """
         self.site = self.site.to(device=device).contiguous()
         self.kind = self.kind.to(device=device).contiguous()
@@ -133,6 +171,7 @@ def apply_within(
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_j: torch.Tensor,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
         """
         Applies the Hamiltonian to the given vector.
@@ -145,21 +184,105 @@ def apply_within(
             A complex64 tensor of shape [batch_size_i] representing the input amplitudes on the given configurations.
         configs_j : torch.Tensor
             A uint8 tensor of shape [batch_size_j, n_qubytes] representing the output configurations.
+        devices : list[torch.device] | None
+            A list of devices to use for computation. If None, uses the device of configs_i.
 
         Returns
         -------
         torch.Tensor
             A tensor of shape [batch_size_j] representing the output amplitudes on the given configurations.
         """
-        self._prepare_data(configs_i.device)
-        _apply_within = getattr(
-            self._load_module(configs_i.device.type, configs_i.size(1), self.particle_cut),
-            "apply_within",
-        )
-        psi_j = torch.view_as_complex(
-            _apply_within(configs_i, torch.view_as_real(psi_i), configs_j, self.site, self.kind, self.coef)
-        )
-        return psi_j
+        if devices is None or len(devices) <= 1:
+            # Single device case: use original implementation
+            device = devices[0] if devices else configs_i.device
+            site, kind, coef = self._prepare_data_for_device(device)
+            configs_i_dev = configs_i.to(device=device)
+            psi_i_dev = psi_i.to(device=device)
+            configs_j_dev = configs_j.to(device=device)
+            _apply_within = getattr(
+                self._load_module(device.type, configs_i_dev.size(1), self.particle_cut),
+                "apply_within",
+            )
+            psi_j = torch.view_as_complex(
+                _apply_within(configs_i_dev, torch.view_as_real(psi_i_dev), configs_j_dev, site, kind, coef)
+            )
+            return psi_j
+
+        # Multi-device case: parallel execution using ThreadPoolExecutor + CUDA streams
+        n_devices = len(devices)
+        batch_size_i = configs_i.size(0)
+        chunk_size = batch_size_i // n_devices
+        remainder = batch_size_i % n_devices
+
+        # Pre-prepare Hamiltonian data for all devices
+        for device in devices:
+            self._prepare_data_for_device(device)
+
+        result_device = devices[0]
+        hamiltonian_self = self  # Capture self for use in worker function
+
+        def compute_chunk(idx: int, device: torch.device) -> torch.Tensor | None:
+            """Worker function to compute a chunk on a specific device."""
+            start_idx = idx * chunk_size + min(idx, remainder)
+            end_idx = start_idx + chunk_size + (1 if idx < remainder else 0)
+
+            if start_idx >= end_idx:
+                return None
+
+            site, kind, coef = hamiltonian_self._prepare_data_for_device(device)
+            stream = _get_stream(device)
+
+            if stream is not None:
+                # CUDA device: use thread + stream for true parallel execution
+                with torch.cuda.device(device.index):
+                    with torch.cuda.stream(stream):
+                        configs_i_chunk = configs_i[start_idx:end_idx].to(device=device, non_blocking=True)
+                        psi_i_chunk = psi_i[start_idx:end_idx].to(device=device, non_blocking=True)
+                        configs_j_dev = configs_j.to(device=device, non_blocking=True)
+
+                        _apply_within = getattr(
+                            hamiltonian_self._load_module(device.type, configs_i_chunk.size(1), hamiltonian_self.particle_cut),
+                            "apply_within",
+                        )
+                        psi_j_chunk = torch.view_as_complex(
+                            _apply_within(configs_i_chunk, torch.view_as_real(psi_i_chunk), configs_j_dev, site, kind, coef)
+                        )
+                        # Synchronize within thread to ensure result is complete
+                        stream.synchronize()
+                        return psi_j_chunk.to(device=result_device)
+            else:
+                # CPU device: execute directly
+                configs_i_chunk = configs_i[start_idx:end_idx].to(device=device)
+                psi_i_chunk = psi_i[start_idx:end_idx].to(device=device)
+                configs_j_dev = configs_j.to(device=device)
+
+                _apply_within = getattr(
+                    hamiltonian_self._load_module(device.type, configs_i_chunk.size(1), hamiltonian_self.particle_cut),
+                    "apply_within",
+                )
+                psi_j_chunk = torch.view_as_complex(
+                    _apply_within(configs_i_chunk, torch.view_as_real(psi_i_chunk), configs_j_dev, site, kind, coef)
+                )
+                return psi_j_chunk.to(device=result_device)
+
+        # Parallel execution using ThreadPoolExecutor
+        results: list[torch.Tensor] = []
+        with ThreadPoolExecutor(max_workers=n_devices) as executor:
+            futures = [executor.submit(compute_chunk, idx, device) for idx, device in enumerate(devices)]
+            for future in futures:
+                result = future.result()
+                if result is not None:
+                    results.append(result)
+
+        # Accumulate results
+        if len(results) == 0:
+            return torch.zeros(configs_j.size(0), dtype=torch.complex64, device=result_device)
+
+        final_result = torch.zeros(configs_j.size(0), dtype=torch.complex64, device=result_device)
+        for psi_j_chunk in results:
+            final_result = final_result + psi_j_chunk
+
+        return final_result
 
     def find_relative(
         self,
@@ -167,6 +290,7 @@ def find_relative(
         psi_i: torch.Tensor,
         count_selected: int,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
         """
         Find relative configurations to the given configurations.
@@ -181,6 +305,8 @@ def find_relative(
             The number of selected configurations to be returned.
         configs_exclude : torch.Tensor, optional
             A uint8 tensor of shape [batch_size_exclude, n_qubytes] representing the configurations to be excluded from the result, by default None
+        devices : list[torch.device] | None
+            A list of devices to use for computation. If None, uses the device of configs_i.
 
         Returns
         -------
@@ -190,21 +316,108 @@ def find_relative(
         """
         if configs_exclude is None:
             configs_exclude = configs_i
-        self._prepare_data(configs_i.device)
-        _find_relative = getattr(
-            self._load_module(configs_i.device.type, configs_i.size(1), self.particle_cut),
-            "find_relative",
-        )
-        configs_j = _find_relative(
-            configs_i, torch.view_as_real(psi_i), count_selected, self.site, self.kind, self.coef, configs_exclude
-        )
-        return configs_j
+
+        if devices is None or len(devices) <= 1:
+            device = devices[0] if devices else configs_i.device
+            site, kind, coef = self._prepare_data_for_device(device)
+            configs_i_dev = configs_i.to(device=device)
+            psi_i_dev = psi_i.to(device=device)
+            configs_exclude_dev = configs_exclude.to(device=device)
+            _find_relative = getattr(
+                self._load_module(device.type, configs_i_dev.size(1), self.particle_cut),
+                "find_relative",
+            )
+            configs_j = _find_relative(
+                configs_i_dev, torch.view_as_real(psi_i_dev), count_selected, site, kind, coef, configs_exclude_dev
+            )
+            return configs_j
+
+        # Multi-device case: parallel execution using ThreadPoolExecutor + CUDA streams
+        n_devices = len(devices)
+        batch_size_i = configs_i.size(0)
+        chunk_size = batch_size_i // n_devices
+        remainder = batch_size_i % n_devices
+
+        # Pre-prepare Hamiltonian data for all devices
+        for device in devices:
+            self._prepare_data_for_device(device)
+
+        result_device = devices[0]
+        hamiltonian_self = self
+
+        def compute_chunk(idx: int, device: torch.device) -> torch.Tensor | None:
+            """Worker function to compute find_relative on a specific device."""
+            start_idx = idx * chunk_size + min(idx, remainder)
+            end_idx = start_idx + chunk_size + (1 if idx < remainder else 0)
+
+            if start_idx >= end_idx:
+                return None
+
+            site, kind, coef = hamiltonian_self._prepare_data_for_device(device)
+            stream = _get_stream(device)
+
+            if stream is not None:
+                with torch.cuda.device(device.index):
+                    with torch.cuda.stream(stream):
+                        configs_i_chunk = configs_i[start_idx:end_idx].to(device=device, non_blocking=True)
+                        psi_i_chunk = psi_i[start_idx:end_idx].to(device=device, non_blocking=True)
+                        configs_exclude_dev = configs_exclude.to(device=device, non_blocking=True)
+
+                        _find_relative = getattr(
+                            hamiltonian_self._load_module(device.type, configs_i_chunk.size(1), hamiltonian_self.particle_cut),
+                            "find_relative",
+                        )
+                        configs_j_chunk = _find_relative(
+                            configs_i_chunk, torch.view_as_real(psi_i_chunk), count_selected, site, kind, coef, configs_exclude_dev
+                        )
+                        stream.synchronize()
+                        return configs_j_chunk.to(device=result_device)
+            else:
+                configs_i_chunk = configs_i[start_idx:end_idx].to(device=device)
+                psi_i_chunk = psi_i[start_idx:end_idx].to(device=device)
+                configs_exclude_dev = configs_exclude.to(device=device)
+
+                _find_relative = getattr(
+                    hamiltonian_self._load_module(device.type, configs_i_chunk.size(1), hamiltonian_self.particle_cut),
+                    "find_relative",
+                )
+                configs_j_chunk = _find_relative(
+                    configs_i_chunk, torch.view_as_real(psi_i_chunk), count_selected, site, kind, coef, configs_exclude_dev
+                )
+                return configs_j_chunk.to(device=result_device)
+
+        # Parallel execution using ThreadPoolExecutor
+        results: list[torch.Tensor] = []
+        with ThreadPoolExecutor(max_workers=n_devices) as executor:
+            futures = [executor.submit(compute_chunk, idx, device) for idx, device in enumerate(devices)]
+            for future in futures:
+                result = future.result()
+                if result is not None:
+                    results.append(result)
+
+        # Handle empty results
+        if len(results) == 0:
+            return torch.empty(0, configs_i.size(1), dtype=configs_i.dtype, device=result_device)
+
+        # Merge with sorted=True to preserve importance ordering
+        all_configs = torch.cat(results, dim=0)
+        unique_configs = torch.unique(all_configs, sorted=True, dim=0)
+
+        # Exclude configs that appear in configs_exclude
+        configs_exclude_dev = configs_exclude.to(device=result_device)
+        exclude_mask = (unique_configs.unsqueeze(1) == configs_exclude_dev.unsqueeze(0)).all(dim=-1).any(dim=-1)
+        filtered_configs = unique_configs[~exclude_mask]
+
+        if filtered_configs.size(0) > count_selected:
+            filtered_configs = filtered_configs[:count_selected]
+        return filtered_configs
 
     def list_relative(
         self,
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> tuple[torch.Tensor, torch.Tensor]:
         """
         List all unique relative configurations and their accumulated amplitudes.
@@ -217,6 +430,8 @@ def list_relative(
             Input amplitudes (complex64).
         configs_exclude : torch.Tensor, optional
             Configurations to exclude from the result. Defaults to configs_i.
+        devices : list[torch.device] | None
+            A list of devices to use for computation. If None, uses the device of configs_i.
 
         Returns
         -------
@@ -226,17 +441,111 @@ def list_relative(
         """
         if configs_exclude is None:
             configs_exclude = configs_i
-        self._prepare_data(configs_i.device)
-        _list_relative = getattr(
-            self._load_module(configs_i.device.type, configs_i.size(1), self.particle_cut),
-            "list_relative",
-        )
-        configs_j, psi_j_real = _list_relative(
-            configs_i, torch.view_as_real(psi_i), self.site, self.kind, self.coef, configs_exclude
-        )
-        return configs_j, torch.view_as_complex(psi_j_real)
-
-    def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
+
+        if devices is None or len(devices) <= 1:
+            device = devices[0] if devices else configs_i.device
+            site, kind, coef = self._prepare_data_for_device(device)
+            configs_i_dev = configs_i.to(device=device)
+            psi_i_dev = psi_i.to(device=device)
+            configs_exclude_dev = configs_exclude.to(device=device)
+            _list_relative = getattr(
+                self._load_module(device.type, configs_i_dev.size(1), self.particle_cut),
+                "list_relative",
+            )
+            configs_j, psi_j_real = _list_relative(
+                configs_i_dev, torch.view_as_real(psi_i_dev), site, kind, coef, configs_exclude_dev
+            )
+            return configs_j, torch.view_as_complex(psi_j_real)
+
+        # Multi-device case: parallel execution using ThreadPoolExecutor + CUDA streams
+        n_devices = len(devices)
+        batch_size_i = configs_i.size(0)
+        chunk_size = batch_size_i // n_devices
+        remainder = batch_size_i % n_devices
+
+        # Pre-prepare Hamiltonian data for all devices
+        for device in devices:
+            self._prepare_data_for_device(device)
+
+        result_device = devices[0]
+        hamiltonian_self = self
+
+        def compute_chunk(idx: int, device: torch.device) -> tuple[torch.Tensor, torch.Tensor] | None:
+            """Worker function to compute list_relative on a specific device."""
+            start_idx = idx * chunk_size + min(idx, remainder)
+            end_idx = start_idx + chunk_size + (1 if idx < remainder else 0)
+
+            if start_idx >= end_idx:
+                return None
+
+            site, kind, coef = hamiltonian_self._prepare_data_for_device(device)
+            stream = _get_stream(device)
+
+            if stream is not None:
+                with torch.cuda.device(device.index):
+                    with torch.cuda.stream(stream):
+                        configs_i_chunk = configs_i[start_idx:end_idx].to(device=device, non_blocking=True)
+                        psi_i_chunk = psi_i[start_idx:end_idx].to(device=device, non_blocking=True)
+                        configs_exclude_dev = configs_exclude.to(device=device, non_blocking=True)
+
+                        _list_relative = getattr(
+                            hamiltonian_self._load_module(device.type, configs_i_chunk.size(1), hamiltonian_self.particle_cut),
+                            "list_relative",
+                        )
+                        configs_j_chunk, psi_j_chunk_real = _list_relative(
+                            configs_i_chunk, torch.view_as_real(psi_i_chunk), site, kind, coef, configs_exclude_dev
+                        )
+                        psi_j_chunk = torch.view_as_complex(psi_j_chunk_real)
+                        stream.synchronize()
+                        return (configs_j_chunk.to(device=result_device), psi_j_chunk.to(device=result_device))
+            else:
+                configs_i_chunk = configs_i[start_idx:end_idx].to(device=device)
+                psi_i_chunk = psi_i[start_idx:end_idx].to(device=device)
+                configs_exclude_dev = configs_exclude.to(device=device)
+
+                _list_relative = getattr(
+                    hamiltonian_self._load_module(device.type, configs_i_chunk.size(1), hamiltonian_self.particle_cut),
+                    "list_relative",
+                )
+                configs_j_chunk, psi_j_chunk_real = _list_relative(
+                    configs_i_chunk, torch.view_as_real(psi_i_chunk), site, kind, coef, configs_exclude_dev
+                )
+                psi_j_chunk = torch.view_as_complex(psi_j_chunk_real)
+                return (configs_j_chunk.to(device=result_device), psi_j_chunk.to(device=result_device))
+
+        # Parallel execution using ThreadPoolExecutor
+        results_configs: list[torch.Tensor] = []
+        results_psi: list[torch.Tensor] = []
+        with ThreadPoolExecutor(max_workers=n_devices) as executor:
+            futures = [executor.submit(compute_chunk, idx, device) for idx, device in enumerate(devices)]
+            for future in futures:
+                result = future.result()
+                if result is not None:
+                    results_configs.append(result[0])
+                    results_psi.append(result[1])
+
+        # Handle empty results
+        if len(results_configs) == 0:
+            return torch.empty(0, configs_i.size(1), dtype=configs_i.dtype, device=result_device), \
+                   torch.empty(0, dtype=psi_i.dtype, device=result_device)
+
+        # Merge and deduplicate
+        all_configs = torch.cat(results_configs, dim=0)
+        all_psi = torch.cat(results_psi, dim=0)
+
+        unique_configs, inverse_indices = torch.unique(all_configs, return_inverse=True, dim=0)
+        unique_psi = torch.zeros(unique_configs.size(0), dtype=all_psi.dtype, device=result_device)
+        unique_psi.scatter_add_(0, inverse_indices, all_psi)
+
+        # Exclude configs that appear in configs_exclude
+        configs_exclude_dev = configs_exclude.to(device=result_device)
+        exclude_mask = (unique_configs.unsqueeze(1) == configs_exclude_dev.unsqueeze(0)).all(dim=-1).any(dim=-1)
+        filtered_configs = unique_configs[~exclude_mask]
+        filtered_psi = unique_psi[~exclude_mask]
+
+        return filtered_configs, filtered_psi
+
+    def diagonal_term(self, configs: torch.Tensor, devices: list[torch.device] | None = None) -> torch.Tensor:
         """
         Get the diagonal term of the Hamiltonian for the given configurations.
 
@@ -244,16 +553,88 @@ def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
         ----------
         configs : torch.Tensor
             A uint8 tensor of shape [batch_size, n_qubytes] representing the input configurations.
+        devices : list[torch.device] | None
+            A list of devices to use for computation. If None, uses the device of configs.
 
         Returns
         -------
         torch.Tensor
             A complex64 tensor of shape [batch_size] representing the diagonal term of the Hamiltonian for the given configurations.
         """
-        self._prepare_data(configs.device)
-        _diagonal_term = getattr(
-            self._load_module(configs.device.type, configs.size(1), self.particle_cut),
-            "diagonal_term",
-        )
-        psi_result = torch.view_as_complex(_diagonal_term(configs, self.site, self.kind, self.coef))
-        return psi_result
+        if devices is None or len(devices) <= 1:
+            device = devices[0] if devices else configs.device
+            site, kind, coef = self._prepare_data_for_device(device)
+            configs_dev = configs.to(device=device)
+            _diagonal_term = getattr(
+                self._load_module(device.type, configs_dev.size(1), self.particle_cut),
+                "diagonal_term",
+            )
+            psi_result = torch.view_as_complex(_diagonal_term(configs_dev, site, kind, coef))
+            return psi_result
+
+        # Multi-device case: parallel execution using ThreadPoolExecutor + CUDA streams
+        n_devices = len(devices)
+        batch_size = configs.size(0)
+        chunk_size = batch_size // n_devices
+        remainder = batch_size % n_devices
+
+        # Pre-prepare Hamiltonian data for all devices
+        for device in devices:
+            self._prepare_data_for_device(device)
+
+        result_device = devices[0]
+        hamiltonian_self = self
+
+        def compute_chunk(idx: int, device: torch.device) -> tuple[torch.Tensor, int] | None:
+            """Worker function to compute diagonal_term on a specific device."""
+            start_idx = idx * chunk_size + min(idx, remainder)
+            end_idx = start_idx + chunk_size + (1 if idx < remainder else 0)
+
+            if start_idx >= end_idx:
+                return None
+
+            site, kind, coef = hamiltonian_self._prepare_data_for_device(device)
+            stream = _get_stream(device)
+
+            if stream is not None:
+                with torch.cuda.device(device.index):
+                    with torch.cuda.stream(stream):
+                        configs_chunk = configs[start_idx:end_idx].to(device=device, non_blocking=True)
+
+                        _diagonal_term = getattr(
+                            hamiltonian_self._load_module(device.type, configs_chunk.size(1), hamiltonian_self.particle_cut),
+                            "diagonal_term",
+                        )
+                        psi_chunk = torch.view_as_complex(_diagonal_term(configs_chunk, site, kind, coef))
+                        stream.synchronize()
+                        return (psi_chunk.to(device=result_device), start_idx)
+            else:
+                configs_chunk = configs[start_idx:end_idx].to(device=device)
+
+                _diagonal_term = getattr(
+                    hamiltonian_self._load_module(device.type, configs_chunk.size(1), hamiltonian_self.particle_cut),
+                    "diagonal_term",
+                )
+                psi_chunk = torch.view_as_complex(_diagonal_term(configs_chunk, site, kind, coef))
+                return (psi_chunk.to(device=result_device), start_idx)
+
+        # Parallel execution using ThreadPoolExecutor
+        pending_results: list[tuple[torch.Tensor, int]] = []
+        with ThreadPoolExecutor(max_workers=n_devices) as executor:
+            futures = [executor.submit(compute_chunk, idx, device) for idx, device in enumerate(devices)]
+            for future in futures:
+                result = future.result()
+                if result is not None:
+                    pending_results.append(result)
+
+        # Sort by start_idx to maintain order
+        pending_results.sort(key=lambda x: x[1])
+
+        results: list[torch.Tensor] = []
+        for psi_chunk, start_idx in pending_results:
+            results.append(psi_chunk)
+
+        if len(results) == 0:
+            return torch.empty(0, dtype=torch.complex64, device=result_device)
+
+        return torch.cat(results, dim=0)
\ No newline at end of file
diff --git a/qmp/models/fcidump.py b/qmp/models/fcidump.py
index a5cbd853..cfc4e519 100644
--- a/qmp/models/fcidump.py
+++ b/qmp/models/fcidump.py
@@ -174,8 +174,14 @@ def __init__(self, args: ModelConfig) -> None:
                 self.ref_energy = 0
         logging.info("Reference energy for model is %.10f", self.ref_energy)
 
-    def apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j)
+    def apply_within(
+        self,
+        configs_i: torch.Tensor,
+        psi_i: torch.Tensor,
+        configs_j: torch.Tensor,
+        devices: list[torch.device] | None = None,
+    ) -> torch.Tensor:
+        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j, devices)
 
     def find_relative(
         self,
@@ -183,19 +189,21 @@ def find_relative(
         psi_i: torch.Tensor,
         count_selected: int,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
-        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude)
+        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude, devices)
 
     def list_relative(
         self,
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> tuple[torch.Tensor, torch.Tensor]:
-        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude)
+        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude, devices)
 
-    def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.diagonal_term(configs)
+    def diagonal_term(self, configs: torch.Tensor, devices: list[torch.device] | None = None) -> torch.Tensor:
+        return self.hamiltonian.diagonal_term(configs, devices)
 
     def show_config(self, config: torch.Tensor) -> str:
         string = "".join(f"{i:08b}"[::-1] for i in config.cpu().numpy())
diff --git a/qmp/models/hubbard.py b/qmp/models/hubbard.py
index c89828cc..28f8a711 100644
--- a/qmp/models/hubbard.py
+++ b/qmp/models/hubbard.py
@@ -118,8 +118,14 @@ def __init__(self, args: ModelConfig):
         self.hamiltonian: Hamiltonian = Hamiltonian(hamiltonian_dict, kind="fermi")
         logging.info("Internal Hamiltonian representation successfully created.")
 
-    def apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j)
+    def apply_within(
+        self,
+        configs_i: torch.Tensor,
+        psi_i: torch.Tensor,
+        configs_j: torch.Tensor,
+        devices: list[torch.device] | None = None,
+    ) -> torch.Tensor:
+        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j, devices)
 
     def find_relative(
         self,
@@ -127,19 +133,21 @@ def find_relative(
         psi_i: torch.Tensor,
         count_selected: int,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
-        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude)
+        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude, devices)
 
     def list_relative(
         self,
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> tuple[torch.Tensor, torch.Tensor]:
-        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude)
+        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude, devices)
 
-    def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.diagonal_term(configs)
+    def diagonal_term(self, configs: torch.Tensor, devices: list[torch.device] | None = None) -> torch.Tensor:
+        return self.hamiltonian.diagonal_term(configs, devices)
 
     def show_config(self, config: torch.Tensor) -> str:
         string = "".join(f"{i:08b}"[::-1] for i in config.cpu().numpy())
diff --git a/qmp/models/ising.py b/qmp/models/ising.py
index be2f6367..0dd0d99a 100644
--- a/qmp/models/ising.py
+++ b/qmp/models/ising.py
@@ -184,8 +184,14 @@ def __init__(self, args: ModelConfig) -> None:
         self.hamiltonian: Hamiltonian = Hamiltonian(hamiltonian_dict, kind="bose2")
         logging.info("Internal Hamiltonian representation successfully created.")
 
-    def apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j)
+    def apply_within(
+        self,
+        configs_i: torch.Tensor,
+        psi_i: torch.Tensor,
+        configs_j: torch.Tensor,
+        devices: list[torch.device] | None = None,
+    ) -> torch.Tensor:
+        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j, devices)
 
     def find_relative(
         self,
@@ -193,19 +199,21 @@ def find_relative(
         psi_i: torch.Tensor,
         count_selected: int,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
-        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude)
+        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude, devices)
 
     def list_relative(
         self,
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> tuple[torch.Tensor, torch.Tensor]:
-        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude)
+        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude, devices)
 
-    def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.diagonal_term(configs)
+    def diagonal_term(self, configs: torch.Tensor, devices: list[torch.device] | None = None) -> torch.Tensor:
+        return self.hamiltonian.diagonal_term(configs, devices)
 
     def show_config(self, config: torch.Tensor) -> str:
         string = "".join(f"{i:08b}"[::-1] for i in config.cpu().numpy())
diff --git a/qmp/models/openfermion.py b/qmp/models/openfermion.py
index 270e9ee0..b2c4f277 100644
--- a/qmp/models/openfermion.py
+++ b/qmp/models/openfermion.py
@@ -58,8 +58,14 @@ def __init__(self, args: ModelConfig) -> None:
         )
         logging.info("Internal Hamiltonian representation successfully created.")
 
-    def apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j)
+    def apply_within(
+        self,
+        configs_i: torch.Tensor,
+        psi_i: torch.Tensor,
+        configs_j: torch.Tensor,
+        devices: list[torch.device] | None = None,
+    ) -> torch.Tensor:
+        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j, devices)
 
     def find_relative(
         self,
@@ -67,19 +73,21 @@ def find_relative(
         psi_i: torch.Tensor,
         count_selected: int,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
-        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude)
+        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude, devices)
 
     def list_relative(
         self,
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> tuple[torch.Tensor, torch.Tensor]:
-        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude)
+        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude, devices)
 
-    def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.diagonal_term(configs)
+    def diagonal_term(self, configs: torch.Tensor, devices: list[torch.device] | None = None) -> torch.Tensor:
+        return self.hamiltonian.diagonal_term(configs, devices)
 
     def show_config(self, config: torch.Tensor) -> str:
         string = "".join(f"{i:08b}"[::-1] for i in config.cpu().numpy())
diff --git a/qmp/models/optimized_basis.py b/qmp/models/optimized_basis.py
index dbc62493..e7bf60c9 100644
--- a/qmp/models/optimized_basis.py
+++ b/qmp/models/optimized_basis.py
@@ -65,8 +65,14 @@ def __init__(self, args: ModelConfig) -> None:
         if self.orbit_unitary is not None:
             logging.info("Orbit unitary matrix loaded (shape: %s)", self.orbit_unitary.shape)
 
-    def apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j)
+    def apply_within(
+        self,
+        configs_i: torch.Tensor,
+        psi_i: torch.Tensor,
+        configs_j: torch.Tensor,
+        devices: list[torch.device] | None = None,
+    ) -> torch.Tensor:
+        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j, devices)
 
     def find_relative(
         self,
@@ -74,19 +80,21 @@ def find_relative(
         psi_i: torch.Tensor,
         count_selected: int,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
-        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude)
+        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude, devices)
 
     def list_relative(
         self,
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> tuple[torch.Tensor, torch.Tensor]:
-        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude)
+        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude, devices)
 
-    def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.diagonal_term(configs)
+    def diagonal_term(self, configs: torch.Tensor, devices: list[torch.device] | None = None) -> torch.Tensor:
+        return self.hamiltonian.diagonal_term(configs, devices)
 
     def show_config(self, config: torch.Tensor) -> str:
         string = "".join(f"{i:08b}"[::-1] for i in config.cpu().numpy())
diff --git a/qmp/models/pyscf.py b/qmp/models/pyscf.py
index 7fe76d82..8516e66c 100644
--- a/qmp/models/pyscf.py
+++ b/qmp/models/pyscf.py
@@ -165,8 +165,14 @@ def __init__(self, args: ModelConfig) -> None:
             self.ref_energy,
         )
 
-    def apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j)
+    def apply_within(
+        self,
+        configs_i: torch.Tensor,
+        psi_i: torch.Tensor,
+        configs_j: torch.Tensor,
+        devices: list[torch.device] | None = None,
+    ) -> torch.Tensor:
+        return self.hamiltonian.apply_within(configs_i, psi_i, configs_j, devices)
 
     def find_relative(
         self,
@@ -174,19 +180,21 @@ def find_relative(
         psi_i: torch.Tensor,
         count_selected: int,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
-        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude)
+        return self.hamiltonian.find_relative(configs_i, psi_i, count_selected, configs_exclude, devices)
 
     def list_relative(
         self,
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> tuple[torch.Tensor, torch.Tensor]:
-        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude)
+        return self.hamiltonian.list_relative(configs_i, psi_i, configs_exclude, devices)
 
-    def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
-        return self.hamiltonian.diagonal_term(configs)
+    def diagonal_term(self, configs: torch.Tensor, devices: list[torch.device] | None = None) -> torch.Tensor:
+        return self.hamiltonian.diagonal_term(configs, devices)
 
     def show_config(self, config: torch.Tensor) -> str:
         string = "".join(f"{i:08b}"[::-1] for i in config.cpu().numpy())
diff --git a/qmp/utility/context.py b/qmp/utility/context.py
index 44049a66..87e96235 100644
--- a/qmp/utility/context.py
+++ b/qmp/utility/context.py
@@ -16,7 +16,7 @@
 from .optimizer import migrate_optimizer
 
 
-DACITE_CAST = [pathlib.Path, torch.device, tuple]
+DACITE_CAST = [pathlib.Path, torch.device, tuple, list]
 
 
 @dataclasses.dataclass
@@ -31,8 +31,11 @@ class RuntimeContext:
     random_seed: int | None = None
     # The interval to save the checkpoint
     checkpoint_interval: int = 5
-    # The device to run on
+    # The device to run on (for network and checkpoint storage)
     device: torch.device = torch.device(type="cuda", index=0)
+    # The devices to use for Hamiltonian computation (multi-GPU support)
+    # If None, defaults to [device]
+    devices: list[torch.device] | None = None
     # The dtype of the network, leave empty to skip modifying the dtype
     dtype: str | torch.dtype | None = None
     # The maximum absolute step for the process, leave empty to loop forever
@@ -55,6 +58,9 @@ def __post_init__(self) -> None:
                     raise ValueError(f"Unsupported dtype: {self.dtype}")
         if self.max_absolute_step is not None and self.max_relative_step is not None:
             raise ValueError("Both max_absolute_step and max_relative_step are set, please set only one of them.")
+        # Initialize devices to default if not provided
+        if self.devices is None:
+            self.devices = [self.device]
 
     def folder(self) -> pathlib.Path:
         """
diff --git a/qmp/utility/model_dict.py b/qmp/utility/model_dict.py
index 2df8bb73..2bb12896 100644
--- a/qmp/utility/model_dict.py
+++ b/qmp/utility/model_dict.py
@@ -125,7 +125,13 @@ def __init__(self, config: ModelConfig) -> None:
             The config of model.
         """
 
-    def apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j: torch.Tensor) -> torch.Tensor:
+    def apply_within(
+        self,
+        configs_i: torch.Tensor,
+        psi_i: torch.Tensor,
+        configs_j: torch.Tensor,
+        devices: list[torch.device] | None = None,
+    ) -> torch.Tensor:
         """
         Applies the Hamiltonian to the given vector.
 
@@ -137,6 +143,8 @@ def apply_within(self, configs_i: torch.Tensor, psi_i: torch.Tensor, configs_j:
             The amplitudes of the configurations.
         configs_j : torch.Tensor
             The configurations subspace for the result of the Hamiltonian application.
+        devices : list[torch.device] | None
+            A list of devices to use for computation. If None, uses the default device.
 
         Returns
         -------
@@ -150,6 +158,7 @@ def find_relative(
         psi_i: torch.Tensor,
         count_selected: int,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> torch.Tensor:
         """
         Find relative configurations to the given configurations.
@@ -164,6 +173,8 @@ def find_relative(
             The number of relative configurations to find.
         configs_exclude : torch.Tensor, optional
             The configurations to exclude from the result.
+        devices : list[torch.device] | None
+            A list of devices to use for computation. If None, uses the default device.
 
         Returns
         -------
@@ -176,6 +187,7 @@ def list_relative(
         configs_i: torch.Tensor,
         psi_i: torch.Tensor,
         configs_exclude: torch.Tensor | None = None,
+        devices: list[torch.device] | None = None,
     ) -> tuple[torch.Tensor, torch.Tensor]:
         """
         List all unique relative configurations and their accumulated amplitudes.
@@ -188,6 +200,8 @@ def list_relative(
             Input amplitudes (complex64).
         configs_exclude : torch.Tensor, optional
             Configurations to exclude from the result. Defaults to configs_i.
+        devices : list[torch.device] | None
+            A list of devices to use for computation. If None, uses the default device.
 
         Returns
         -------
@@ -196,7 +210,7 @@ def list_relative(
             and psi_j are their summed amplitudes from all connected paths.
         """
 
-    def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
+    def diagonal_term(self, configs: torch.Tensor, devices: list[torch.device] | None = None) -> torch.Tensor:
         """
         Calculate the diagonal term for the given configurations.
 
@@ -204,6 +218,8 @@ def diagonal_term(self, configs: torch.Tensor) -> torch.Tensor:
         ----------
         configs : torch.Tensor
             The configurations to calculate the diagonal term for.
+        devices : list[torch.device] | None
+            A list of devices to use for computation. If None, uses the default device.
 
         Returns
         -------