added chem subpackage

pyproject.toml

@@ -38,6 +38,12 @@ dependencies = [
     "tensornetwork>=0.4.6",
     "torch>=2.7.1",
     "matplotlib",
+    "pandas",
+    "pyscf",
+    "openfermion",
+    "pylatexenc",
+    "noisyopt",
+    "renormalizer==0.0.11",
 ]
 
 [project.optional-dependencies]

src/tyxonq/chem/__init__.py

@@ -0,0 +1,65 @@
+__version__ = "0.1.0"
+__author__ = "TyxonQ"
+
+# Forked from https://github.com/tencent-quantum-lab/TenCirChem 
+# Reconstruction IS IN PROGRESS
+
+
+#  Copyright (c) 2023. The TenCirChem Developers. All Rights Reserved.
+#
+#  This file is distributed under ACADEMIC PUBLIC LICENSE
+#  and WITHOUT ANY WARRANTY. See the LICENSE file for details.
+
+import os
+import logging
+
+os.environ["JAX_ENABLE_X64"] = "True"
+# for debugging
+# os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
+
+# disable CUDA 11.1 warning
+os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"
+
+logger = logging.getLogger("tyxonq")
+logger.setLevel(logging.FATAL)
+
+os.environ["RENO_LOG_LEVEL"] = "100"
+
+logger = logging.getLogger("tencirchem")
+logger.setLevel(logging.WARNING)
+
+# finish logger stuff
+del logger
+
+from .utils.backend import set_backend, set_dtype
+
+# by default use float64 rather than float32
+set_dtype("complex128")
+
+# static module
+# as an external interface
+from pyscf import M
+
+from .static.ucc import UCC
+from .static.uccsd import UCCSD, ROUCCSD
+from .static.kupccgsd import KUPCCGSD
+from .static.puccd import PUCCD
+from .static.hea import HEA, parity, binary, get_noise_conf
+
+# dynamic module
+# as an external interface
+from renormalizer import Op, BasisSHO, BasisHalfSpin, BasisSimpleElectron, BasisMultiElectron, Model, Mpo
+from renormalizer.model import OpSum
+
+from .utils.misc import get_dense_operator
+from .dynamic.time_evolution import TimeEvolution
+
+
+def clear_cache():
+    from .utils.backend import ALL_JIT_LIBS
+    from .static.evolve_civector import CI_OPERATOR_CACHE, CI_OPERATOR_BATCH_CACHE
+
+    for l in ALL_JIT_LIBS:
+        l.clear()
+    CI_OPERATOR_CACHE.clear()
+    CI_OPERATOR_BATCH_CACHE.clear()

src/tyxonq/chem/constants.py

@@ -0,0 +1,20 @@
+#  Copyright (c) 2023. The TenCirChem Developers. All Rights Reserved.
+#
+#  This file is distributed under ACADEMIC PUBLIC LICENSE
+#  and WITHOUT ANY WARRANTY. See the LICENSE file for details.
+
+import numpy as np
+
+a = np.array([[0, 1], [0, 0]])
+# a^\dagger
+ad = a.T
+# a^\dagger_i a_j
+ad_a = np.kron(ad, a)
+# a^\dagger_i a_j - a_j a^\dagger_i
+ad_a_hc = ad_a - ad_a.T
+adad_aa = np.kron(np.kron(np.kron(ad, ad), a), a)
+adad_aa_hc = adad_aa - adad_aa.T
+ad_a_hc2 = ad_a_hc @ ad_a_hc
+adad_aa_hc2 = adad_aa_hc @ adad_aa_hc
+
+DISCARD_EPS = 1e-12

src/tyxonq/chem/dynamic/__init__.py

@@ -0,0 +1,10 @@
+#  Copyright (c) 2023. The TenCirChem Developers. All Rights Reserved.
+#
+#  This file is distributed under ACADEMIC PUBLIC LICENSE
+#  and WITHOUT ANY WARRANTY. See the LICENSE file for details.
+
+
+from ..dynamic.model import sbm, pyrazine
+from ..dynamic.transform import qubit_encode_op, qubit_encode_op_grouped, qubit_encode_basis
+from ..dynamic.time_derivative import get_ansatz, get_jacobian_func, get_deriv, regularized_inversion
+from ..dynamic.time_evolution import TimeEvolution

