Merge pull request #111 from zpparks314/ibm_ro_buffer_work

Ibm ro buffer work

python/CMakeLists.txt

@@ -48,6 +48,10 @@ else()
    set_target_properties(_pyxacc PROPERTIES LINK_FLAGS "-shared")
 endif()
 
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/benchmark/manage.py.in
+                    ${CMAKE_CURRENT_SOURCE_DIR}/benchmark/manage.py)
+
 install(TARGETS _pyxacc DESTINATION ${CMAKE_INSTALL_PREFIX})
 install(FILES xacc.py DESTINATION ${CMAKE_INSTALL_PREFIX})
 install(DIRECTORY plugin_generator DESTINATION ${CMAKE_INSTALL_PREFIX})
+install(DIRECTORY benchmark DESTINATION ${CMAKE_INSTALL_PREFIX})

python/benchmark/README.md

+# xacc-python-plugin
+
+XACC Python plugin installer for benchmark and decorator algorithms
+
+A simple project and script to install XACC Python plugins to the XACC plugin directory
+
+To install a plugin package:
+```bash
+$ python3 manage.py -i name_of_package
+```
+
+To print the lists of available packages to install:
+```bash
+$ python3 manage.py --list
+```
+
+Adding New Plugin Packages
+--------------------------
+The project is divided into master directories and plugin package directories.
+The plugin packages are found within the master directories.
+
+The master directories must contain an `install.ini` file.
+The `install.ini` files are structured as
+```bash
+name_of_package: /{dir_of_package}
+```
+
+A new master directory can be added to the installer by editing the `MASTER_DIRS` field in `manage.py`.
+
+If the XACC Python plugin directory exists somewhere besides `~/.xacc/py-plugins`, the
+`XACC_PYTHON_PLUGIN_PATH` field in `manage.py` should be edited.
+
+Within the plugin package directories are the XACC plugins (`.py`) that get installed.
+If the plugin package has multiple directories, another `install.ini` file is required in that package and
+must be structured in the same manner as before.
+
+If the plugin package does not have extra directories, nothing else needs to be done.

python/benchmark/manage.py.in

+# Installation management script for XACC benchmark plugin installation
+#
+import sys
+import argparse
+import os
+import configparser
+from shutil import copy
+import xacc
+MASTER_DIRS = ['vqe']
+
+MASTER_PACKAGES = {}
+
+TOP_PATH = os.path.dirname(os.path.realpath(__file__))
+
+XACC_PYTHON_PLUGIN_PATH = "${CMAKE_INSTALL_PREFIX}/py-plugins/"
+
+PLUGIN_INSTALLATIONS = {}
+
+
+def parse_args(args):
+    parser = argparse.ArgumentParser(description="Installation manager for XACC benchmark plugins.",
+                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+                                     fromfile_prefix_chars='@')
+    parser.add_argument("-i", "--install", type=str, help="Install a plugin package to the XACC plugin directory.", required=False)
+    parser.add_argument("-l", "--list", help="List all available plugin packages for installation.", required=False, action='store_true')
+    parser.add_argument("-p", "--path", help="Set the XACC Python Plugin path (default = /root/.xacc/py-plugins)", required=False)
+
+    opts = parser.parse_args(args)
+    return opts
+
+def install_package(install_name):
+    try:
+        package_path = PLUGIN_INSTALLATIONS[install_name]
+    except KeyError as ex:
+        xacc.info("There is no '{}' XACC Python plugin package available.".format(install_name))
+        exit(1)
+
+    install_directive = os.path.join(package_path+"/install.ini") if os.path.isfile(package_path+"/install.ini") else None
+    plugin_files = []
+    if not install_directive:
+        plugin_files += [package_path+"/"+f for f in os.listdir(
+            package_path) if os.path.isfile(os.path.join(package_path, f)) and f.endswith(".py")]
+    else:
+        plugin_dict, plugin_list = read_install_directive(install_directive, package_path)
+        for k,v in plugin_dict.items():
+            mini_package_path = v
+            plugin_files += [v+"/"+f for f in os.listdir(v) if os.path.isfile(os.path.join(v, f)) and f.endswith(".py")]
+    n_plugins = len(plugin_files)
+    for plugin in plugin_files:
+        copy(os.path.join(plugin), XACC_PYTHON_PLUGIN_PATH)
+
+    xacc.info("Installed {} plugins from the '{}' package to the {} directory.".format(n_plugins, install_name, XACC_PYTHON_PLUGIN_PATH))
+
+def set_plugin_path(path):
+    global XACC_PYTHON_PLUGIN_PATH
+    XACC_PYTHON_PLUGIN_PATH = path
+
+def read_install_directive(install_file, parent):
+    config = configparser.RawConfigParser()
+    config.read(install_file)
+    results = {}
+    packages_here = []
+    for section in config.sections():
+        for installation in config.items(section):
+            name, folder = installation
+            packages_here.append(name)
+            folder = parent+folder
+            results.update(dict([(name, folder)]))
+    return results, packages_here
+
+def get_packages():
+
+    for mdir in MASTER_DIRS:
+        plugins_path = os.path.join(TOP_PATH, mdir)
+        install_directives = plugins_path+"/install.ini"
+        package_dict, package_list = read_install_directive(install_directives, plugins_path)
+        PLUGIN_INSTALLATIONS.update(package_dict)
+
+        if mdir not in MASTER_PACKAGES:
+            MASTER_PACKAGES[mdir] = package_list
+        else:
+            MASTER_PACKAGES[mdir] += package_list
+
+def main(argv=None):
+    opts = parse_args(sys.argv[1:])
+    get_packages()
+
+    if opts.path:
+        set_plugin_path(opts.path)
+
+    if opts.install:
+        install_package(opts.install)
+
+    if opts.list:
+        xacc.info("Available XACC Python plugin packages:")
+        for k, v in MASTER_PACKAGES.items():
+            xacc.info("{:5}: {!s}".format(k, v))
+
+if __name__ == "__main__":
+    sys.exit(main())
\ No newline at end of file

