From c96bda246115d8d52276efa37c3a7f131194915e Mon Sep 17 00:00:00 2001
From: Hugh <keminghu@outlook.com>
Date: Fri, 24 Apr 2026 15:38:42 +0800
Subject: [PATCH 1/6] add class FermionMapBoson

---
 src/deepquantum/photonic/ansatz.py | 308 +++++++++++++++++++++++++++++
 1 file changed, 308 insertions(+)

diff --git a/src/deepquantum/photonic/ansatz.py b/src/deepquantum/photonic/ansatz.py
index 4482f054..59755eb6 100644
--- a/src/deepquantum/photonic/ansatz.py
+++ b/src/deepquantum/photonic/ansatz.py
@@ -1,6 +1,7 @@
 """Ansatze: various photonic quantum circuits"""
 
 import copy
+from itertools import combinations
 from typing import Any
 
 import networkx as nx
@@ -227,3 +228,310 @@ def graph_density(graph: nx.Graph, samples: dict) -> dict:
             samples_[key] = [temp_prob, density]
         sort_samples = sort_dict_fock_basis(samples_, 1)
         return sort_samples
+
+
+class FermionMapBoson:
+    """A class to map molecular Fermionic Hamiltonians to Bosonic Qumode representations.
+
+    This class utilizes OpenFermion and PySCF to perform electronic structure
+    calculations and applies the Dhar-Mandal-Suryanarayana (DMS) mapping to
+    transform the second-quantized Hamiltonian into a Bosonic Fock space.
+    It supports Active Space approximations to reduce the simulation dimensionality.
+
+    Args:
+        config (dict): A configuration dictionary containing the following keys:
+            'geometry' (list): Molecular structure, e.g., [('H', (0,0,0)), ('H', (0,0,0.74))].
+            'basis' (str): Quantum chemistry basis set (e.g., 'sto-3g', '6-31g').
+            'multiplicity' (int): Spin multiplicity (2S + 1). Usually 1 for closed-shell.
+            'charge' (int): Net charge of the molecule.
+            'n_ele' (int): Total number of electrons in the full system.
+            'n_orbit' (int): Total number of spin-orbitals in the full space.
+            'occupied_indices' (list): Indices of spatial orbitals to be frozen (occupied).
+                                        Electrons in these orbitals do not participate in
+                                        excitations but contribute to the energy constant.
+            'active_indices' (list): Indices of spatial orbitals to be treated as active.
+                                       These form the primary Hilbert space for VQE.
+    """
+
+    def __init__(self, config: dict = None) -> None:
+        self.geometry = config['geometry']
+        self.basis = config['basis']
+        self.multiplicity = config['multiplicity']
+        self.charge = config['charge']
+        self.occupied_indices = config['occupied_indices']
+        self.active_indices = config['active_indices']
+        self.n = config['n_ele'] - 2 * len(self.occupied_indices)
+        self.m = 2 * len(self.active_indices)
+
+    def construct_h_fermion(self):
+        from openfermion import get_fermion_operator
+        from openfermion.chem import MolecularData
+        from openfermionpyscf import run_pyscf
+
+        molecule = MolecularData(self.geometry, self.basis, self.multiplicity, self.charge)
+        molecule = run_pyscf(molecule, run_scf=True, run_fci=True)
+        hamiltonian = molecule.get_molecular_hamiltonian(
+            occupied_indices=self.occupied_indices, active_indices=self.active_indices
+        )
+        fermion_op = get_fermion_operator(hamiltonian)
+        h_matrix, basis_f = self.compute_hamiltonian_matrix(fermion_op, self.n, self.m)
+        self.molecule = molecule
+        self.hamiltonian = hamiltonian
+        self.constant = hamiltonian.constant
+        self.fermion_op = fermion_op
+        self.basis_f = basis_f
+        self.h_matrix = h_matrix
+        return h_matrix
+
+    def mapping(self):
+        self.h_fock, self.map_dic = self.get_dms_mapping(self.h_matrix, self.n, self.m)
+        return self.h_fock
+
+    def fci_energy(self):
+        return self.molecule.fci_energy
+
+    @staticmethod
+    def apply_annihilation(state, k):
+        """湮灭算符 f_k 作用在 Slater 行列式上
+
+        输入:
+            state: tuple，有序轨道列表 (p1,...,pN)，p1<p2<...<pN
+            k:     int，要湮灭的轨道
+
+        输出:
+            (new_state, sign) 或 None（结果为零）
+
+            new_state: tuple，(N-1)粒子态
+            sign:      +1 或 -1
+        """
+        state = list(state)
+
+        if k not in state:
+            return None  # 轨道未被占据，结果为零
+
+        m = state.index(k)  # k 在列表中的位置（0-indexed）
+        sign = (-1) ** m  # 费米子符号
+        new_state = tuple(state[:m] + state[m + 1 :])
+
+        return (new_state, sign)
+
+    @staticmethod
+    def apply_creation(state, k):
+        """产生算符 f^†_k 作用在 Slater 行列式上
+
+        输入:
+            state: tuple，有序轨道列表
+            k:     int，要产生的轨道
+
+        输出:
+            (new_state, sign) 或 None（结果为零）
+        """
+        state = list(state)
+
+        if k in state:
+            return None  # 轨道已被占据，Pauli不相容
+
+        # 找到插入位置（保持有序）
+        n = sum(1 for p in state if p < k)  # k 前面有 n 个轨道
+        sign = (-1) ** n
+        new_state = tuple(sorted(state + [k]))
+
+        return (new_state, sign)
+
+    def matrix_element_one_body(self, bra, ket, p, q):
+        """计算单体算符矩阵元 ⟨bra|f^†_p f_q|ket⟩
+
+        输入:
+            bra, ket: tuple，N粒子Slater行列式
+            p, q:     int，轨道指标
+
+        输出:
+            矩阵元值（实数或复数）
+        """
+        # Step 1: f_q 作用在 ket 上
+        result = self.apply_annihilation(ket, q)
+        if result is None:
+            return 0.0
+        int1, sign1 = result
+
+        # Step 2: f†_p 作用在中间态上
+        result = self.apply_creation(int1, p)
+        if result is None:
+            return 0.0
+        final_state, sign2 = result
+
+        # Step 3: 计算内积
+        if final_state == bra:
+            return float(sign1 * sign2)
+        else:
+            return 0.0
+
+    def matrix_element_two_body(self, bra, ket, p, q, r, s):
+        """计算双体算符矩阵元 ⟨bra|f†_p f†_q f_r f_s|ket⟩
+
+        注意算符作用顺序：从右到左，先 f_s，再 f_r，再 f†_q，再 f†_p
+
+        输入:
+            bra, ket: tuple，N粒子Slater行列式
+            p, q, r, s: int，轨道指标
+
+        输出:
+            矩阵元值
+        """
+        # Step 1: f_s 作用在 ket 上
+        result = self.apply_annihilation(ket, s)
+        if result is None:
+            return 0.0
+        int1, sign_s = result
+
+        # Step 2: f_r 作用在 int1 上
+        result = self.apply_annihilation(int1, r)
+        if result is None:
+            return 0.0
+        int2, sign_r = result
+
+        # Step 3: f†_q 作用在 int2 上
+        result = self.apply_creation(int2, q)
+        if result is None:
+            return 0.0
+        int3, sign_q = result
+
+        # Step 4: f†_p 作用在 int3 上
+        result = self.apply_creation(int3, p)
+        if result is None:
+            return 0.0
+        final_state, sign_p = result
+
+        # Step 5: 计算内积
+        if final_state == bra:
+            return float(sign_p * sign_q * sign_r * sign_s)
+        else:
+            return 0.0
+
+    @staticmethod
+    def extract_integrals(fermion_op):
+        """从 FermionOperator 中提取单体积分 h[p,q] 和双体积分 v[p,q,r,s]
+
+        OpenFermion 的 FermionOperator 格式：
+            FermionOperator('0^ 1', 0.5)  → 0.5 * f†_0 f_1
+            FermionOperator('0^ 1^ 2 3', 0.25) → 0.25 * f†_0 f†_1 f_2 f_3
+
+        输出:
+            h: dict，{(p,q): coefficient}，单体项
+            v: dict，{(p,q,r,s): coefficient}，双体项
+            constant: float，常数项（核排斥能等）
+        """
+        h = {}  # 单体积分
+        v = {}  # 双体积分
+        constant = 0.0
+
+        for term, coeff in fermion_op.terms.items():
+            if len(term) == 0:
+                # 常数项
+                constant += coeff.real
+
+            elif len(term) == 2:
+                # 单体项 f†_p f_q
+                # term = ((p, 1), (q, 0))，其中1=产生，0=湮灭
+                # 确保顺序是：产生算符在前
+
+                # 找产生和湮灭算符
+                creators = [idx for idx, dag in term if dag == 1]
+                annihilators = [idx for idx, dag in term if dag == 0]
+
+                if len(creators) == 1 and len(annihilators) == 1:
+                    p, q = creators[0], annihilators[0]
+                    h[(p, q)] = h.get((p, q), 0.0) + coeff.real
+
+            elif len(term) == 4:
+                # 双体项 f†_p f†_q f_r f_s
+                creators = [idx for idx, dag in term if dag == 1]
+                annihilators = [idx for idx, dag in term if dag == 0]
+
+                if len(creators) == 2 and len(annihilators) == 2:
+                    p, q = creators[0], creators[1]
+                    # 注意：湮灭算符顺序需要与算符定义一致
+                    # OpenFermion 中 f†_p f†_q f_r f_s 的 annihilators 顺序是 r, s
+                    r, s = annihilators[0], annihilators[1]
+                    v[(p, q, r, s)] = v.get((p, q, r, s), 0.0) + coeff.real
+
+        return h, v, constant
+
+    def compute_hamiltonian_matrix(self, fermion_op, n, m):
+        """直接计算费米子哈密顿量在 n 粒子子空间的矩阵
+
+        输入:
+            fermion_op: OpenFermion FermionOperator
+            n: 电子数
+            m: 自旋轨道数
+
+        输出:
+            h_matrix: np.ndarray，dim × dim 的厄米矩阵
+            basis:    list of tuples，基矢列表
+        """
+        # Step 1: 枚举基矢
+        basis = list(combinations(range(m), m))
+        dim = len(basis)
+
+        # Step 2: 提取积分系数
+        h_integrals, v_integrals, constant = self.extract_integrals(fermion_op)
+
+        # Step 3: 构造矩阵
+        h_matrix = np.zeros((dim, dim), dtype=complex)
+
+        # 常数项：贡献到所有对角元
+        # for i in range(dim):
+        #     h_matrix[i, i] += constant
+
+        # 单体项
+        for (p, q), h_pq in h_integrals.items():
+            if abs(h_pq) < 1e-8:
+                continue
+
+            for j, ket in enumerate(basis):
+                for i, bra in enumerate(basis):
+                    me = self.matrix_element_one_body(bra, ket, p, q)
+                    if me != 0.0:
+                        h_matrix[i, j] += h_pq * me
+
+        # 双体项
+        for (p, q, r, s), v_pqrs in v_integrals.items():
+            if abs(v_pqrs) < 1e-8:
+                continue
+
+            for j, ket in enumerate(basis):
+                for i, bra in enumerate(basis):
+                    me = self.matrix_element_two_body(bra, ket, p, q, r, s)
+                    if me != 0.0:
+                        # h_matrix[i, j] += 0.5 * v_pqrs * me
+                        h_matrix[i, j] += 1 * v_pqrs * me
+
+        return h_matrix.real, basis
+
+    @staticmethod
+    def get_dms_mapping(h_f, n, m):
+        # 1. 生成所有费米子组合 (按升序)
+        f_basis = list(combinations(range(m), n))
+
+        mapping = {}
+        for i, p in enumerate(f_basis):
+            # 注意：p 是有序元组，如 (0, 1)
+            # p[0] 是 p1, p[1] 是 p2...
+            q = [0] * n
+            # 应用公式 (注意索引对应关系)
+            q[n - 1] = p[0]  # q_N = p1
+            for j in range(1, n):  # 处理 q1 到 q_{N-1}
+                # j 从 1 到 N-1 对应公式中的下标
+                # p 数组索引从 0 开始，所以 p_{N-j+1} 对应 p[N-j], p_{N-j} 对应 p[N-j-1]
+                q[j - 1] = p[n - j] - p[n - j - 1] - 1
+
+            mapping[i] = ['F:', p, 'B:', tuple(q)]
+
+        idx = h_f.nonzero()
+        fock_mapping = []
+        for i in range(len(idx[0])):
+            k = idx[0][i]
+            j = idx[1][i]
+            temp = (h_f[k, j].item(), mapping[k][-1], mapping[j][-1])
+            fock_mapping.append(temp)
+        return fock_mapping, mapping