src/tyxonq/chem/dynamic/model/__init__.py

@@ -0,0 +1,4 @@
+#  Copyright (c) 2023. The TenCirChem Developers. All Rights Reserved.
+#
+#  This file is distributed under ACADEMIC PUBLIC LICENSE
+#  and WITHOUT ANY WARRANTY. See the LICENSE file for details.

src/tyxonq/chem/dynamic/model/pyrazine.py

@@ -0,0 +1,140 @@
+#  Copyright (c) 2023. The TenCirChem Developers. All Rights Reserved.
+#
+#  This file is distributed under ACADEMIC PUBLIC LICENSE
+#  and WITHOUT ANY WARRANTY. See the LICENSE file for details.
+
+
+import logging
+from itertools import permutations as permut
+from itertools import product
+
+
+import numpy as np
+from renormalizer.model import Op
+from renormalizer.model.basis import BasisSHO, BasisMultiElectron
+from renormalizer.utils.constant import ev2au
+
+
+logger = logging.getLogger(__name__)
+
+
+r"""
+Bi-linear vibronic coupling model for Pyrazine, 4 modes
+See: Raab, Worth, Meyer, Cederbaum.  J.Chem.Phys. 110 (1999) 936
+The parameters are from heidelberg mctdh package pyr4+.op
+"""
+
+
+# frequencies
+w10a = 0.1139 * ev2au
+w6a = 0.0739 * ev2au
+w1 = 0.1258 * ev2au
+w9a = 0.1525 * ev2au
+
+# energy-gap
+delta = 0.42300 * ev2au
+
+# linear, on-diagonal coupling coefficients
+# H(1,1)
+_6a_s1_s1 = 0.09806 * ev2au
+_1_s1_s1 = 0.05033 * ev2au
+_9a_s1_s1 = 0.14521 * ev2au
+# H(2,2)
+_6a_s2_s2 = -0.13545 * ev2au
+_1_s2_s2 = 0.17100 * ev2au
+_9a_s2_s2 = 0.03746 * ev2au
+
+# quadratic, on-diagonal coupling coefficients
+# H(1,1)
+_10a_10a_s1_s1 = -0.01159 * ev2au
+_6a_6a_s1_s1 = 0.00000 * ev2au
+_1_1_s1_s1 = 0.00000 * ev2au
+_9a_9a_s1_s1 = 0.00000 * ev2au
+# H(2,2)
+_10a_10a_s2_s2 = -0.01159 * ev2au
+_6a_6a_s2_s2 = 0.00000 * ev2au
+_1_1_s2_s2 = 0.00000 * ev2au
+_9a_9a_s2_s2 = 0.00000 * ev2au
+
+# bilinear, on-diagonal coupling coefficients
+# H(1,1)
+_6a_1_s1_s1 = 0.00108 * ev2au
+_1_9a_s1_s1 = -0.00474 * ev2au
+_6a_9a_s1_s1 = 0.00204 * ev2au
+# H(2,2)
+_6a_1_s2_s2 = -0.00298 * ev2au
+_1_9a_s2_s2 = -0.00155 * ev2au
+_6a_9a_s2_s2 = 0.00189 * ev2au
+
+# linear, off-diagonal coupling coefficients
+_10a_s1_s2 = 0.20804 * ev2au
+
+# bilinear, off-diagonal coupling coefficients
+# H(1,2) and H(2,1)
+_1_10a_s1_s2 = 0.00553 * ev2au
+_6a_10a_s1_s2 = 0.01000 * ev2au
+_9a_10a_s1_s2 = 0.00126 * ev2au
+
+e_list = ["s1", "s2"]
+v_list = ["10a", "6a", "9a", "1"]
+
+
+def get_ham_terms():
+    ham_terms = []
+    for e_isymbol, e_jsymbol in product(e_list, repeat=2):
+        e_op = Op(r"a^\dagger a", [e_isymbol, e_jsymbol])
+        for v_isymbol, v_jsymbol in product(v_list, repeat=2):
+            # linear
+            if v_isymbol == v_jsymbol:
+                # if one of the permutation is defined, then the `e_idx_tuple` term should
+                # be defined as required by Hermitian Hamiltonian
+                for eterm1, eterm2 in permut([e_isymbol, e_jsymbol], 2):
+                    factor = globals().get(f"_{v_isymbol}_{eterm1}_{eterm2}")
+                    if factor is not None:
+                        factor *= np.sqrt(eval(f"w{v_isymbol}"))
+                        ham_terms.append(e_op * Op("x", v_isymbol) * factor)
+                        logger.debug(f"term: {v_isymbol}_{e_isymbol}_{e_jsymbol}")
+                        break
+                else:
+                    logger.debug(f"no term: {v_isymbol}_{e_isymbol}_{e_jsymbol}")
+
+            # quadratic
+            # use product to guarantee `break` breaks the whole loop
+            it = product(permut([v_isymbol, v_jsymbol], 2), permut([e_isymbol, e_jsymbol], 2))
+            for (vterm1, vterm2), (eterm1, eterm2) in it:
+                factor = globals().get(f"_{vterm1}_{vterm2}_{eterm1}_{eterm2}")
+
+                if factor is not None:
+                    factor *= np.sqrt(eval(f"w{v_isymbol}") * eval(f"w{v_jsymbol}"))
+                    if vterm1 == vterm2:
+                        v_op = Op("x^2", vterm1)
+                    else:
+                        v_op = Op("x", vterm1) * Op("x", vterm2)
+                    ham_terms.append(e_op * v_op * factor)
+                    logger.debug(f"term: {v_isymbol}_{v_jsymbol}_{e_isymbol}_{e_jsymbol}")
+                    break
+            else:
+                logger.debug(f"no term: {v_isymbol}_{v_jsymbol}_{e_isymbol}_{e_jsymbol}")
+
+    # electronic coupling
+    ham_terms.append(Op(r"a^\dagger a", "s1", -delta, [0, 0]))
+    ham_terms.append(Op(r"a^\dagger a", "s2", delta, [0, 0]))
+
+    # vibrational kinetic and potential
+    for v_isymbol in v_list:
+        ham_terms.extend([Op("p^2", v_isymbol, 0.5), Op("x^2", v_isymbol, 0.5 * eval("w" + v_isymbol) ** 2)])
+
+    ham_terms = [term for term in ham_terms if term.factor != 0]
+
+    return ham_terms
+
+
+def get_basis(nlevels):
+    if isinstance(nlevels, int):
+        nlevels = [nlevels] * len(v_list)
+    if not isinstance(nlevels, list):
+        raise TypeError(f"`nlevels` must be int or list. Got {type(nlevels)}")
+    basis = [BasisMultiElectron(e_list, [0, 0])]
+    for v_isymbol, n in zip(v_list, nlevels):
+        basis.append(BasisSHO(v_isymbol, globals()[f"w{v_isymbol}"], n))
+    return basis