python/benchmark/vqe/chemistry/ansatz_generators/deuteron_h2_ansatz.py

+from pelix.ipopo.decorators import ComponentFactory, Property, Requires, Provides, \
+    Validate, Invalidate, Instantiate
+from ansatzgenerator import AnsatzGenerator
+import ast
+import xacc
+
+@ComponentFactory("dh2_ansatz_generator_factory")
+@Provides("ansatz_generator")
+@Property("_ansatz_generator", "ansatz_generator", "dh2")
+@Property("_name", "name", "dh2")
+@Instantiate("dh2_ansatz_generator_instance")
+class DH2(AnsatzGenerator):
+
+    def generate(self, inputParams, nQubits):
+        dh2_function = xacc.gate.createFunction('dh2', [x for x in range(nQubits)])
+        ry = xacc.gate.create("Ry", [1],['t1'])
+        cx = xacc.gate.create("CNOT", [1, 0])
+        dh2_function.addInstruction(ry)
+        dh2_function.addInstruction(cx)
+        return dh2_function

python/benchmark/vqe/chemistry/ansatz_generators/deuteron_h3_ansatz.py

+from pelix.ipopo.decorators import ComponentFactory, Property, Requires, Provides, \
+    Validate, Invalidate, Instantiate
+from ansatzgenerator import AnsatzGenerator
+import ast
+import xacc
+
+@ComponentFactory("dh3_ansatz_generator_factory")
+@Provides("ansatz_generator")
+@Property("_ansatz_generator", "ansatz_generator", "dh3")
+@Property("_name", "name", "dh3")
+@Instantiate("dh3_ansatz_generator_instance")
+class DeuteronH3(AnsatzGenerator):
+
+    def generate(self, inputParams, nQubits):
+        function = xacc.gate.createFunction('ob', [x for x in range(nQubits)])
+        x = xacc.gate.create('X', [0])
+        ry = xacc.gate.create("Ry", [1], ['t1'])
+        ry2 = xacc.gate.create('Ry', [2], ['t2'])
+        cx = xacc.gate.create("CNOT", [2, 0])
+        cx2 = xacc.gate.create("CNOT", [0, 1])
+        ry3 = xacc.gate.create('Ry', [1], ['-t1'])
+        cx3 = xacc.gate.create("CNOT", [0, 1])
+        cx4 = xacc.gate.create("CNOT", [1, 0])
+        function.addInstruction(x)
+        function.addInstruction(ry)
+        function.addInstruction(ry2)
+        function.addInstruction(cx)
+        function.addInstruction(cx2)
+        function.addInstruction(ry3)
+        function.addInstruction(cx3)
+        function.addInstruction(cx4)
+        return function