From 08834c639304f2acb5de6a91a38f0774ba07bec8 Mon Sep 17 00:00:00 2001
From: kh <keminghu@outlook.com>
Date: Fri, 24 Apr 2026 16:04:45 +0800
Subject: [PATCH 2/6] update

---
 src/deepquantum/photonic/ansatz.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/deepquantum/photonic/ansatz.py b/src/deepquantum/photonic/ansatz.py
index 59755eb6..e73d4ea7 100644
--- a/src/deepquantum/photonic/ansatz.py
+++ b/src/deepquantum/photonic/ansatz.py
@@ -470,7 +470,7 @@ def compute_hamiltonian_matrix(self, fermion_op, n, m):
             basis:    list of tuples，基矢列表
         """
         # Step 1: 枚举基矢
-        basis = list(combinations(range(m), m))
+        basis = list(combinations(range(m), n))
         dim = len(basis)
 
         # Step 2: 提取积分系数

From 0148f7439845795fee52bd5fc82fd2069b0dae1e Mon Sep 17 00:00:00 2001
From: Hugh <keminghu@outlook.com>
Date: Mon, 27 Apr 2026 13:42:41 +0800
Subject: [PATCH 3/6] update docstring

---
 src/deepquantum/photonic/ansatz.py | 97 ++++++++++++++----------------
 1 file changed, 44 insertions(+), 53 deletions(-)

diff --git a/src/deepquantum/photonic/ansatz.py b/src/deepquantum/photonic/ansatz.py
index 59755eb6..5f8988c8 100644
--- a/src/deepquantum/photonic/ansatz.py
+++ b/src/deepquantum/photonic/ansatz.py
@@ -2,7 +2,7 @@
 
 import copy
 from itertools import combinations
-from typing import Any
+from typing import Any, TYPE_CHECKING
 
 import networkx as nx
 import numpy as np
@@ -14,6 +14,9 @@
 from .qmath import sort_dict_fock_basis, takagi
 from .state import FockState
 
+if TYPE_CHECKING:
+    from openfermion import FermionOperator
+
 
 class Clements(QumodeCircuit):
     """Clements circuit."""
@@ -291,18 +294,13 @@ def fci_energy(self):
         return self.molecule.fci_energy
 
     @staticmethod
-    def apply_annihilation(state, k):
-        """湮灭算符 f_k 作用在 Slater 行列式上
-
-        输入:
-            state: tuple，有序轨道列表 (p1,...,pN)，p1<p2<...<pN
-            k:     int，要湮灭的轨道
+    def apply_annihilation(state: tuple, k: int):
+        """Apply the annihilation operator :math:`f_k` to a Slater determinant.
 
-        输出:
-            (new_state, sign) 或 None（结果为零）
+        Args:
+            state: An ordered list of occupied orbitals (p1, ..., pN) where p1 < p2 < ... < pN.
+            k: The index of the orbital to be annihilated.
 
-            new_state: tuple，(N-1)粒子态
-            sign:      +1 或 -1
         """
         state = list(state)
 