src/tyxonq/chem/dynamic/model/sbm.py

@@ -0,0 +1,27 @@
+#  Copyright (c) 2023. The TenCirChem Developers. All Rights Reserved.
+#
+#  This file is distributed under ACADEMIC PUBLIC LICENSE
+#  and WITHOUT ANY WARRANTY. See the LICENSE file for details.
+
+
+from renormalizer import BasisHalfSpin, BasisSHO, Op
+
+
+def get_ham_terms(epsilon, delta, nmode, omega_list, g_list):
+    terms = [Op("sigma_z", "spin", epsilon), Op("sigma_x", "spin", delta)]
+    for i in range(nmode):
+        g = g_list[i]
+        omega = omega_list[i]
+        terms.extend([Op(r"b^\dagger b", f"v{i}", omega), Op("sigma_z", "spin", g) * Op(r"b^\dagger+b", f"v{i}")])
+    return terms
+
+
+def get_basis(omega_list, nlevels):
+    if isinstance(nlevels, int):
+        nlevels = [nlevels] * len(omega_list)
+    if not isinstance(nlevels, list):
+        raise TypeError(f"`nlevels` must be int or list. Got {type(nlevels)}")
+    basis = [BasisHalfSpin("spin")]
+    for i in range(len(omega_list)):
+        basis.append(BasisSHO(f"v{i}", omega=omega_list[i], nbas=nlevels[i]))
+    return basis