Feat/modular julia pardiso (#310)

* feat(postprocess): Add functionality to export system data for Pardiso/Julia band structure calculations and include related tests.

* feat: simplify band calculation by using solve_eigen_at_k directly

- Remove redundant solve_eigen_dense function and fixed band window logic
- Replace custom eigenvalue solving with existing solve_eigen_at_k function
- Update example notebook with new imports and execution outputs

* docs: improve dptb to Pardiso example notebook

- Fix root_dir path from Pardiso_teach to To_pardiso
- Add sys.path manipulation for proper module imports
- Expand Julia installation instructions with detailed steps
- Add MKL and Pardiso environment variable configuration
- Include troubleshooting section for common issues
- Clear unnecessary cell outputs and update execution counts

* Refactor Julia backend: Modular architecture with Dense/Pardiso solver support

* feat(pardiso): enforce MKL Pardiso and add Apple Silicon compatibility check

Replace generic PardisoSolver with MKLPardisoSolver for consistent performance and add explicit architecture check for Apple Silicon systems to prevent runtime failures with helpful guidance to use numpy solver instead.

* feat(postprocess): implement modular Julia backend with HDF5 output

- Implement `dos_calculation.jl` for Density of States calculations
- Update `main.jl` to switch between `band` and `dos` tasks
- Refactor `band_calculation.jl` to export `bandstructure.h5` using native HDF5
- Add incremental `bands.dat` text output for real-time tracking
- Add verification notebook `examples/To_pardiso/dptb_to_Pardiso_new.ipynb`

* Refactor: Modular Julia backend and Optimized JSON Schema

- Refactored  to use direct ASE integration in .
- Implemented modular Julia backend structure in .
- Optimized : compressed species, removed redundant orbital arrays, aligned formatting.
- Updated example notebook  with new API usage.
- Improved logging: Output logs to both console and .
- Renamed output directory to .

* feat: Add `pdso` command-line entry point for running the Julia/Pardiso solver backend.

* Refactor: Rename Julia backend to 'pardiso', add legacy support, and fix spin handling

- Renamed `dptb/postprocess/julia` to `dptb/postprocess/pardiso`.
- Updated `io.jl` to support legacy `.dat` files via `load_structure_dat`.
- Fixed `load_structure_json` to correctly account for spin degeneracy.
- Updated `pdso.py` and tests to reflect directory rename.

* refactor: Make JSON-based Pardiso export the default `to_pardiso` method and rename the legacy text export to `to_pardiso_debug`.

* update unit test

* feat: Add a comprehensive Pardiso tutorial notebook and enhance the `pdso` entrypoint with ill-conditioned state projection parameters.

* feat: Add Julia installation scripts and documentation for the Pardiso backend, including platform support details.

* refactor: Rename solver functions for clarity, encapsulate Pardiso solver in a module, and add ill-conditioned projection and improved eigenvector handling to the dense solver. (#2)

* fix: address pardiso review comments after merge

* fix: address follow-up pardiso review comments

* fix: stabilize pardiso solver cleanup scope

* chore: remove CLAUDE.md from tracking and add to .gitignore

* docs: extend pardiso examples

* docs: align pardiso documentation

---------

Co-authored-by: YJQ-7 <yjq7472@163.com>

Author: QG-phy

.gitignore

@@ -186,3 +186,4 @@ poetry.toml
 
 # JetBrains
 .idea/
+CLAUDE.md

README.md

@@ -180,6 +180,45 @@ Installing **DeePTB** is straightforward with UV, a fast Python package manager.
   > [!TIP]
   > For easier installation with automatic GPU/CPU detection, use the **From Source** method above instead.
 
+- **Julia Backend** (Optional - for High-Performance Pardiso Solver)
+  
+  > [!NOTE]
+  > **Platform Support**: Pardiso backend currently supports **Linux only**.
+  > - **macOS**: Not supported (Intel MKL limitations)
+  > - **Windows**: Use WSL2 (Windows Subsystem for Linux)
+  
+  If you want to use the Pardiso backend for accelerated band structure calculations:
+  
+  **Automated Installation** (Recommended):
+  ```bash
+  ./install_julia.sh
+  ```
+  
+  **Manual Installation**:
+  1. Install Julia:
+     ```bash
+     # Linux (macOS can install Julia, but the Pardiso backend is not supported)
+     curl -fsSL https://install.julialang.org | sh
+     ```
+  2. Install required packages:
+     ```bash
+     julia install_julia_packages.jl
+     ```
+  
+  **Verify Installation**:
+  ```bash
+  julia -e 'using Pardiso; println("Pardiso available: ", Pardiso.MKL_PARDISO_LOADED[])'
+  ```
+  
+  **Usage**:
+  ```bash
+  dptb pdso band.json -i model.pth -stu structure.vasp -o ./output
+  ```
+  
+  For more details, see:
+  - [Pardiso Backend README](dptb/postprocess/pardiso/README.md)
+  - [Example Tutorial](examples/To_pardiso/README.md)
+
 ## Test code 
 
 To ensure the code is correctly installed, please run the unit tests first:

docs/pardiso_architecture.md

@@ -0,0 +1,212 @@
+"""
+Proposed architecture for Pardiso integration.
+
+This document outlines a better design for integrating high-performance
+eigensolvers (like Pardiso) into DeePTB's postprocessing workflow.
+
+## Current Issues
+
+1. **Tight coupling**: Julia script is standalone, requires manual invocation
+2. **Data redundancy**: Full H(R), S(R) exported to files
+3. **Fragmented workflow**: Python → files → Julia → files → Python
+4. **Limited reusability**: Only works for band calculations
+
+## Proposed Architecture
+
+### Phase 1: Modular Julia Backend (Current PR)
+
+```
+dptb/postprocess/pardiso/
+├── io/
+│   └── io.jl                # Read structure.json/legacy .dat files and H5 matrices
+├── solvers/
+│   ├── dense_solver.jl      # Dense LAPACK solver
+│   └── pardiso_solver.jl    # Pardiso eigenvalue solver
+├── tasks/
+│   ├── band_calculation.jl  # Band structure task
+│   └── dos_calculation.jl   # DOS task
+├── utils/
+│   ├── hamiltonian.jl       # H(R) -> H(k)
+│   └── kpoints.jl           # K-point generation
+└── main.jl                  # Entry point (calls tasks)
+```
+
+**Benefits**:
+- Modular: Easy to add new tasks (DOS, optical, etc.)
+- Testable: Each module can be unit tested
+- Maintainable: Clear separation of concerns
+
+### Phase 2: Python Solver Abstraction
+
+```python
+# dptb/postprocess/unified/solvers/base.py
+from abc import ABC, abstractmethod
+
+class EigenSolver(ABC):
+    \"\"\"Abstract base class for eigenvalue solvers.\"\"\"
+
+    @abstractmethod
+    def solve_at_k(self, H_k: np.ndarray, S_k: np.ndarray,
+                   num_bands: int, **kwargs) -> Tuple[np.ndarray, np.ndarray]:
+        \"\"\"
+        Solve generalized eigenvalue problem H|ψ⟩ = E S|ψ⟩ at single k-point.
+
+        Returns
+        -------
+        eigenvalues : np.ndarray (num_bands,)
+        eigenvectors : np.ndarray (norb, num_bands)
+        \"\"\"
+        pass
+
+    @abstractmethod
+    def solve_batch(self, kpoints: np.ndarray, H_R: dict, S_R: dict,
+                    num_bands: int, **kwargs) -> np.ndarray:
+        \"\"\"
+        Solve for multiple k-points efficiently.
+
+        Returns
+        -------
+        eigenvalues : np.ndarray (nk, num_bands)
+        \"\"\"
+        pass
+
+
+# dptb/postprocess/unified/solvers/numpy_solver.py
+class NumpySolver(EigenSolver):
+    \"\"\"Default numpy/scipy solver (current implementation).\"\"\"
+
+    def solve_at_k(self, H_k, S_k, num_bands, **kwargs):
+        if S_k is None:
+            evals, evecs = np.linalg.eigh(H_k)
+        else:
+            evals, evecs = scipy.linalg.eigh(H_k, S_k)
+        return evals[:num_bands], evecs[:, :num_bands]
+
+    def solve_batch(self, kpoints, H_R, S_R, num_bands, **kwargs):
+        eigenvalues = []
+        for kpt in kpoints:
+            H_k = self._construct_hk(kpt, H_R)
+            S_k = self._construct_hk(kpt, S_R) if S_R else None
+            evals, _ = self.solve_at_k(H_k, S_k, num_bands)
+            eigenvalues.append(evals)
+        return np.array(eigenvalues)
+
+
+# dptb/postprocess/unified/solvers/pardiso_solver.py
+class PardisoSolver(EigenSolver):
+    \"\"\"Julia/Pardiso backend solver for large systems.\"\"\"
+
+    def __init__(self, julia_script_path=None, use_pyjulia=False):
+        self.use_pyjulia = use_pyjulia
+        if use_pyjulia:
+            self._init_pyjulia()
+        else:
+            self.julia_script = julia_script_path or self._default_script_path()
+
+    def solve_batch(self, kpoints, H_R, S_R, num_bands, **kwargs):
+        if self.use_pyjulia:
+            # Direct Julia call (no file I/O)
+            return self._solve_via_pyjulia(kpoints, H_R, S_R, num_bands)
+        else:
+            # File-based workflow (current approach)
+            return self._solve_via_files(kpoints, H_R, S_R, num_bands)
+
+    def _solve_via_files(self, kpoints, H_R, S_R, num_bands):
+        # Export data
+        temp_dir = self._export_data(H_R, S_R, kpoints)
+        # Call Julia script
+        result = subprocess.run(['julia', self.julia_script, temp_dir])
+        # Load results
+        eigenvalues = np.load(os.path.join(temp_dir, 'eigenvalues.npy'))
+        return eigenvalues
+
+
+# Usage in BandAccessor
+class BandAccessor:
+    def __init__(self, tbsystem, solver='numpy'):
+        self.tbsystem = tbsystem
+        self.solver = get_solver(solver)  # Factory function
+
+    def compute(self, **kwargs):
+        hr, sr = self.tbsystem.calculator.get_hr(self.tbsystem.data)
+        eigenvalues = self.solver.solve_batch(
+            self.kpoints, hr, sr, self.num_bands
+        )
+        self.eigenvalues = eigenvalues
+        return self
+```
+
+**Benefits**:
+- Unified interface: All properties (band, DOS, optical) use same solver
+- Pluggable: Easy to add new solvers (cuSOLVER, ELPA, etc.)
+- Backward compatible: Default numpy solver maintains current behavior
+
+### Phase 3: Backend System (Future)
+
+```python
+# User-facing API
+import dptb
+dptb.set_backend('pardiso')  # Global setting
+
+tbsys = TBSystem(data, calculator)
+bands = tbsys.band.compute(kpath)  # Automatically uses Pardiso
+
+# Or per-calculation
+bands = tbsys.band.compute(kpath, solver='pardiso')
+```
+
+## Implementation Roadmap
+
+### Immediate (Current PR)
+1. ✅ Optimize data format (JSON instead of text files)
+2. ✅ Use OrbitalMapper.atom_norb (avoid recomputation)
+3. ✅ Modularize Julia script
+4. ✅ Add Python export tests and lightweight Julia-load coverage
+
+### Short-term (Next PR)
+1. Create `dptb/postprocess/unified/solvers/` module
+2. Implement `EigenSolver` ABC
+3. Refactor current code into `NumpySolver`
+4. Implement `PardisoSolver` (file-based)
+5. Integrate into `BandAccessor`
+
+### Medium-term
+1. Add PyJulia support (optional dependency)
+2. Implement `PardisoSolver._solve_via_pyjulia()`
+3. Add solver benchmarks
+4. Documentation and examples
+
+### Long-term
+1. Add more solvers (cuSOLVER for GPU, ELPA for HPC)
+2. Auto-selection based on system size
+3. Distributed computing support
+
+## Design Principles
+
+1. **Separation of concerns**: I/O, solving, post-processing are separate
+2. **Pluggable architecture**: Easy to add new solvers
+3. **Backward compatibility**: Default behavior unchanged
+4. **Performance**: Minimize data copying and file I/O
+5. **Testability**: Each component can be unit tested
+6. **Documentation**: Clear examples and API docs
+
+## Example Usage
+
+```python
+# Current workflow (after Phase 1)
+tbsys = TBSystem(data, calculator)
+tbsys.to_pardiso_json("pardiso_input")
+# Run `dptb pdso` or call dptb/postprocess/pardiso/main.jl manually.
+# Band results are written to bandstructure.h5/bands.dat.
+
+# After Phase 2
+tbsys = TBSystem(data, calculator)
+bands = tbsys.band.compute(kpath, solver='pardiso')
+bands.plot()
+
+# After Phase 3
+dptb.set_backend('pardiso')
+tbsys = TBSystem(data, calculator)
+bands = tbsys.band.compute(kpath)  # Automatically fast!
+```
+"""