@@ -317,14 +315,12 @@ def apply_annihilation(state, k):
 
     @staticmethod
     def apply_creation(state, k):
-        """产生算符 f^†_k 作用在 Slater 行列式上
+        """Apply the creation operator :math:`f^†_k` to a Slater determinant.
 
-        输入:
-            state: tuple，有序轨道列表
-            k:     int，要产生的轨道
+        Args:
+            state: An ordered list of occupied orbitals (p1, ..., pN) where p1 < p2 < ... < pN.
+            k: The index of the orbital to be created.
 
-        输出:
-            (new_state, sign) 或 None（结果为零）
         """
         state = list(state)
 
@@ -338,15 +334,14 @@ def apply_creation(state, k):
 
         return (new_state, sign)
 
-    def matrix_element_one_body(self, bra, ket, p, q):
-        """计算单体算符矩阵元 ⟨bra|f^†_p f_q|ket⟩
-
-        输入:
-            bra, ket: tuple，N粒子Slater行列式
-            p, q:     int，轨道指标
+    def matrix_element_one_body(self, bra: tuple, ket: tuple, p: int, q: int):
+        """Calculate the matrix element of a one-body operator: :math:`⟨bra|f^†_p f_q|ket⟩`.
 
-        输出:
-            矩阵元值（实数或复数）
+        Args:
+            bra: N-particle Slater determinant (bra state).
+            ket: N-particle Slater determinant (ket state).
+            p: Index of the creation orbital.
+            q: Index of the annihilation orbital.
         """
         # Step 1: f_q 作用在 ket 上
         result = self.apply_annihilation(ket, q)
@@ -366,17 +361,19 @@ def matrix_element_one_body(self, bra, ket, p, q):
         else:
             return 0.0
 
-    def matrix_element_two_body(self, bra, ket, p, q, r, s):
-        """计算双体算符矩阵元 ⟨bra|f†_p f†_q f_r f_s|ket⟩
+    def matrix_element_two_body(self, bra: tuple, ket: tuple, p: int, q: int, r: int, s: int):
+        """Calculate the matrix element of a two-body operator: :math:`⟨bra|f†_p f†_q f_r f_s|ket⟩`.
 
-        注意算符作用顺序：从右到左，先 f_s，再 f_r，再 f†_q，再 f†_p
+        The operators are applied from right to left: :math:`f_s`, then :math:`f_r`,
+        then :math:`f†_q`, then :math:`f†_p`.
 
