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

Issue #116 Fixed #119

Merged
merged 2 commits into from
Jun 9, 2024
Merged

Issue #116 Fixed #119

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
21df587
#116 PariserParr module
giovanni-br Jun 8, 2024
a864ff5
Merge branch 'theochem:main' into issue-117-fix
giovanni-br Jun 9, 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
160 changes: 155 additions & 5 deletions examples/rauk.ipynb
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading