Skip to content

Commit

Permalink
Merge pull request #119 from giovanni-br/main
Browse files Browse the repository at this point in the history 
Issue #116 Fixed
  • Loading branch information
@giovanni-br
giovanni-br authored Jun 9, 2024
2 parents f2090ad + a864ff5 commit 6c74a58
Show file tree
Hide file tree
Showing 4 changed files with 404 additions and 19 deletions.
160 changes: 155 additions & 5 deletions examples/rauk.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,10 +213,10 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[-0.41380839 -0.77906404 0. -1.87261684]\n",
" [-0.77906404 -0.47655051 -1.95444655 0. ]\n",
" [ 0. -1.95444655 -0.64027609 -1.64472969]\n",
" [-1.87261684 0. -1.64472969 -0.29956945]]\n"
"[[-0.41380839 -0.77906404 0. -0.62420561]\n",
" [-0.77906404 -0.47655051 -0.97722328 0. ]\n",
" [ 0. -0.97722328 -0.64027609 -0.82236485]\n",
" [-0.62420561 0. -0.82236485 -0.29956945]]\n"
]
}
],
Expand All @@ -237,11 +237,161 @@
"system = [('C1', 'Cl2', 1.1), ('Cl2', 'F3', 1.0),\n",
" ('F3', 'Si4', 1.0), ('Si4', 'C1', 1.0)]\n",
"hubbard = HamHub(system, u_onsite=np.array([1, 1]),\n",
" orbital_overlap={'C,Cl': 1, 'Cl,F': 2, 'F,Si': 2, 'Si,C': 3})\n",
" orbital_overlap=np.array([1,1,1,1]))\n",
"\n",
"print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: Second Body Terms: PariserParr Module"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Zero energy: 0\n",
"One body integrals in spatial basis: \n",
" [[-0.41380839 -0.77906404 0. -0.62420561]\n",
" [-0.77906404 -0.47655051 -0.97722328 0. ]\n",
" [ 0. -0.97722328 -0.64027609 -0.82236485]\n",
" [-0.62420561 0. -0.82236485 -0.29956945]]\n",
"Shape of two body integral in spatial basis: (4, 4, 4, 4)\n",
"------------------------------------------------------------\n",
"One body integrals in spin basis: \n",
" [[-0.41380839 -0.77906404 0. -0.62420561 0. 0.\n",
" 0. 0. ]\n",
" [-0.77906404 -0.47655051 -0.97722328 0. 0. 0.\n",
" 0. 0. ]\n",
" [ 0. -0.97722328 -0.64027609 -0.82236485 0. 0.\n",
" 0. 0. ]\n",
" [-0.62420561 0. -0.82236485 -0.29956945 0. 0.\n",
" 0. 0. ]\n",
" [ 0. 0. 0. 0. -0.41380839 -0.77906404\n",
" 0. -0.62420561]\n",
" [ 0. 0. 0. 0. -0.77906404 -0.47655051\n",
" -0.97722328 0. ]\n",
" [ 0. 0. 0. 0. 0. -0.97722328\n",
" -0.64027609 -0.82236485]\n",
" [ 0. 0. 0. 0. -0.62420561 0.\n",
" -0.82236485 -0.29956945]]\n",
"Shape of two body integral in spinorbital basis: (8, 8, 8, 8)\n"
]
}
],
"source": [
"# Example: generating XXZ Heisenberg model \n",
"# Returning electron integrals in a spatial orbital basis\n",
"# Assuming 4-fold symmetry\n",
"# Returning output as dense matrix\n",
"\n",
"# In the future we are planning to support different types of bonds for different atoms.\n",
"system = [('C1', 'Cl2', 1.1), ('Cl2', 'F3', 1.0),\n",
" ('F3', 'Si4', 1.0), ('Si4', 'C1', 1.0)]\n",
"ham = HamHub(system, u_onsite=np.array([1, 1,1,1]),\n",
" orbital_overlap=np.array([1,1,1,1]))\n",
"\n",
"e0 = ham.generate_zero_body_integral()\n",
"h1 = ham.generate_one_body_integral(dense=True, basis='spatial basis') \n",
"h2 = ham.generate_two_body_integral(dense=True, basis='spatial basis', sym=4)\n",
"\n",
"print(\"Zero energy: \", e0)\n",
"print(\"One body integrals in spatial basis: \\n\", h1)\n",
"print(\"Shape of two body integral in spatial basis: \", h2.shape)\n",
"print(\"-\"*60)\n",
"\n",
"# Example: generating XXZ Heisenberg model in spin orbital basis\n",
"# Assuming 4-fold symmetry\n",
"# Returning output as dense matrix\n",
"\n",
"h1 = ham.generate_one_body_integral(dense=True, basis='spinorbital basis')\n",
"h2 = ham.generate_two_body_integral(dense=True, basis='spinorbital basis', sym=4)\n",
"\n",
"print(\"One body integrals in spin basis: \\n\", h1)\n",
"print(\"Shape of two body integral in spinorbital basis: \", h2.shape)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Zero energy: 0\n",
"One body integrals in spatial basis: \n",
" [[ 0. -1. 0. 0. 0. -1.]\n",
" [-1. 0. -1. 0. 0. 0.]\n",
" [ 0. -1. 0. -1. 0. 0.]\n",
" [ 0. 0. -1. 0. -1. 0.]\n",
" [ 0. 0. 0. -1. 0. -1.]\n",
" [-1. 0. 0. 0. -1. 0.]]\n",
"Shape of two body integral in spatial basis: (6, 6, 6, 6)\n",
"------------------------------------------------------------\n",
"One body integrals in spin basis: \n",
" [[ 0. -1. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0.]\n",
" [-1. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [ 0. -1. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [ 0. 0. -1. 0. -1. 0. 0. 0. 0. 0. 0. 0.]\n",
" [ 0. 0. 0. -1. 0. -1. 0. 0. 0. 0. 0. 0.]\n",
" [-1. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0.]\n",
" [ 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. -1.]\n",
" [ 0. 0. 0. 0. 0. 0. -1. 0. -1. 0. 0. 0.]\n",
" [ 0. 0. 0. 0. 0. 0. 0. -1. 0. -1. 0. 0.]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. -1. 0.]\n",
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. -1.]\n",
" [ 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. -1. 0.]]\n",
"Shape of two body integral in spinorbital basis: (12, 12, 12, 12)\n"
]
}
],
"source": [
"# Example: generating 6 site Hubbard model \n",
"# Returning electron integrals in a spatial orbital basis\n",
"# Assuming 8-fold symmetry\n",
"# Returning output as dense matrix\n",
"\n",
"connectivity = np.array([[0, 1, 0, 0, 0, 1],\n",
" [1, 0, 1, 0, 0, 0],\n",
" [0, 1, 0, 1, 0, 0],\n",
" [0, 0, 1, 0, 1, 0],\n",
" [0, 0, 0, 1, 0, 1],\n",
" [1, 0, 0, 0, 1, 0]])\n",
"hubbard = HamHub(connectivity,\n",
" alpha=0, \n",
" beta=-1, \n",
" u_onsite=np.array([1, 1, 1, 1, 1, 1]))\n",
"\n",
"e0 = hubbard.generate_zero_body_integral()\n",
"h1 = hubbard.generate_one_body_integral(dense=True, basis='spatial basis') \n",
"h2 = hubbard.generate_two_body_integral(dense=True, basis='spatial basis', sym=8)\n",
"\n",
"print(\"Zero energy: \", e0)\n",
"print(\"One body integrals in spatial basis: \\n\", h1)\n",
"print(\"Shape of two body integral in spatial basis: \", h2.shape)\n",
"print(\"-\"*60)\n",
"\n",
"# Example: generating Hubbard model in spin orbital basis\n",
"# Assuming 8-fold symmetry\n",
"# Returning output as dense matrix\n",
"\n",
"h1 = hubbard.generate_one_body_integral(dense=True, basis='spinorbital basis')\n",
"h2 = hubbard.generate_two_body_integral(dense=True, basis='spinorbital basis', sym=4)\n",
"\n",
"print(\"One body integrals in spin basis: \\n\", h1)\n",
"print(\"Shape of two body integral in spinorbital basis: \", h2.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down
22 changes: 19 additions & 3 deletions moha/hamiltonians.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from typing import Union

from moha.rauk.rauk import assign_rauk_parameters
from moha.rauk.PariserParr import compute_overlap
from moha.rauk.PariserParr import compute_overlap, compute_gamma, compute_u
import warnings

warnings.simplefilter('ignore',
Expand Down Expand Up @@ -40,7 +40,9 @@ def __init__(
sym=1,
atom_dictionary=None,
bond_dictionary=None,
orbital_overlap=None
orbital_overlap=None,
affinity_dct=None,
Rxy_list=None
):
r"""
Initialize Pariser-Parr-Pople Hamiltonian.
Expand Down Expand Up @@ -103,6 +105,8 @@ def __init__(
self.bond_dictionary = bond_dictionary
self.atom_dictionary = atom_dictionary
self.orbital_overlap = orbital_overlap
self.affinity_dct = affinity_dct
self.Rxy_list = Rxy_list

def generate_zero_body_integral(self):
r"""Generate zero body integral.
Expand Down Expand Up @@ -217,12 +221,24 @@ def generate_two_body_integral(self, basis: str, dense: bool, sym=1):
Nv = 2 * n_sp
v = lil_matrix((Nv * Nv, Nv * Nv))

if self.u_onsite is None or (
self.affinity_dct is None and not isinstance(
self.connectivity, np.ndarray)):

self.u_onsite, self.Rxy_list = compute_u(
self.connectivity, self.atom_dictionary, self.affinity_dct)

if self.u_onsite is not None:
for p in range(n_sp):
i, j = convert_indices(Nv, p, p + n_sp, p, p + n_sp)
v[i, j] = self.u_onsite[p]

if self.gamma is not None:
if self.gamma is None and not isinstance(
self.connectivity, np.ndarray):
self.gamma = compute_gamma(
self.u_onsite, self.Rxy_list, self.connectivity)

elif self.gamma is not None:
if basis == "spinorbital basis" and \
self.gamma.shape != (n_sp, n_sp):
raise TypeError("Gamma matrix has wrong shape")
Expand Down
97 changes: 86 additions & 11 deletions moha/rauk/PariserParr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,18 +251,20 @@ def compute_overlap(
bond_dictionary[bond_key_forward] = beta_xy
bond_dictionary[bond_key_reverse] = beta_xy
else:
for tpl in connectivity:
for i, tpl in enumerate(connectivity):
atom1, atom2, dist = tpl[0], tpl[1], tpl[2]
atom1_name, _ = get_atom_type(atom1)
atom2_name, _ = get_atom_type(atom2)
bond_key_forward = ','.join([atom1_name, atom2_name])
bond_key_reverse = ','.join([atom2_name, atom1_name])
Sxy = orbital_overlap[bond_key_forward]

Sxy = orbital_overlap[i]

beta_xy = populate_PP_dct(
dist, atom1_name, atom2_name, ionization, Sxy)
bond_dictionary[bond_key_forward] = beta_xy
bond_dictionary[bond_key_reverse] = beta_xy
i += 1

one_body = build_one_body(
connectivity,
Expand All @@ -272,20 +274,93 @@ def compute_overlap(
return one_body


def calculate_gamma(Uxy_bar, Rxy):
def compute_gamma(u_onsite, Rxy_list, connectivity):
r"""
Calculate the gamma values for each pair of sites.
Parameters
----------
u_onsite (list of float): List of potential energies for sites.
Rxy_list (list of float): List of distances
connectivity (list of tuples): Each tuple contains indices of two
sites (atom1, atom2), indicating that a
gamma value should be computed for
this pair.
Returns
----------
np.ndarray: A matrix of computed gamma values for each pair.
Non-connected pairs have a gamma value of zero.
"""
Calculate the gamma value based on Uxy and Rxy.
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)

if site1 < num_sites and site2 < num_sites:
Ux = u_onsite[site1]
Uy = u_onsite[site2]
Rxy = Rxy_list[site1]
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

return gamma_matrix


def compute_u(connectivity, atom_dictionary, affinity_dictionary):
r"""
Calculate the onsite potential energy (U) for each site.
Parameters
----------
Uxy_bar (float): Represents the potential energy
Rxy (float): Represents the distance or a related measure.
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.
affinity_dictionary (dict): Dictionary mapping atom types to their
electron affinities.
Returns
----------
float: Computed gamma value based on the given parameters.
tuple: Contains two elements:
1. List of float: Onsite potential energies calculated
for each site based on their ionization energy
and electron affinity.
2. List of float: Distances corresponding to each
connectivity tuple.
"""
# Example formula, needs actual formula to be replaced here
# This is just a placeholder formula
gamma = Uxy_bar / (Uxy_bar * Rxy + np.exp(-1 / 2 * Uxy_bar**2 * Rxy ** 2))
return gamma
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 affinity_dictionary is None:
affinity_path = Path(__file__).parent / "affinity.json"
affinity_dictionary = json.load(open(affinity_path, "rb"))

u_onsite = []
Rxy_list = []

for tpl in connectivity:
atom1, atom2, dist = tpl[0], tpl[1], tpl[2]
atom1_name, _ = get_atom_type(atom1)
atom2_name, _ = get_atom_type(atom2)

ionization = atom_dictionary[atom1_name]
affinity = affinity_dictionary[atom1_name]
U_x = ionization - affinity

u_onsite.append(U_x)
Rxy_list.append(dist)

return u_onsite, Rxy_list
Loading

0 comments on commit 6c74a58

Please sign in to comment.