python/benchmark/vqe/chemistry/ansatz_generators/hartree_fock_ansatz.py

+from pelix.ipopo.decorators import ComponentFactory, Property, Requires, Provides, \
+    Validate, Invalidate, Instantiate
+from ansatzgenerator import AnsatzGenerator
+import ast
+import xacc
+
+@ComponentFactory("hf_ansatz_generator_factory")
+@Provides("ansatz_generator")
+@Property("_ansatz_generator", "ansatz_generator", "hf")
+@Property("_name", "name", "hf")
+@Instantiate("hf_ansatz_generator_instance")
+class HF(AnsatzGenerator):
+    """
+        AnsatzGenerator for calculating the Hartree-Fock state of molecules
+
+        Required input configurations | (examples):
+            x-gates | ([0,1])
+    """
+
+    def generate(self, inputParams, nQubits):
+        hf_function = xacc.gate.createFunction('hf', [x for x in range(nQubits)])
+        x_gate_qubits = ast.literal_eval(inputParams['x-gates'])
+        for q in x_gate_qubits:
+            x_gate = xacc.gate.create('X', [q])
+            hf_function.addInstruction(x_gate)
+        return hf_function

python/benchmark/vqe/chemistry/ansatz_generators/hwe.py

+from pelix.ipopo.decorators import ComponentFactory, Property, Requires, Provides, \
+    Validate, Invalidate, Instantiate
+from ansatzgenerator import AnsatzGenerator
+import ast
+import xacc
+
+@ComponentFactory("hwe_ansatz_generator_factory")
+@Provides("ansatz_generator")
+@Property("_ansatz_generator", "ansatz_generator", "hwe")
+@Property("_name", "name", "hwe")
+@Instantiate("hwe_ansatz_generator_instance")
+class HWE(AnsatzGenerator):
+    """
+        AnsatzGenerator for the hardware-efficient ansatz to calculate molecular energies
+
+        Required input configurations | (examples)
+            layers | (1)
+            connectivity | [[0,1],[1,2],[2,3]]
+    """
+
+    def generate(self, inputParams, nQubits):
+        ir_generator = xacc.getIRGenerator(inputParams['name'])
+        function = ir_generator.generate([int(inputParams['layers']), nQubits, inputParams['connectivity']])
+        return function

python/benchmark/vqe/chemistry/ansatz_generators/ucc1_ansatz.py

+from pelix.ipopo.decorators import ComponentFactory, Property, Requires, Provides, \
+    Validate, Invalidate, Instantiate
+from ansatzgenerator import AnsatzGenerator
+import ast
+import xacc
+
+@ComponentFactory("ucc1_ansatz_generator_factory")
+@Provides("ansatz_generator")
+@Property("_ansatz_generator", "ansatz_generator", "ucc1")
+@Property("_name", "name", "ucc1")
+@Instantiate("ucc1_ansatz_generator_instance")
+class UCC1(AnsatzGenerator):
+    """
+        AnsatzGenerator for the UCC1 ansatz for calculating molecular energies
+
+        Required input configurations:
+            x-gates | ([0, 2])
+    """
+
+    def generate(self, inputParams, nQubits):
+        ucc1_function = xacc.gate.createFunction('reduced-uccsd', [x for x in range(nQubits)])
+        x_gate_qubits = ast.literal_eval(inputParams['x-gates'])
+        for q in x_gate_qubits:
+            x_gate = xacc.gate.create('X', [q])
+            ucc1_function.addInstruction(x_gate)
+        rx1 = xacc.gate.create("Rx", [0], [1.57079632])
+        h1 = xacc.gate.create("H", [1])
+        h2 = xacc.gate.create("H", [2])
+        h3 = xacc.gate.create("H", [3])
+        cn1 = xacc.gate.create("CNOT", [0,1])
+        cn2 = xacc.gate.create("CNOT", [1,2])
+        cn3 = xacc.gate.create("CNOT", [2,3])
+        rz = xacc.gate.create("Rz", [3], ['t1'])
+        cn4 = xacc.gate.create("CNOT", [2,3])
+        cn5 = xacc.gate.create("CNOT", [1,2])
+        cn6 = xacc.gate.create("CNOT", [0,1])
+        rx2 = xacc.gate.create("Rx", [0], [-1.57079632])
+        h4 = xacc.gate.create("H", [1])
+        h5 = xacc.gate.create("H", [2])
+        h6 = xacc.gate.create("H", [3])
+        ucc1_function.addInstruction(rx1)
+        ucc1_function.addInstruction(h1)
+        ucc1_function.addInstruction(h2)
+        ucc1_function.addInstruction(h3)
+        ucc1_function.addInstruction(cn1)
+        ucc1_function.addInstruction(cn2)
+        ucc1_function.addInstruction(cn3)
+        ucc1_function.addInstruction(rz)
+        ucc1_function.addInstruction(cn4)
+        ucc1_function.addInstruction(cn5)
+        ucc1_function.addInstruction(cn6)
+        ucc1_function.addInstruction(rx2)
+        ucc1_function.addInstruction(h4)
+        ucc1_function.addInstruction(h5)
+        ucc1_function.addInstruction(h6)
+        return ucc1_function