-        输入:
-            bra, ket: tuple，N粒子Slater行列式
-            p, q, r, s: int，轨道指标
-
-        输出:
-            矩阵元值
+        Args:
+            bra: N-particle Slater determinant (bra state).
+            ket: N-particle Slater determinant (ket state).
+            p: Indices for the creation operators.
+            q: Indices for the creation operators.
+            r: Indices for the annihilation operators.
+            s: Indices for the annihilation operators.
         """
         # Step 1: f_s 作用在 ket 上
         result = self.apply_annihilation(ket, s)
@@ -409,17 +406,15 @@ def matrix_element_two_body(self, bra, ket, p, q, r, s):
             return 0.0
 
     @staticmethod
-    def extract_integrals(fermion_op):
-        """从 FermionOperator 中提取单体积分 h[p,q] 和双体积分 v[p,q,r,s]
+    def extract_integrals(fermion_op: 'FermionOperator'):
+        """Extract one-body integrals :math:`h[p,q]` and two-body integrals :math:`v[p,q,r,s]` from aFermionOperator.
 
-        OpenFermion 的 FermionOperator 格式：
-            FermionOperator('0^ 1', 0.5)  → 0.5 * f†_0 f_1
-            FermionOperator('0^ 1^ 2 3', 0.25) → 0.25 * f†_0 f†_1 f_2 f_3
+        The FermionOperator format handles terms like:
+            FermionOperator(0.5, '0^ 1')       -> 0.5 * f†_0 f_1
+            FermionOperator(0.25, '0^ 1^ 2 3') -> 0.25 * f†_0 f†_1 f_2 f_3
 
-        输出:
-            h: dict，{(p,q): coefficient}，单体项
-            v: dict，{(p,q,r,s): coefficient}，双体项
-            constant: float，常数项（核排斥能等）
+        Args:
+            fermion_op: The operator object from OpenFermion.
         """
         h = {}  # 单体积分
         v = {}  # 双体积分
@@ -457,17 +452,13 @@ def extract_integrals(fermion_op):
 
         return h, v, constant
 
-    def compute_hamiltonian_matrix(self, fermion_op, n, m):
-        """直接计算费米子哈密顿量在 n 粒子子空间的矩阵
-
-        输入:
-            fermion_op: OpenFermion FermionOperator
-            n: 电子数
-            m: 自旋轨道数
+    def compute_hamiltonian_matrix(self, fermion_op: 'FermionOperator', n: int, m: int):
+        """Directly compute the Hamiltonian matrix in the n-particle subspace.
 
-        输出:
-            h_matrix: np.ndarray，dim × dim 的厄米矩阵
-            basis:    list of tuples，基矢列表
+        Args:
+            fermion_op: The input OpenFermion Hamiltonian.
+            n: Number of particles (electrons).
+            m: Number of spin-orbitals.
         """
         # Step 1: 枚举基矢
         basis = list(combinations(range(m), m))

From 441a2f1ec990bff5dafdf7dbb211f8e762198f84 Mon Sep 17 00:00:00 2001
From: Hugh <keminghu@outlook.com>
Date: Mon, 27 Apr 2026 14:01:43 +0800
Subject: [PATCH 4/6] update

---
 src/deepquantum/photonic/ansatz.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/deepquantum/photonic/ansatz.py b/src/deepquantum/photonic/ansatz.py
index f6318753..b684390c 100644
--- a/src/deepquantum/photonic/ansatz.py
+++ b/src/deepquantum/photonic/ansatz.py
@@ -298,7 +298,7 @@ def apply_annihilation(state: tuple, k: int):
         """Apply the annihilation operator :math:`f_k` to a Slater determinant.
 
         Args:
-            state: An ordered list of occupied orbitals (p1, ..., pN) where p1 < p2 < ... < pN.
+            state: An ordered list of occupied orbitals :math:`(p1, ..., pN)` where :math:`p1 < p2 < ... < pN`.
             k: The index of the orbital to be annihilated.
 
         """
@@ -318,7 +318,7 @@ def apply_creation(state, k):
         """Apply the creation operator :math:`f^†_k` to a Slater determinant.
 
         Args:
