diff --git a/examples/Real_examples1.ipynb b/examples/Real_examples1.ipynb index a72f64b..39f2492 100644 --- a/examples/Real_examples1.ipynb +++ b/examples/Real_examples1.ipynb @@ -4,10 +4,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Tutorial: Benzene and Biphenyl\n", + "# Tutorial: Benzene\n", "\n", "## Introduction\n", - "In this tutorial, we'll focus on generating integrals for the PPP model considering real organic molecules: Benzene and Biphenyl " + "In this tutorial, we'll focus on generating integrals for the PPP model considering real organic molecule -- Benzene.\n", + "This example is used to showcase the simplicity of usage a ModelHamiltoian package to build a Hamiltonian that corresponds to the one found in the [literature](https://www.sciencedirect.com/science/article/pii/S0301010499000178). " ] }, { @@ -25,43 +26,52 @@ "source": [ "## Example: Defining PPP Hamiltonian\n", "\n", - "In order to define a PPP hamiltonina, it's necessary to provide the parameter `gamma`. Performing a Full Configuration Iteraction of the hamiltonian, varying $\\beta$\n" + "There are several ways to build a PPP Hamiltonian. The most straightforward way is to evaluate $gamma$ matrix and provide it to the `ModelHamiltonian` class. \n" ] }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "65.30151\n", "Beta: 0.0, FCI Energy: -1.4210854715202004e-14\n", - "65.30151\n", - "Beta: -1.0, FCI Energy: -2.779514577625193\n", - "65.30151\n", + "Beta: -1.0, FCI Energy: -2.779514577625207\n", "Beta: -2.0, FCI Energy: -9.084603708594514\n", - "65.30151\n", "Beta: -3.0, FCI Energy: -16.484889549995387\n", - "65.30151\n", - "Beta: -4.0, FCI Energy: -24.19501406479904\n", - "65.30151\n", - "Beta: -5.0, FCI Energy: -32.024692914537425\n" + "Beta: -4.0, FCI Energy: -24.195014064799068\n", + "Beta: -5.0, FCI Energy: -32.024692914537454\n" ] } ], "source": [ - "import sys \n", - "sys.path.insert(0, '../')\n", + "# import sys \n", + "# sys.path.insert(0, '../')\n", "import numpy as np \n", "import sys \n", "import numpy as np \n", "from moha import HamPPP \n", "from pyscf import gto, scf, fci \n", "\n", - "RXY = np.array([0.00504725, 0.09930802, 0.16008021])\n", + "def build_gamma(norb=6):\n", + " # Define the gamma matrix here (example shape) \n", + " gamma_matrix = np.zeros((norb, norb))\n", + "\n", + " gamma_a_a1 = 5.27760 \n", + " gamma_a_a2 = 3.86206 \n", + " gamma_a_a3 = 3.48785 \n", + " \n", + " # Fill the matrix according to the given rules \n", + " for a in range(norb): \n", + " gamma_matrix[a, a-1] = gamma_a_a1 \n", + " gamma_matrix[a, a-2] = gamma_a_a2 \n", + " gamma_matrix[a, a-3] = gamma_a_a3 \n", + " \n", + " gamma_matrix = np.maximum(gamma_matrix, np.transpose(gamma_matrix))\n", + " return gamma_matrix\n", " \n", "# Define the benzene system with a 6-membered ring structure \n", "system = [('C1', 'C2', 1), ('C2', 'C3', 1), \n", @@ -70,40 +80,19 @@ "\n", "\n", "#Define U_onsite, alpha and gamma matrix \n", - "u_onsite = np.array([10.84, 10.84, 10.84, 10.84, 10.84, 10.84])/2\n", + "u_onsite = 10.84 * np.ones(6) / 2\n", "alpha = 0 \n", " \n", - "# Define the gamma matrix here (example shape) \n", - "norb = 6 # Number of orbitals for benzene (6 carbons) \n", - "gamma_matrix = np.zeros((norb, norb))\n", - "\n", - "gamma_a_a1 = 5.27760 \n", - "gamma_a_a2 = 3.86206 \n", - "gamma_a_a3 = 3.48785 \n", - " \n", - "# Fill the matrix according to the given rules \n", - "for a in range(norb): \n", - " gamma_matrix[a-1, a] = gamma_a_a1 \n", - " gamma_matrix[a-2, a] = gamma_a_a2 \n", - " gamma_matrix[a-3, a] = gamma_a_a3 \n", - " \n", - " gamma_matrix[a, a-1] = gamma_a_a1 \n", - " gamma_matrix[a, a-2] = gamma_a_a2 \n", - " gamma_matrix[a, a-3] = gamma_a_a3\n", - "\n", - " gamma_matrix[a, a] = 10.84\n", - " \n", - "gamma_lib = gamma_matrix - np.diag(np.diag(gamma_matrix))\n", + "gamma_matrix = build_gamma()\n", "# Loop to vary beta from 0 to -5\n", "for beta in np.linspace(0, -5, num=6):\n", " # Generate the PPP object\n", - " ppp_hamiltonian = HamPPP(system, alpha=alpha, beta=beta, u_onsite=u_onsite, gamma=gamma_lib, charges=np.array([1 for _ in range(6)]))\n", + " ppp_hamiltonian = HamPPP(system, alpha=alpha, beta=beta, u_onsite=u_onsite, gamma=gamma_matrix, charges=np.ones(6))\n", " \n", " # Generate the 0, 1, 2 body integral elements\n", " e01 = ppp_hamiltonian.generate_zero_body_integral() \n", " h11 = ppp_hamiltonian.generate_one_body_integral(dense=True, basis='spatial basis') \n", - " h2 = ppp_hamiltonian.generate_two_body_integral(dense=True, basis='spatial basis',sym = 4) \n", - " print(e01)\n", + " h2 = ppp_hamiltonian.generate_two_body_integral(dense=True, basis='spatial basis', sym=4) \n", " \n", " # Convert h2 to chemists notation \n", " h21_ch = 1*np.transpose(h2, (0, 2, 1, 3)) \n", @@ -121,23 +110,29 @@ "source": [ "# PPP Benzene: Building gamma\n", "\n", + "What if we want to build $gamma$ matrix for a molecule? The power and flexibility of the `ModelHamiltonian` package allows us to do that.\n", "\n", - "MoHa package offers a functionality to compute `gamma` using the module `PariserParr` and the function `compute_gamma`. \n", - "\n", - "In order to do that, user can provide all the paramters of the function or let the library compute with the default values:\n", - "\n", - "* If user does not provide `U_xy`(potential energies for sites) or `R_xy`matrix with distances of the sites; Those values are going to be computed using the distances provided by `connectivity`, using the formulas below, based on the deafult values of `affinity_dct` and `atom_dict`:\n", "\n", - "$$ U_x = \\alpha_x - A_x $$\n", + "MoHa package offers a functionality to compute $gamma$ using the module `PariserParr` and the function `compute_gamma`. \n", + "Computation of $gamma$ matrix requires two ingridients: \n", + "* `u_onsite` - potential on each site/atoms\n", + "* `R_xy` - distances between sites/atoms\n", "\n", + "Note: you don't need to provide any of those parameters. If `u_onsite` is not provided, the matrix $\\overline{U}_{XY}$ will be populated using the formula:\n", "$$\n", "\\overline{U}_{XY} = \\frac{1}{2}(U_X + U_Y)\n", - "$$\n", + "$$ \n", + "\n", + "where \n", + "$$ U_x = I_x - A_x $$\n", + "Values of `I_x` and `A_x` ionisation potential and electron affinity of the atom `x` respectively. User can provide those values in the form of a dictionary `affinity_dct` and `atom_dict` or use the default values.\n", "\n", - "User can also provide this dictionaries(affinity_dct` and `atom_dict`) to the function and proceed in the same way.\n", "\n", + "If `R_xy` is not provided, it will be automatically populated from the distances provided at `connectivity` parameter.\n", "\n", - "Finally, `gamma` is computed using the formula:\n", + "\n", + "\n", + "Finally, $gamma$ is computed using the formula:\n", "\n", "$$\n", "\\gamma_{XY} = \\frac{\\overline{U}_{XY}}{\\overline{U}_{XY} R_{XY} + e^{-\\frac{1}{2} \\overline{U}_{XY}^2 R_{XY}^2}}\n", @@ -153,167 +148,59 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Defining the Benzene system following th same sttructure defined in the previous tutorials, we decide to move on with Full Configuration Iteraction of the hamiltonian, defining `gamma`\n", + "Defining the Benzene system following th same structure defined in the previous tutorials, we decide to move on with Full Configuration Iteraction of the hamiltonian, building $gamma$ matrix using the `PariserParr` module.\n", "\n", "Varying $\\beta$ from 0 to -5, in this piece of code we prove that we can achieve the spectrum of energies reported in Bendazzoli et al." ] }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "65.30151\n", - "Beta: 0.0, FCI Energy: -1.4210854715202004e-14\n", - "65.30151\n", - "Beta: -1.0, FCI Energy: -2.779514577625193\n", - "65.30151\n", - "Beta: -2.0, FCI Energy: -9.084603708594514\n", - "65.30151\n", - "Beta: -3.0, FCI Energy: -16.484889549995387\n", - "65.30151\n", - "Beta: -4.0, FCI Energy: -24.19501406479904\n", - "65.30151\n", - "Beta: -5.0, FCI Energy: -32.024692914537425\n" + "Beta: 0.0, FCI Energy: 8.881784197001252e-16\n", + "Beta: -1.0, FCI Energy: -1.6000550795826625\n", + "Beta: -2.0, FCI Energy: -5.956135566002675\n", + "Beta: -3.0, FCI Energy: -12.105432809022926\n", + "Beta: -4.0, FCI Energy: -19.126580752154776\n", + "Beta: -5.0, FCI Energy: -26.54429608422473\n" ] } ], "source": [ - "import sys \n", - "sys.path.insert(0, '../')\n", - "import numpy as np \n", - "import sys \n", - "import numpy as np \n", - "from moha import HamPPP \n", - "from pyscf import gto, scf, fci \n", - "\n", - "RXY = np.array([0.00504725, 0.09930802, 0.16008021])\n", - " \n", - "# Define the benzene system with a 6-membered ring structure \n", - "system = [('C1', 'C2', 1), ('C2', 'C3', 1), \n", - " ('C3', 'C4', 1), ('C4', 'C5', 1), \n", - " ('C5', 'C6', 1), ('C6', 'C1', 1)] \n", - "\n", - "\n", - "#Define U_onsite, alpha and gamma matrix \n", - "u_onsite = np.array([10.84, 10.84, 10.84, 10.84, 10.84, 10.84])/2\n", - "alpha = 0 \n", - " \n", - "# Define the gamma matrix here (example shape) \n", - "norb = 6 # Number of orbitals for benzene (6 carbons) \n", - "gamma_matrix = np.zeros((norb, norb))\n", - "\n", - "gamma_a_a1 = 5.27760 \n", - "gamma_a_a2 = 3.86206 \n", - "gamma_a_a3 = 3.48785 \n", - " \n", - "# Fill the matrix according to the given rules \n", - "for a in range(norb): \n", - " gamma_matrix[a-1, a] = gamma_a_a1 \n", - " gamma_matrix[a-2, a] = gamma_a_a2 \n", - " gamma_matrix[a-3, a] = gamma_a_a3 \n", - " \n", - " gamma_matrix[a, a-1] = gamma_a_a1 \n", - " gamma_matrix[a, a-2] = gamma_a_a2 \n", - " gamma_matrix[a, a-3] = gamma_a_a3\n", - "\n", - " gamma_matrix[a, a] = 10.84\n", - " \n", - "gamma_lib = gamma_matrix - np.diag(np.diag(gamma_matrix))\n", - "# Loop to vary beta from 0 to -5\n", - "for beta in np.linspace(0, -5, num=6):\n", - " # Generate the PPP object\n", - " ppp_hamiltonian = HamPPP(system, alpha=alpha, beta=beta, u_onsite=u_onsite, gamma=gamma_lib, charges=np.array([1 for _ in range(6)]))\n", - " \n", - " # Generate the 0, 1, 2 body integral elements\n", - " e01 = ppp_hamiltonian.generate_zero_body_integral() \n", - " h11 = ppp_hamiltonian.generate_one_body_integral(dense=True, basis='spatial basis') \n", - " h2 = ppp_hamiltonian.generate_two_body_integral(dense=True, basis='spatial basis',sym = 4) \n", - " print(e01)\n", - " \n", - " # Convert h2 to chemists notation \n", - " h21_ch = 1*np.transpose(h2, (0, 2, 1, 3)) \n", - " norb = h11.shape[0] \n", - " nelec = norb # Assuming one electron per site \n", - "\n", - " # Perform the FCI calculation \n", - " e_fci, ci = fci.direct_spin1.kernel(1*h11, 2*h21_ch, norb, nelec, ecore=e01, nroots=1)\n", - " print(f'Beta: {beta}, FCI Energy: {e_fci}')\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Pariser Parr Atomic dictionary\n", - "\n", - "The previous approach requires to build the $\\gamma$ matrix directly, but MoHa package supports a method of computing the matrix based on the distance($R_{XY}$) beetween the atoms and the potential($U_{XY}$). \n", - "\n", - "\n", - "If the value of `gamma` is not provided, it is computed using the formula:\n", - "\n", - "$$\n", - "\\gamma_{XY} = \\frac{\\overline{U}_{XY}}{\\overline{U}_{XY} R_{XY} + e^{-\\frac{1}{2} \\overline{U}_{XY}^2 R_{XY}^2}}\n", - "$$\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 31, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Beta: 0.0, FCI Energy: 7.105427357601002e-15\n", - "Beta: -1.0, FCI Energy: -2.950976576402695\n", - "Beta: -2.0, FCI Energy: -9.425258557059209\n", - "Beta: -3.0, FCI Energy: -16.75428207213497\n", - "Beta: -4.0, FCI Energy: -24.36705165013518\n", - "Beta: -5.0, FCI Energy: -32.1115278385854\n" - ] - } - ], - "source": [ - "import sys\n", - "import numpy as np\n", - "from moha import HamPPP\n", "from moha.rauk import PariserParr\n", - "from pyscf import gto, scf, fci\n", "\n", - "# Define the benzene system with a 6-membered ring structure\n", + "def build_Rxy_matrix(norb=6):\n", + " d = np.zeros((norb, norb))\n", + " for mu in range(6):\n", + " for nu in range(6):\n", + " d[mu, nu] = 1.4 * np.sin(np.pi * abs(mu - nu) / norb) / np.sin(np.pi / norb)\n", + " return d\n", "\n", + "# Define the benzene system with a 6-membered ring structure\n", "system = [('C1', 'C2', 1), ('C2', 'C3', 1), \n", " ('C3', 'C4', 1), ('C4', 'C5', 1), \n", " ('C5', 'C6', 1), ('C6', 'C1', 1)] \n", "\n", - "# Define the distances\n", - "RXY = [0.00504725, 0.09930802, 0.16008021]\n", + "Rxy_matrix = build_Rxy_matrix()\n", "\n", - "# Build the adjacency matrix directly\n", - "Rxy_matrix = np.array([\n", - " [0.0, RXY[0], RXY[1], RXY[2], RXY[1], RXY[0]],\n", - " [RXY[0], 0.0, RXY[0], RXY[1], RXY[2], RXY[1]],\n", - " [RXY[1], RXY[0], 0.0, RXY[0], RXY[1], RXY[2]],\n", - " [RXY[2], RXY[1], RXY[0], 0.0, RXY[0], RXY[1]],\n", - " [RXY[1], RXY[2], RXY[1], RXY[0], 0.0, RXY[0]],\n", - " [RXY[0], RXY[1], RXY[2], RXY[1], RXY[0], 0.0]\n", - "])\n", + "# compute u_onsite for the given system\n", + "u, _ = PariserParr.compute_u(system, ionization_dictionary={'C': 11.16}, # ionization energy of carbon in eV\n", + " affinity_dictionary={'C': 0.03}) # electron affinity of carbon in eV\n", "\n", - "u_onsite = np.array([10.84, 10.84, 10.84, 10.84, 10.84, 10.84]) / 2\n", - "gamma_matrix = PariserParr.compute_gamma(system, U_xy=u_onsite, Rxy_matrix=Rxy_matrix)\n", + "# compute gamma for the given system\n", + "gamma_matrix_ppp = PariserParr.compute_gamma(system, Rxy_matrix=Rxy_matrix,\n", + " u_onsite=u/2)\n", "\n", "alpha = 0\n", "\n", "# Loop to vary beta from 0 to -5\n", "for beta in np.linspace(0, -5, num=6):\n", - " ppp_hamiltonian = HamPPP(system, alpha=alpha, gamma = gamma_matrix, beta=beta, u_onsite=u_onsite, charges=np.array([1 for _ in range(6)]))\n", + " ppp_hamiltonian = HamPPP(system, alpha=alpha, gamma = gamma_matrix_ppp, beta=beta, u_onsite=u/2, charges=np.ones(6))\n", " \n", " # Generate the 0, 1, 2 body integral elements\n", " e0 = ppp_hamiltonian.generate_zero_body_integral() \n", @@ -343,14 +230,14 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Alpha: -0.414, Beta: -0.0533, FCI Energy: -5.468422237553814\n" + "Alpha: -0.414, Beta: -0.0533, FCI Energy: -29.85842553617428\n" ] } ], @@ -366,9 +253,8 @@ " ('C5', 'C6', 1), ('C6', 'B1', 1)] \n", " \n", "u_onsite = np.array([10.84, 10.84, 10.84, 10.84, 10.84, 10.84]) / 2\n", + "gamma_matrix = PariserParr.compute_gamma(system, u_onsite=u_onsite)\n", "\n", - "gamma_matrix = PariserParr.compute_gamma(system, U_xy=u_onsite)\n", - " \n", "norb = 6 # Number of orbitals for benzene (6 carbons) \n", "ppp_hamiltonian = HamPPP(system, u_onsite=u_onsite, gamma=gamma_matrix, charges=np.array([1 for _ in range(6)]))\n", "\n", @@ -415,7 +301,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.9.7" } }, "nbformat": 4, diff --git a/examples/toml/huckel.toml b/examples/toml/huckel.toml index 7d07fc2..a2c9934 100644 --- a/examples/toml/huckel.toml +++ b/examples/toml/huckel.toml @@ -13,3 +13,6 @@ nelec = 6 [model] hamiltonian = "Huckel" basis = "spatial basis" +alpha = 1.0 +beta = -2.0 +u_onsite = 1.0 \ No newline at end of file diff --git a/examples/toml/make_ham.py b/examples/toml/make_ham.py new file mode 100644 index 0000000..b1e4e3e --- /dev/null +++ b/examples/toml/make_ham.py @@ -0,0 +1,9 @@ +import moha +from moha.toml_tools import from_toml + +ham = from_toml("hubbard.toml") +zero_body = ham.generate_zero_body_integral() +one_body = ham.generate_one_body_integral(basis='spatial basis', dense=True) +two_body = ham.generate_two_body_integral(basis='spatial basis', dense=True) + +print(zero_body, one_body.shape, two_body.shape) \ No newline at end of file diff --git a/moha/hamiltonians.py b/moha/hamiltonians.py index 0e73f40..e51317e 100644 --- a/moha/hamiltonians.py +++ b/moha/hamiltonians.py @@ -421,7 +421,6 @@ def __init__( alpha=alpha, beta=beta, u_onsite=None, - gamma=None, sym=sym, atom_dictionary=atom_dictionary, bond_dictionary=bond_dictionary, diff --git a/moha/rauk/PariserParr.py b/moha/rauk/PariserParr.py index a01b6b7..c376d35 100644 --- a/moha/rauk/PariserParr.py +++ b/moha/rauk/PariserParr.py @@ -274,7 +274,7 @@ def compute_overlap( return one_body -def compute_gamma(connectivity, U_xy=None, Rxy_matrix=None, +def compute_gamma(connectivity, u_onsite=None, Rxy_matrix=None, atom_dictionary=None, affinity_dct=None, adjacency=None): r""" Calculate the gamma values for each pair of sites. @@ -293,35 +293,30 @@ def compute_gamma(connectivity, U_xy=None, Rxy_matrix=None, np.ndarray: A matrix of computed gamma values for each pair. Non-connected pairs have a gamma value of zero. """ - if Rxy_matrix is None and U_xy is None: - U_xy, Rxy_matrix = compute_u( + if Rxy_matrix is None and u_onsite is None: + u_onsite, Rxy_matrix = compute_u( connectivity, atom_dictionary, affinity_dct) elif Rxy_matrix is None: _, Rxy_matrix = compute_u( connectivity, atom_dictionary, affinity_dct) - num_sites = len(U_xy) + elif u_onsite is None: + u_onsite, _ = compute_u( + connectivity, atom_dictionary, affinity_dct) + num_sites = len(u_onsite) gamma_matrix = np.zeros((num_sites, num_sites)) - for tpl in connectivity: - atom1, atom2 = tpl[0], tpl[1] - atom1_name, site1 = get_atom_type(atom1) - atom2_name, site2 = get_atom_type(atom2) - # Get_atom_type returns site index starting from 1 - # So we need to subtract 1 to get the correct index - site1, site2 = site1 - 1, site2 - 1 - Ux = U_xy[site1] - Uy = U_xy[site2] - Rxy = Rxy_matrix[site1][site2] - Uxy_bar = 0.5 * (Ux + Uy) - gamma = Uxy_bar / \ - (Uxy_bar * Rxy + np.exp(-0.5 * Uxy_bar**2 * Rxy**2)) - gamma_matrix[site1][site2] = gamma - gamma_matrix[site2][site1] = gamma # Ensure symmetry + Uxy_bar = np.zeros((num_sites, num_sites)) + for i, ux in enumerate(u_onsite): + for j, uy in enumerate(u_onsite): + if i != j: + Uxy_bar[i, j] = 0.5 * (ux + uy) + gamma_matrix = Uxy_bar / \ + (Uxy_bar * Rxy_matrix + np.exp(-0.5 * Uxy_bar**2 * Rxy_matrix**2)) return gamma_matrix -def compute_u(connectivity, atom_dictionary, affinity_dictionary): +def compute_u(connectivity, ionization_dictionary, affinity_dictionary): r""" Calculate the onsite potential energy (U) for each site. @@ -329,10 +324,10 @@ def compute_u(connectivity, atom_dictionary, affinity_dictionary): ---------- connectivity (list of tuples): Each tuple contains indices of two sites (atom1, atom2) and the distance between them. - atom_dictionary (dict): Dictionary mapping atom types to their - ionization energies. + ionization_dictionary (dict): Dictionary mapping atom types to their + ionization energies. affinity_dictionary (dict): Dictionary mapping atom types to their - electron affinities. + electron affinities. Returns ---------- @@ -343,15 +338,14 @@ def compute_u(connectivity, atom_dictionary, affinity_dictionary): 2. List of float: Distances corresponding to each connectivity tuple. """ - if atom_dictionary is None: - atom_dictionary = {} - hx_dictionary_path = Path(__file__).parent / "hx_dictionary.json" - hx_dictionary = json.load(open(hx_dictionary_path, "rb")) - alpha_c = -0.414 # Value for sp2 orbital of Carbon atom. - beta_c = -0.0533 # Value for sp2 orbitals of Carbon atom. - for key, value in hx_dictionary.items(): - hx_value = value * abs(beta_c) - atom_dictionary[key] = alpha_c + hx_value + if ionization_dictionary is None: + ionization_dictionary = {} + ionization_dictionary_path = Path(__file__).parent / "ionization.json" + ionization_dictionary = json.load( + open(ionization_dictionary_path, "rb") + ) + for key, value in ionization_dictionary.items(): + ionization_dictionary[key] = value if affinity_dictionary is None: affinity_path = Path(__file__).parent / "affinity.json" @@ -372,7 +366,7 @@ def compute_u(connectivity, atom_dictionary, affinity_dictionary): # So we need to subtract 1 to get the correct index site1, site2 = site1 - 1, site2 - 1 - ionization = atom_dictionary[atom1_name] + ionization = ionization_dictionary[atom1_name] affinity = affinity_dictionary[atom1_name] U_x = ionization - affinity @@ -385,10 +379,10 @@ def compute_u(connectivity, atom_dictionary, affinity_dictionary): if len(connectivity) != max_site: tpl = connectivity[-1] atom_name, _ = get_atom_type(tpl[1]) - ionization = atom_dictionary[atom_name] + ionization = ionization_dictionary[atom_name] affinity = affinity_dictionary[atom_name] U_x = ionization - affinity u_onsite.append(U_x) - return u_onsite, Rxy_matrix + return np.array(u_onsite), Rxy_matrix diff --git a/moha/toml_tools/tools.py b/moha/toml_tools/tools.py index d42a14e..3054fd7 100644 --- a/moha/toml_tools/tools.py +++ b/moha/toml_tools/tools.py @@ -308,33 +308,34 @@ def build_moha(data): if data["model"]["hamiltonian"] == "ppp": charge_arr = charge * np.ones(norb) u_onsite_arr = u_onsite * np.ones(norb) - ham = moha.HamPPP(connectivity=adjacency, alpha=alpha, beta=beta, - u_onsite=u_onsite_arr, charges=charge_arr) + ham = moha.HamPPP(adjacency=adjacency, alpha=alpha, beta=beta, + u_onsite=u_onsite_arr, charges=charge_arr, + gamma=np.zeros((norb, norb))) return ham # Huckel elif data["model"]["hamiltonian"] == "huckel": - ham = moha.HamHuck(connectivity=adjacency, alpha=alpha, beta=beta) + ham = moha.HamHuck(adjacency=adjacency, alpha=alpha, beta=beta) return ham # Hubbard elif data["model"]["hamiltonian"] == "hubbard": u_onsite_arr = u_onsite * np.ones(norb) - ham = moha.HamHub(connectivity=adjacency, + ham = moha.HamHub(adjacency=adjacency, alpha=alpha, beta=beta, u_onsite=u_onsite_arr) return ham # -- Spin models --# # Heisenberg elif data["model"]["hamiltonian"] == "heisenberg": - ham = moha.HamHeisenberg(connectivity=adjacency, + ham = moha.HamHeisenberg(adjacency=adjacency, mu=mu, J_eq=J_eq, J_ax=J_ax) return ham # Ising elif data["model"]["hamiltonian"] == "ising": - ham = moha.HamIsing(connectivity=adjacency, mu=mu, J_ax=J_ax) + ham = moha.HamIsing(adjacency=adjacency, mu=mu, J_ax=J_ax) return ham # Richardson-Gaudin elif data["model"]["hamiltonian"] == "rg": - ham = moha.HamRG(connectivity=adjacency, mu=mu, J_eq=J_eq) + ham = moha.HamRG(adjacency=adjacency, mu=mu, J_eq=J_eq) return ham else: raise ValueError("Model hamiltonian " + data["model"]["hamiltonian"] +