python/benchmark/vqe/chemistry/ansatz_generators/ucc3_ansatz.py

+from pelix.ipopo.decorators import ComponentFactory, Property, Requires, Provides, \
+    Validate, Invalidate, Instantiate
+from ansatzgenerator import AnsatzGenerator
+import ast
+import xacc
+
+@ComponentFactory("ucc3_ansatz_generator_factory")
+@Provides("ansatz_generator")
+@Property("_ansatz_generator", "ansatz_generator", "ucc3")
+@Property("_name", "name", "ucc3")
+@Instantiate("ucc3_ansatz_generator_instance")
+class UCC3(AnsatzGenerator):
+
+    def generate(self, inputParams, nQubits):
+        ucc3_function = xacc.gate.createFunction('ucc3', [x for x in range(nQubits)])
+        x0 = xacc.gate.create("X", [0])
+        x1 = xacc.gate.create("X", [2])
+        rx0 = xacc.gate.create("Rx", [0], [1.57079632])
+        h0 = xacc.gate.create("H", [1])
+        cn0 = xacc.gate.create("CNOT", [0,1])
+        rz0 = xacc.gate.create("Rz", [1], ['t1'])
+        rx1 = xacc.gate.create("Rx", [2], [1.57079632])
+        h1 = xacc.gate.create("H", [3])
+        cn1 = xacc.gate.create("CNOT", [2,3])
+        rz1 = xacc.gate.create("Rz", [3], ['t2'])
+        cn2 = xacc.gate.create("CNOT", [2,3])
+        rx2 = xacc.gate.create("Rx", [2], [-1.57079632])
+        h2 = xacc.gate.create("H", [2])
+        cn3 = xacc.gate.create("CNOT", [1,2])
+        cn4 = xacc.gate.create("CNOT", [2,3])
+        rz2 = xacc.gate.create("Rz", [3], ['t3'])
+        cn5 = xacc.gate.create("CNOT", [2,3])
+        cn6 = xacc.gate.create("CNOT", [1,2])
+        cn7 = xacc.gate.create("CNOT", [0,1])
+        rx3 = xacc.gate.create("Rx", [0], [-1.57079632])
+        h4 = xacc.gate.create("H", [1])
+        h5 = xacc.gate.create("H", [2])
+        h6 = xacc.gate.create("H", [3])
+        ucc3_function.addInstruction(x0)
+        ucc3_function.addInstruction(x1)
+        ucc3_function.addInstruction(rx0)
+        ucc3_function.addInstruction(h0)
+        ucc3_function.addInstruction(cn0)
+        ucc3_function.addInstruction(rz0)
+        ucc3_function.addInstruction(rx1)
+        ucc3_function.addInstruction(h1)
+        ucc3_function.addInstruction(cn1)
+        ucc3_function.addInstruction(rz1)
+        ucc3_function.addInstruction(cn2)
+        ucc3_function.addInstruction(rx2)
+        ucc3_function.addInstruction(h2)
+        ucc3_function.addInstruction(cn3)
+        ucc3_function.addInstruction(cn4)
+        ucc3_function.addInstruction(rz2)
+        ucc3_function.addInstruction(cn5)
+        ucc3_function.addInstruction(cn6)
+        ucc3_function.addInstruction(cn7)
+        ucc3_function.addInstruction(rx3)
+        ucc3_function.addInstruction(h4)
+        ucc3_function.addInstruction(h5)
+        ucc3_function.addInstruction(h6)
+        return ucc3_function

python/benchmark/vqe/chemistry/ansatz_generators/uccsd.py

