Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Final run #142

Merged
merged 8 commits into from
Aug 27, 2024
Merged

Final run #142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
dc78fb4
add example of toml usage
RichRick1 Aug 16, 2024
80f8a17
fix toml for heisenberg model
RichRick1 Aug 16, 2024
117e47a
update examples
RichRick1 Aug 26, 2024
82c266e
fix pycodestyle
RichRick1 Aug 26, 2024
2d23a96
fix test
RichRick1 Aug 26, 2024
96dc771
fix PariserParr module
RichRick1 Aug 27, 2024
af82620
fix tutorial
RichRick1 Aug 27, 2024
d657f3a
fix pycodestyle
RichRick1 Aug 27, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
258 changes: 72 additions & 186 deletions examples/Real_examples1.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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). "
]
},
{
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand Down Expand Up @@ -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"
]
}
],
Expand All @@ -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",
Expand Down Expand Up @@ -415,7 +301,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.9.7"
}
},
"nbformat": 4,
Expand Down
Loading
Loading