-            state: An ordered list of occupied orbitals (p1, ..., pN) where p1 < p2 < ... < pN.
+            state: An ordered list of occupied orbitals :math:`(p1, ..., pN)` where :math:`p1 < p2 < ... < pN`.
             k: The index of the orbital to be created.
 
         """
@@ -407,7 +407,7 @@ def matrix_element_two_body(self, bra: tuple, ket: tuple, p: int, q: int, r: int
 
     @staticmethod
     def extract_integrals(fermion_op: 'FermionOperator'):
-        """Extract one-body integrals :math:`h[p,q]` and two-body integrals :math:`v[p,q,r,s]` from aFermionOperator.
+        """Extract one-body integrals :math:`h[p,q]` and two-body integrals :math:`v[p,q,r,s]` from a FermionOperator.
 
         The FermionOperator format handles terms like:
             FermionOperator(0.5, '0^ 1')       -> 0.5 * f†_0 f_1

From 4936a661263e123946b41f869327cce37202f441 Mon Sep 17 00:00:00 2001
From: Hugh <keminghu@outlook.com>
Date: Tue, 26 May 2026 14:37:23 +0800
Subject: [PATCH 5/6] move FermoinMapBoson class to qchem.py

---
 src/deepquantum/photonic/ansatz.py | 301 +----------------------------
 src/deepquantum/photonic/qchem.py  | 298 ++++++++++++++++++++++++++++
 2 files changed, 299 insertions(+), 300 deletions(-)
 create mode 100644 src/deepquantum/photonic/qchem.py

diff --git a/src/deepquantum/photonic/ansatz.py b/src/deepquantum/photonic/ansatz.py
index b684390c..4482f054 100644
--- a/src/deepquantum/photonic/ansatz.py
+++ b/src/deepquantum/photonic/ansatz.py
@@ -1,8 +1,7 @@
 """Ansatze: various photonic quantum circuits"""
 
 import copy
-from itertools import combinations
-from typing import Any, TYPE_CHECKING
+from typing import Any
 
 import networkx as nx
 import numpy as np
@@ -14,9 +13,6 @@
 from .qmath import sort_dict_fock_basis, takagi
 from .state import FockState
 
-if TYPE_CHECKING:
-    from openfermion import FermionOperator
-
 
 class Clements(QumodeCircuit):
     """Clements circuit."""
@@ -231,298 +227,3 @@ def graph_density(graph: nx.Graph, samples: dict) -> dict:
             samples_[key] = [temp_prob, density]
         sort_samples = sort_dict_fock_basis(samples_, 1)
         return sort_samples
-
-
-class FermionMapBoson:
-    """A class to map molecular Fermionic Hamiltonians to Bosonic Qumode representations.
-
-    This class utilizes OpenFermion and PySCF to perform electronic structure
-    calculations and applies the Dhar-Mandal-Suryanarayana (DMS) mapping to
-    transform the second-quantized Hamiltonian into a Bosonic Fock space.
-    It supports Active Space approximations to reduce the simulation dimensionality.
-
-    Args:
-        config (dict): A configuration dictionary containing the following keys:
-            'geometry' (list): Molecular structure, e.g., [('H', (0,0,0)), ('H', (0,0,0.74))].
-            'basis' (str): Quantum chemistry basis set (e.g., 'sto-3g', '6-31g').
-            'multiplicity' (int): Spin multiplicity (2S + 1). Usually 1 for closed-shell.
-            'charge' (int): Net charge of the molecule.
-            'n_ele' (int): Total number of electrons in the full system.
-            'n_orbit' (int): Total number of spin-orbitals in the full space.
-            'occupied_indices' (list): Indices of spatial orbitals to be frozen (occupied).
-                                        Electrons in these orbitals do not participate in
-                                        excitations but contribute to the energy constant.
-            'active_indices' (list): Indices of spatial orbitals to be treated as active.
-                                       These form the primary Hilbert space for VQE.
-    """
-
-    def __init__(self, config: dict = None) -> None:
-        self.geometry = config['geometry']
-        self.basis = config['basis']
-        self.multiplicity = config['multiplicity']
-        self.charge = config['charge']
-        self.occupied_indices = config['occupied_indices']
-        self.active_indices = config['active_indices']
-        self.n = config['n_ele'] - 2 * len(self.occupied_indices)
-        self.m = 2 * len(self.active_indices)
-
-    def construct_h_fermion(self):
-        from openfermion import get_fermion_operator
-        from openfermion.chem import MolecularData
-        from openfermionpyscf import run_pyscf
-
-        molecule = MolecularData(self.geometry, self.basis, self.multiplicity, self.charge)
-        molecule = run_pyscf(molecule, run_scf=True, run_fci=True)
-        hamiltonian = molecule.get_molecular_hamiltonian(
-            occupied_indices=self.occupied_indices, active_indices=self.active_indices
-        )
-        fermion_op = get_fermion_operator(hamiltonian)
-        h_matrix, basis_f = self.compute_hamiltonian_matrix(fermion_op, self.n, self.m)
-        self.molecule = molecule
-        self.hamiltonian = hamiltonian
-        self.constant = hamiltonian.constant
-        self.fermion_op = fermion_op
-        self.basis_f = basis_f
-        self.h_matrix = h_matrix
-        return h_matrix
-
-    def mapping(self):
-        self.h_fock, self.map_dic = self.get_dms_mapping(self.h_matrix, self.n, self.m)
-        return self.h_fock
-
-    def fci_energy(self):
-        return self.molecule.fci_energy
-
-    @staticmethod
-    def apply_annihilation(state: tuple, k: int):
-        """Apply the annihilation operator :math:`f_k` to a Slater determinant.
-
-        Args:
-            state: An ordered list of occupied orbitals :math:`(p1, ..., pN)` where :math:`p1 < p2 < ... < pN`.
-            k: The index of the orbital to be annihilated.
-
-        """
-        state = list(state)
-
-        if k not in state:
-            return None  # 轨道未被占据，结果为零
-
-        m = state.index(k)  # k 在列表中的位置（0-indexed）
-        sign = (-1) ** m  # 费米子符号
-        new_state = tuple(state[:m] + state[m + 1 :])
-
-        return (new_state, sign)
-
-    @staticmethod
-    def apply_creation(state, k):
-        """Apply the creation operator :math:`f^†_k` to a Slater determinant.
-
-        Args:
-            state: An ordered list of occupied orbitals :math:`(p1, ..., pN)` where :math:`p1 < p2 < ... < pN`.
-            k: The index of the orbital to be created.
-
-        """
-        state = list(state)
-
-        if k in state:
-            return None  # 轨道已被占据，Pauli不相容
-
-        # 找到插入位置（保持有序）
-        n = sum(1 for p in state if p < k)  # k 前面有 n 个轨道
-        sign = (-1) ** n
-        new_state = tuple(sorted(state + [k]))
-
-        return (new_state, sign)
-
-    def matrix_element_one_body(self, bra: tuple, ket: tuple, p: int, q: int):
-        """Calculate the matrix element of a one-body operator: :math:`⟨bra|f^†_p f_q|ket⟩`.
-
-        Args:
-            bra: N-particle Slater determinant (bra state).
-            ket: N-particle Slater determinant (ket state).
-            p: Index of the creation orbital.
-            q: Index of the annihilation orbital.
-        """
-        # Step 1: f_q 作用在 ket 上
-        result = self.apply_annihilation(ket, q)
-        if result is None:
-            return 0.0
-        int1, sign1 = result
-
-        # Step 2: f†_p 作用在中间态上
-        result = self.apply_creation(int1, p)
-        if result is None:
-            return 0.0
-        final_state, sign2 = result
-
-        # Step 3: 计算内积
-        if final_state == bra:
-            return float(sign1 * sign2)
-        else:
-            return 0.0
-
-    def matrix_element_two_body(self, bra: tuple, ket: tuple, p: int, q: int, r: int, s: int):
-        """Calculate the matrix element of a two-body operator: :math:`⟨bra|f†_p f†_q f_r f_s|ket⟩`.
-
-        The operators are applied from right to left: :math:`f_s`, then :math:`f_r`,
-        then :math:`f†_q`, then :math:`f†_p`.
-
-        Args:
-            bra: N-particle Slater determinant (bra state).
-            ket: N-particle Slater determinant (ket state).
-            p: Indices for the creation operators.
-            q: Indices for the creation operators.
-            r: Indices for the annihilation operators.
-            s: Indices for the annihilation operators.
-        """
-        # Step 1: f_s 作用在 ket 上
-        result = self.apply_annihilation(ket, s)
-        if result is None:
-            return 0.0
-        int1, sign_s = result
-
-        # Step 2: f_r 作用在 int1 上
-        result = self.apply_annihilation(int1, r)
-        if result is None:
-            return 0.0
-        int2, sign_r = result
-
-        # Step 3: f†_q 作用在 int2 上
-        result = self.apply_creation(int2, q)
-        if result is None:
-            return 0.0
-        int3, sign_q = result
-
-        # Step 4: f†_p 作用在 int3 上
-        result = self.apply_creation(int3, p)
-        if result is None:
-            return 0.0
-        final_state, sign_p = result
-
-        # Step 5: 计算内积
-        if final_state == bra:
-            return float(sign_p * sign_q * sign_r * sign_s)
-        else:
-            return 0.0
-
-    @staticmethod
-    def extract_integrals(fermion_op: 'FermionOperator'):
-        """Extract one-body integrals :math:`h[p,q]` and two-body integrals :math:`v[p,q,r,s]` from a FermionOperator.
-
-        The FermionOperator format handles terms like:
-            FermionOperator(0.5, '0^ 1')       -> 0.5 * f†_0 f_1
-            FermionOperator(0.25, '0^ 1^ 2 3') -> 0.25 * f†_0 f†_1 f_2 f_3
-
-        Args:
-            fermion_op: The operator object from OpenFermion.
-        """
-        h = {}  # 单体积分
-        v = {}  # 双体积分
-        constant = 0.0
-
-        for term, coeff in fermion_op.terms.items():
-            if len(term) == 0:
-                # 常数项
-                constant += coeff.real
-
-            elif len(term) == 2:
-                # 单体项 f†_p f_q
-                # term = ((p, 1), (q, 0))，其中1=产生，0=湮灭
-                # 确保顺序是：产生算符在前
-
-                # 找产生和湮灭算符
-                creators = [idx for idx, dag in term if dag == 1]
-                annihilators = [idx for idx, dag in term if dag == 0]
-
-                if len(creators) == 1 and len(annihilators) == 1:
-                    p, q = creators[0], annihilators[0]
-                    h[(p, q)] = h.get((p, q), 0.0) + coeff.real
-
-            elif len(term) == 4:
-                # 双体项 f†_p f†_q f_r f_s
-                creators = [idx for idx, dag in term if dag == 1]
-                annihilators = [idx for idx, dag in term if dag == 0]
-
-                if len(creators) == 2 and len(annihilators) == 2:
-                    p, q = creators[0], creators[1]
-                    # 注意：湮灭算符顺序需要与算符定义一致
-                    # OpenFermion 中 f†_p f†_q f_r f_s 的 annihilators 顺序是 r, s
-                    r, s = annihilators[0], annihilators[1]
-                    v[(p, q, r, s)] = v.get((p, q, r, s), 0.0) + coeff.real
-
-        return h, v, constant
-
-    def compute_hamiltonian_matrix(self, fermion_op: 'FermionOperator', n: int, m: int):
-        """Directly compute the Hamiltonian matrix in the n-particle subspace.
-
-        Args:
-            fermion_op: The input OpenFermion Hamiltonian.
-            n: Number of particles (electrons).
-            m: Number of spin-orbitals.
-        """
-        # Step 1: 枚举基矢
-        basis = list(combinations(range(m), n))
-        dim = len(basis)
-
-        # Step 2: 提取积分系数
-        h_integrals, v_integrals, constant = self.extract_integrals(fermion_op)
-
-        # Step 3: 构造矩阵
-        h_matrix = np.zeros((dim, dim), dtype=complex)
-
-        # 常数项：贡献到所有对角元
-        # for i in range(dim):
-        #     h_matrix[i, i] += constant
-
-        # 单体项
-        for (p, q), h_pq in h_integrals.items():
-            if abs(h_pq) < 1e-8:
-                continue
-
-            for j, ket in enumerate(basis):
-                for i, bra in enumerate(basis):
-                    me = self.matrix_element_one_body(bra, ket, p, q)
-                    if me != 0.0:
-                        h_matrix[i, j] += h_pq * me
-
-        # 双体项
-        for (p, q, r, s), v_pqrs in v_integrals.items():
-            if abs(v_pqrs) < 1e-8:
-                continue
-
-            for j, ket in enumerate(basis):
-                for i, bra in enumerate(basis):
-                    me = self.matrix_element_two_body(bra, ket, p, q, r, s)
-                    if me != 0.0:
-                        # h_matrix[i, j] += 0.5 * v_pqrs * me
-                        h_matrix[i, j] += 1 * v_pqrs * me
-
-        return h_matrix.real, basis
-
-    @staticmethod
-    def get_dms_mapping(h_f, n, m):
-        # 1. 生成所有费米子组合 (按升序)
-        f_basis = list(combinations(range(m), n))
-
-        mapping = {}
-        for i, p in enumerate(f_basis):
-            # 注意：p 是有序元组，如 (0, 1)
-            # p[0] 是 p1, p[1] 是 p2...
-            q = [0] * n
-            # 应用公式 (注意索引对应关系)
-            q[n - 1] = p[0]  # q_N = p1
-            for j in range(1, n):  # 处理 q1 到 q_{N-1}
-                # j 从 1 到 N-1 对应公式中的下标
-                # p 数组索引从 0 开始，所以 p_{N-j+1} 对应 p[N-j], p_{N-j} 对应 p[N-j-1]
-                q[j - 1] = p[n - j] - p[n - j - 1] - 1
-
-            mapping[i] = ['F:', p, 'B:', tuple(q)]
-
-        idx = h_f.nonzero()
-        fock_mapping = []
-        for i in range(len(idx[0])):
-            k = idx[0][i]
-            j = idx[1][i]
-            temp = (h_f[k, j].item(), mapping[k][-1], mapping[j][-1])
-            fock_mapping.append(temp)
-        return fock_mapping, mapping
diff --git a/src/deepquantum/photonic/qchem.py b/src/deepquantum/photonic/qchem.py
new file mode 100644
index 00000000..d03fc9ee
--- /dev/null
+++ b/src/deepquantum/photonic/qchem.py
@@ -0,0 +1,298 @@
+"""Quantum chemistry helpers for photonic simulations.
+
+This module builds molecular fermionic Hamiltonians with OpenFermion/PySCF and maps
+fixed-particle electronic-structure problems to bosonic Fock-space representations for
+photonic qumode circuits.
+"""
+
+from itertools import combinations
+from typing import Any, TYPE_CHECKING
+import numpy as np
+
+if TYPE_CHECKING:
+    from openfermion import FermionOperator
+
+class FermionMapBoson:
+    """A class to map molecular Fermionic Hamiltonians to Bosonic Qumode representations.
+
+    This class utilizes OpenFermion and PySCF to perform electronic structure
+    calculations and applies the Dhar-Mandal-Suryanarayana (DMS) mapping to
+    transform the second-quantized Hamiltonian into a Bosonic Fock space.
+    It supports Active Space approximations to reduce the simulation dimensionality.
+
+    Args:
+        config (dict): A configuration dictionary containing the following keys:
+            'geometry' (list): Molecular structure, e.g., [('H', (0,0,0)), ('H', (0,0,0.74))].
+            'basis' (str): Quantum chemistry basis set (e.g., 'sto-3g', '6-31g').
+            'multiplicity' (int): Spin multiplicity (2S + 1). Usually 1 for closed-shell.
+            'charge' (int): Net charge of the molecule.
+            'n_ele' (int): Total number of electrons in the full system.
+            'n_orbit' (int): Total number of spin-orbitals in the full space.
+            'occupied_indices' (list): Indices of spatial orbitals to be frozen (occupied).
+                                        Electrons in these orbitals do not participate in
+                                        excitations but contribute to the energy constant.
+            'active_indices' (list): Indices of spatial orbitals to be treated as active.
+                                       These form the primary Hilbert space for VQE.
+    """
+
+    def __init__(self, config: dict = None) -> None:
+        self.geometry = config['geometry']
+        self.basis = config['basis']
+        self.multiplicity = config['multiplicity']
+        self.charge = config['charge']
+        self.occupied_indices = config['occupied_indices']
+        self.active_indices = config['active_indices']
+        self.n = config['n_ele'] - 2 * len(self.occupied_indices)
+        self.m = 2 * len(self.active_indices)
+
+    def construct_h_fermion(self):
+        from openfermion import get_fermion_operator
+        from openfermion.chem import MolecularData
+        from openfermionpyscf import run_pyscf
+
+        molecule = MolecularData(self.geometry, self.basis, self.multiplicity, self.charge)
+        molecule = run_pyscf(molecule, run_scf=True, run_fci=True)
+        hamiltonian = molecule.get_molecular_hamiltonian(
+            occupied_indices=self.occupied_indices, active_indices=self.active_indices
+        )
+        fermion_op = get_fermion_operator(hamiltonian)
+        h_matrix, basis_f = self.compute_hamiltonian_matrix(fermion_op, self.n, self.m)
+        self.molecule = molecule
+        self.hamiltonian = hamiltonian
+        self.constant = hamiltonian.constant
+        self.fermion_op = fermion_op
+        self.basis_f = basis_f
+        self.h_matrix = h_matrix
+        return h_matrix
+
+    def mapping(self):
+        self.h_fock, self.map_dic = self.get_dms_mapping(self.h_matrix, self.n, self.m)
+        return self.h_fock
+
+    def fci_energy(self):
+        return self.molecule.fci_energy
+
+    @staticmethod
+    def apply_annihilation(state: tuple, k: int):
+        """Apply the annihilation operator :math:`f_k` to a Slater determinant.
+
+        Args:
+            state: An ordered list of occupied orbitals :math:`(p1, ..., pN)` where :math:`p1 < p2 < ... < pN`.
+            k: The index of the orbital to be annihilated.
+
+        """
+        state = list(state)
+
+        if k not in state:
+            return None
+
+        m = state.index(k)
+        sign = (-1) ** m
+        new_state = tuple(state[:m] + state[m + 1 :])
+
+        return (new_state, sign)
+
+    @staticmethod
+    def apply_creation(state, k):
+        """Apply the creation operator :math:`f^†_k` to a Slater determinant.
+
+        Args:
+            state: An ordered list of occupied orbitals :math:`(p1, ..., pN)` where :math:`p1 < p2 < ... < pN`.
+            k: The index of the orbital to be created.
+
+        """
+        state = list(state)
+
+        if k in state:
+            return None
+
+        n = sum(1 for p in state if p < k)
+        sign = (-1) ** n
+        new_state = tuple(sorted(state + [k]))
+
+        return (new_state, sign)
+
+    def matrix_element_one_body(self, bra: tuple, ket: tuple, p: int, q: int):
+        """Calculate the matrix element of a one-body operator: :math:`⟨bra|f^†_p f_q|ket⟩`.
+
+        Args:
+            bra: N-particle Slater determinant (bra state).
+            ket: N-particle Slater determinant (ket state).
+            p: Index of the creation orbital.
+            q: Index of the annihilation orbital.
+        """
+        # Step 1: f_q acting on ket
+        result = self.apply_annihilation(ket, q)
+        if result is None:
+            return 0.0
+        int1, sign1 = result
+
+        # Step 2: f†_p acting on int1
+        result = self.apply_creation(int1, p)
+        if result is None:
+            return 0.0
+        final_state, sign2 = result
+
+        # Step 3: inner product
+        if final_state == bra:
+            return float(sign1 * sign2)
+        else:
+            return 0.0
+
+    def matrix_element_two_body(self, bra: tuple, ket: tuple, p: int, q: int, r: int, s: int):
+        """Calculate the matrix element of a two-body operator: :math:`⟨bra|f†_p f†_q f_r f_s|ket⟩`.
+
+        The operators are applied from right to left: :math:`f_s`, then :math:`f_r`,
+        then :math:`f†_q`, then :math:`f†_p`.
+
+        Args:
+            bra: N-particle Slater determinant (bra state).
+            ket: N-particle Slater determinant (ket state).
+            p: Indices for the creation operators.
+            q: Indices for the creation operators.
+            r: Indices for the annihilation operators.
+            s: Indices for the annihilation operators.
+        """
+        # Step 1: f_s acting on ket
+        result = self.apply_annihilation(ket, s)
+        if result is None:
+            return 0.0
+        int1, sign_s = result
+
+        # Step 2: f_r acting on int1
+        result = self.apply_annihilation(int1, r)
+        if result is None:
+            return 0.0
+        int2, sign_r = result
+
+        # Step 3: f†_q acting on int2
+        result = self.apply_creation(int2, q)
+        if result is None:
+            return 0.0
+        int3, sign_q = result
+
+        # Step 4: f†_p acting on int3
+        result = self.apply_creation(int3, p)
+        if result is None:
+            return 0.0
+        final_state, sign_p = result
+
+        # Step 5: inner product
+        if final_state == bra:
+            return float(sign_p * sign_q * sign_r * sign_s)
+        else:
+            return 0.0
+
+    @staticmethod
+    def extract_integrals(fermion_op: 'FermionOperator'):
+        """Extract one-body integrals :math:`h[p,q]` and two-body integrals :math:`v[p,q,r,s]` from a FermionOperator.
+
+        This helper assumes that ``fermion_op`` is already in normal order, with all creation operators placed before
+        annihilation operators. The parsed terms therefore follow :math:`a_p^\dagger a_q` for one-body contributions
+        and :math:`a_p^\dagger a_q^\dagger a_r a_s` for two-body contributions. If the input may contain unordered
+        products, apply OpenFermion's ``normal_ordered`` before calling this method so that fermionic sign changes and
+        constant contractions are handled correctly.
+
+        The FermionOperator format handles terms like:
+            FermionOperator(0.5, '0^ 1')       -> 0.5 * f†_0 f_1
+            FermionOperator(0.25, '0^ 1^ 2 3') -> 0.25 * f†_0 f†_1 f_2 f_3
+
+        Args:
+            fermion_op: The operator object from OpenFermion.
+        """
+        h = {}  # one-body
+        v = {}  # two-body
+        constant = 0.0
+
+        for term, coeff in fermion_op.terms.items():
+            if len(term) == 0:
+                # constant
+                constant += coeff.real
+
+            elif len(term) == 2:
+                # one-body term f†_p f_q
+                # term = ((p, 1), (q, 0))，1=creation，0=annihilation
+
+                creators = [idx for idx, dag in term if dag == 1]
+                annihilators = [idx for idx, dag in term if dag == 0]
+
+                if len(creators) == 1 and len(annihilators) == 1:
+                    p, q = creators[0], annihilators[0]
+                    h[(p, q)] = h.get((p, q), 0.0) + coeff.real
+
+            elif len(term) == 4:
+                # two-body term f†_p f†_q f_r f_s
+                creators = [idx for idx, dag in term if dag == 1]
+                annihilators = [idx for idx, dag in term if dag == 0]
+
+                if len(creators) == 2 and len(annihilators) == 2:
+                    p, q = creators[0], creators[1]
+                    r, s = annihilators[0], annihilators[1]
+                    v[(p, q, r, s)] = v.get((p, q, r, s), 0.0) + coeff.real
+
+        return h, v, constant
+
+    def compute_hamiltonian_matrix(self, fermion_op: 'FermionOperator', n: int, m: int):
+        """Directly compute the Hamiltonian matrix in the n-particle subspace.
+
+        Args:
+            fermion_op: The input OpenFermion Hamiltonian.
+            n: Number of particles (electrons).
+            m: Number of spin-orbitals.
+        """
+        basis = list(combinations(range(m), n))
+        dim = len(basis)
+
+        h_integrals, v_integrals, constant = self.extract_integrals(fermion_op)
+
+        h_matrix = np.zeros((dim, dim), dtype=complex)
+
+        # constant term is not added here
+        # for i in range(dim):
+        #     h_matrix[i, i] += constant
+
+        # one-body term
+        for (p, q), h_pq in h_integrals.items():
+            if abs(h_pq) < 1e-8:
+                continue
+
+            for j, ket in enumerate(basis):
+                for i, bra in enumerate(basis):
+                    me = self.matrix_element_one_body(bra, ket, p, q)
+                    if me != 0.0:
+                        h_matrix[i, j] += h_pq * me
+
+        # two-body term
+        for (p, q, r, s), v_pqrs in v_integrals.items():
+            if abs(v_pqrs) < 1e-8:
+                continue
+
+            for j, ket in enumerate(basis):
+                for i, bra in enumerate(basis):
+                    me = self.matrix_element_two_body(bra, ket, p, q, r, s)
+                    if me != 0.0:
+                        # h_matrix[i, j] += 0.5 * v_pqrs * me
+                        h_matrix[i, j] += 1 * v_pqrs * me
+
+        return h_matrix.real, basis
+
+    @staticmethod
+    def get_dms_mapping(h_f, n, m):
+        f_basis = list(combinations(range(m), n))
+        mapping = {}
+        for i, p in enumerate(f_basis):
+            q = [0] * n
+            q[n - 1] = p[0]  # q_N = p1
+            for j in range(1, n):
+                q[j - 1] = p[n - j] - p[n - j - 1] - 1
+
+            mapping[i] = ['F:', p, 'B:', tuple(q)]
+
+        idx = h_f.nonzero()
+        fock_mapping = []
+        for i in range(len(idx[0])):
+            k = idx[0][i]
+            j = idx[1][i]
+            temp = (h_f[k, j].item(), mapping[k][-1], mapping[j][-1])
+            fock_mapping.append(temp)
+        return fock_mapping, mapping

From cf350d6dd88ccc032ca7cd8b04379f861e9fd4a6 Mon Sep 17 00:00:00 2001
From: Hugh <keminghu@outlook.com>
Date: Tue, 26 May 2026 14:40:55 +0800
Subject: [PATCH 6/6] update

---
 src/deepquantum/photonic/qchem.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/deepquantum/photonic/qchem.py b/src/deepquantum/photonic/qchem.py
index d03fc9ee..e2821a4e 100644
--- a/src/deepquantum/photonic/qchem.py
+++ b/src/deepquantum/photonic/qchem.py
@@ -6,12 +6,14 @@
 """
 
 from itertools import combinations
-from typing import Any, TYPE_CHECKING
+from typing import TYPE_CHECKING
+
 import numpy as np
 
 if TYPE_CHECKING:
     from openfermion import FermionOperator
 
+
 class FermionMapBoson:
     """A class to map molecular Fermionic Hamiltonians to Bosonic Qumode representations.
 
@@ -185,7 +187,7 @@ def matrix_element_two_body(self, bra: tuple, ket: tuple, p: int, q: int, r: int
 
     @staticmethod
     def extract_integrals(fermion_op: 'FermionOperator'):
-        """Extract one-body integrals :math:`h[p,q]` and two-body integrals :math:`v[p,q,r,s]` from a FermionOperator.
+        r"""Extract one-body integrals :math:`h[p,q]` and two-body integrals :math:`v[p,q,r,s]` from a FermionOperator.
 
         This helper assumes that ``fermion_op`` is already in normal order, with all creation operators placed before
         annihilation operators. The parsed terms therefore follow :math:`a_p^\dagger a_q` for one-body contributions