Skip to content

Commit

Permalink
Merge pull request #142 from RichRick1/PariserParr
Browse files Browse the repository at this point in the history
  • Loading branch information
@giovanni-br
giovanni-br authored Aug 27, 2024
2 parents 9b750a8 + d657f3a commit 14d46e9
Show file tree
Hide file tree
Showing 6 changed files with 121 additions and 229 deletions.
258 changes: 72 additions & 186 deletions examples/Real_examples1.ipynb
View file
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

0 comments on commit 14d46e9

Please sign in to comment.