+from pelix.ipopo.decorators import ComponentFactory, Property, Requires, Provides, \
+    Validate, Invalidate, Instantiate
+from ansatzgenerator import AnsatzGenerator
+import ast
+import xacc
+
+@ComponentFactory("uccsd_ansatz_generator_factory")
+@Provides("ansatz_generator")
+@Property("_ansatz_generator", "ansatz_generator", "uccsd")
+@Property("_name", "name", "uccsd")
+@Instantiate("uccsd_ansatz_generator_instance")
+class UCCSD(AnsatzGenerator):
+    """
+        AnsatzGenerator for the UCCSD ansatz for calculating molecular energies
+
+        Required input configurations:
+            n-electrons | (2)
+    """
+
+    def generate(self, inputParams, nQubits):
+        ir_generator = xacc.getIRGenerator(inputParams['name'])
+        function = ir_generator.generate([int(inputParams['n-electrons']), nQubits])
+        return function

python/benchmark/vqe/chemistry/hamiltonian_generators/frozencore.py

+from pelix.ipopo.decorators import ComponentFactory, Property, Requires, Provides, \
+    Validate, Invalidate, Instantiate
+
+from hamiltonian_generator import HamiltonianGenerator
+import ast
+import xaccvqe
+import numpy as np
+
+try:
+    import psi4
+    from molutil import geometry
+    from openfermion.transforms import get_fermion_operator
+    from openfermion.transforms import get_fermion_operator, jordan_wigner
+    from openfermion.ops import InteractionOperator
+except:
+    pass
+
+@ComponentFactory("frozencore_hamiltonian_generator_factory")
+@Provides("hamiltonian_generator")
+@Property("_hamiltonian_generator", "hamiltonian_generator", "frozen-core")
+@Property("_name", "name", "frozen-core")
+@Instantiate("frozencore_hamiltonian_generator_instance")
+class FrozenCoreHamiltonian(HamiltonianGenerator):
+    """
+        HamiltonianGenerator for molecules using the frozen-core approximation of electronic structure
+
+        Required input configurations | (example):
+            basis | (sto-3g)
+            geometry | ('0 1\nNa  0.000000   0.0      0.0\nH   0.0        0.0  1.914388\nsymmetry c1')
+            frozen-spin-orbitals | [0, 1, 2, 3, 4, 10, 11, 12, 13, 14]
+            active-spin-orbitals | [5, 9, 15, 19]
+    """
+
+    def generate(self, inputParams):
+
+        g = str(inputParams['geometry'])
+        basis = inputParams['basis']
+        moleculeGeom = geometry(g)
+        psi4.set_options({'basis':'sto-3g',
+                          'scf_type':'pk',
+                          'mp2_type':'conv',
+                          'e_convergence': 1e-8,
+                          'd_convergence': 1e-8})
+        scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)
+        E_nucl = moleculeGeom.nuclear_repulsion_energy()
+
+        # Get MO coefficients from SCF wavefunction
+        # ==> ERIs <==
+        # Create instance of MintsHelper class:
+
+        mints = psi4.core.MintsHelper(scf_wfn.basisset())
+
+
+
+        nbf = mints.nbf()           # number of basis functions
+        nso = 2 * nbf               # number of spin orbitals
+        nalpha = scf_wfn.nalpha()   # number of alpha electrons
+        nbeta = scf_wfn.nbeta()     # number of beta electrons
+        nocc = nalpha + nbeta       # number of occupied orbitals
+        nvirt = 2 * nbf - nocc      # number of virtual orbitals
+        list_occ_alpha = np.asarray(scf_wfn.occupation_a())
+        list_occ_beta = np.asarray(scf_wfn.occupation_b())
+
+        # Get orbital energies
+        eps_a = np.asarray(scf_wfn.epsilon_a())
+        eps_b = np.asarray(scf_wfn.epsilon_b())
+        eps = np.append(eps_a, eps_b)
+
+        # Get orbital coefficients:
+        Ca = np.asarray(scf_wfn.Ca())
+        Cb = np.asarray(scf_wfn.Cb())
+        C = np.block([
+            [Ca,np.zeros_like(Cb)],
+            [np.zeros_like(Ca),Cb]
+        ])
+
+        # Get the two electron integrals using MintsHelper
+        Ints = np.asarray(mints.ao_eri())
+        def spin_block_tei(I):
+            """
+             Function that spin blocks two-electron integrals
+             Using np.kron, we project I into the space of the 2x2 identity, tranpose the result
+             and project into the space of the 2x2 identity again. This doubles the size of each axis.
+             The result is our two electron integral tensor in the spin orbital form.
+            """
+            identity = np.eye(2)
+            I = np.kron(identity, I)
+            return np.kron(identity, I.T)
+
+        # Spin-block the two electron integral array
+        I_spinblock = spin_block_tei(Ints)
+
+        # Converts chemist's notation to physicist's notation, and antisymmetrize
+        # (pq | rs) ---> <pr | qs>
+        # Physicist's notation
+        tmp = I_spinblock.transpose(0, 2, 1, 3)
+
+        # Antisymmetrize:
+        # <pr||qs> = <pr | qs> - <pr | sq>
+        gao = tmp - tmp.transpose(0, 1, 3, 2)
+        gmo = np.einsum('pQRS, pP -> PQRS',
+                        np.einsum('pqRS, qQ -> pQRS',
+                                  np.einsum('pqrS, rR -> pqRS',
+                                            np.einsum('pqrs, sS -> pqrS', gao, C), C), C), C)
+
+        # -------- 0-body term:
+        Hamiltonian_0body = E_nucl
+
+        # -------- 1-body term:
+        #   Ca*
+        # Build core Hamiltonian
+        T = np.asarray(mints.ao_kinetic())
+        V = np.asarray(mints.ao_potential())
+        H_core_ao = T + V
+
+
+
+        # -- check  which one more efficient (matmul vs einsum)
+        #   H_core_mo = np.matmul(Ca.T,np.matmul(H_core_ao,Ca)))
+        H_core_mo = np.einsum('ij, jk, kl -> il', Ca.T, H_core_ao, Ca)
+        H_core_mo_alpha = H_core_mo
+        H_core_mo_beta = H_core_mo
+
+        # ---- this version breaks is we permuted  SCF eigvecs
+        # Hamiltonian_1body = np.block([
+        #            [  H_core_mo_alpha             ,           np.zeros_like(H_core_mo_alpha)],
+        #            [np.zeros_like(H_core_mo_beta) ,  H_core_mo_beta      ]])
+        #
+        # --- th is version is  safer than above (H_1b is permutted correctly if eigvecs are permutted)
+        Hamiltonian_1body_ao = np.block([[H_core_ao, np.zeros_like(H_core_ao)],
+                                         [np.zeros_like(H_core_ao), H_core_ao]])
+        Hamiltonian_1body = np.einsum('ij, jk, kl -> il', C.T, Hamiltonian_1body_ao, C)
+        Hamiltonian_2body = gmo
+        MSO_frozen_list = ast.literal_eval(inputParams['frozen-spin-orbitals'])
+        MSO_active_list = ast.literal_eval(inputParams['active-spin-orbitals'])
+        n_frozen = len(MSO_frozen_list)
+        n_active = len(MSO_active_list)
+
+        # ----- 0-body frozen-core:
+        Hamiltonian_fc_0body = E_nucl
+        for a in range(n_frozen):
+
+            ia = MSO_frozen_list[a]
+            Hamiltonian_fc_0body += Hamiltonian_1body[ia, ia]
+
+            for b in range(a):
+
+                ib = MSO_frozen_list[b]
+                Hamiltonian_fc_0body += gmo[ia, ib, ia, ib]
+
+
+
+        # --- 1-body frozen-core:
+        Hamiltonian_fc_1body = np.zeros((n_active, n_active))
+        Hamiltonian_fc_1body_tmp = np.zeros((n_active, n_active))
+        for p in range(n_active):
+
+            ip = MSO_active_list[p]
+
+            for q in range(n_active):
+
+                iq = MSO_active_list[q]
+                Hamiltonian_fc_1body[p, q] = Hamiltonian_1body[ip, iq]
+                #Hamiltonian_fc_1body_tmp[p,q] =  Hamiltonian_1body[ip,iq]
+
+                for a in range(n_frozen):
+
+                    ia = MSO_frozen_list[a]
+                    Hamiltonian_fc_1body[p, q] += gmo[ia, ip, ia, iq]
+
+
+
+        # ------- 2-body frozen-core:
+
+        Hamiltonian_fc_2body = np.zeros((n_active, n_active, n_active, n_active))
+        for p in range(n_active):
+