Skip to content

Commit

Permalink
Merge pull request #126 from theochem/giovanni-br-issue-116-fix
Browse files Browse the repository at this point in the history 
Giovanni br issue 116 fix
  • Loading branch information
@RichRick1
RichRick1 authored Jul 2, 2024
2 parents 9d0b5be + 8af1393 commit 5ae0eb3
Show file tree
Hide file tree
Showing 5 changed files with 296 additions and 17 deletions.
164 changes: 158 additions & 6 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. ]\n",
" [-0.77906404 -0.47655051 -2.93166983 0. ]\n",
" [ 0. -2.93166983 -0.64027609 -3.28945939]\n",
" [ 0. 0. -3.28945939 -0.29956945]]\n"
]
}
],
Expand All @@ -237,18 +237,170 @@
"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([[0, 1,0,0],[0, 0,3,0], [0, 0,0,4], [0, 0,0,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. ]\n",
" [-0.77906404 -0.47655051 -2.93166983 0. ]\n",
" [ 0. -2.93166983 -0.64027609 -3.28945939]\n",
" [ 0. 0. -3.28945939 -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. 0. 0.\n",
" 0. 0. ]\n",
" [-0.77906404 -0.47655051 -2.93166983 0. 0. 0.\n",
" 0. 0. ]\n",
" [ 0. -2.93166983 -0.64027609 -3.28945939 0. 0.\n",
" 0. 0. ]\n",
" [ 0. 0. -3.28945939 -0.29956945 0. 0.\n",
" 0. 0. ]\n",
" [ 0. 0. 0. 0. -0.41380839 -0.77906404\n",
" 0. 0. ]\n",
" [ 0. 0. 0. 0. -0.77906404 -0.47655051\n",
" -2.93166983 0. ]\n",
" [ 0. 0. 0. 0. 0. -2.93166983\n",
" -0.64027609 -3.28945939]\n",
" [ 0. 0. 0. 0. 0. 0.\n",
" -3.28945939 -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([[0, 1, 0, 0],[0, 0,3,0], [0, 0,0,4], [0, 0,0,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": {},
"source": [
"## References:\n",
"\n",
"[1]A. Orbital Interaction Theory of Organic Chemistry. [s.l.] John Wiley & Sons, 2004."
"[1]A. Orbital Interaction Theory of Organic Chemistry. [s.l.] John Wiley & Sons, 2004.\n",
"\n",
"[2] https://en.wikipedia.org/wiki/Electron_affinity_(data_page)"
]
}
],
Expand Down
28 changes: 27 additions & 1 deletion moha/hamiltonians.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,21 @@ def generate_one_body_integral(self, basis: str, dense: bool):
# check if elements in connectivity matrix are float
elif np.all([isinstance(elem[2], float)
for elem in self.connectivity]):
if self.orbital_overlap is not None:
# Assuming self.orbital_overlap should be a 2D array
if hasattr(
self,
'orbital_overlap') and self.orbital_overlap.ndim == 2:
if self.n_sites != self.orbital_overlap.shape[
0] or self.n_sites != self.orbital_overlap.shape[
1]:
raise TypeError("Overlap matrix has wrong dimensions")
else:
raise ValueError(
"orbital_overlap is not properly initialized ")

# Proceed with the rest of the function

one_body_term = compute_overlap(
self.connectivity,
self.atom_dictionary,
Expand Down Expand Up @@ -216,13 +231,24 @@ def generate_two_body_integral(self, basis: str, dense: bool, sym=1):
n_sp = self.n_sites
Nv = 2 * n_sp
v = lil_matrix((Nv * Nv, Nv * Nv))
if self.u_onsite is None:
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.Rxy_list = compute_u(
self.connectivity, self.atom_dictionary, self.affinity_dct)
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
10 changes: 6 additions & 4 deletions moha/rauk/PariserParr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,12 +252,14 @@ def compute_overlap(
bond_dictionary[bond_key_reverse] = beta_xy
else:
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)
atom1, atom2, dist = tpl
atom1_name, site1 = get_atom_type(atom1)
atom2_name, site2 = 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]
site1, site2 = site1 - 1, site2 - 1

Sxy = orbital_overlap[site1, site2]

beta_xy = populate_PP_dct(
dist, atom1_name, atom2_name, ionization, Sxy)
Expand Down
91 changes: 91 additions & 0 deletions moha/rauk/affinity.Json
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
{
"H": 0.754195,
"He": -0.5,
"Li": 0.618049,
"Be": -0.5,
"B": 0.279723,
"C": 1.2621226,
"N": -0.07,
"O": 1.46111297,
"F": 3.4011898,
"Ne": -1.2,
"Na": 0.547926,
"Mg": -0.4,
"Al": 0.43283,
"Si": 1.3895212,
"P": 0.746609,
"S": 2.0771042,
"Cl": 3.612725,
"Ar": -1.0,
"K": 0.501459,
"Ca": 0.02455,
"Sc": 0.17938,
"Ti": 0.07554,
"V": 0.52766,
"Cr": 0.675928,
"Mn": -0.5,
"Fe": 0.153236,
"Co": 0.662255,
"Ni": 1.15716,
"Cu": 1.23578,
"Zn": -0.6,
"Ga": 0.301166,
"Ge": 1.2326764,
"As": 0.8048,
"Se": 2.0206047,
"Br": 3.363588,
"Kr": -1.0,
"Rb": 0.485916,
"Sr": 0.05206,
"Y": 0.31129,
"Zr": 0.43328,
"Nb": 0.91740,
"Mo": 0.74723,
"Tc": 0.55,
"Ru": 1.04627,
"Rh": 1.14289,
"Pd": 0.56214,
"Ag": 1.30447,
"Cd": -0.7,
"In": 0.38392,
"Sn": 1.1120702,
"Sb": 1.047401,
"Te": 1.9708757,
"I": 3.059052,
"Xe": -0.8,
"Cs": 0.4715983,
"Ba": 0.14462,
"La": 0.557546,
"Ce": 0.600160,
"Pr": 0.10923,
"Nd": 0.09749,
"Pm": 0.129,
"Sm": 0.162,
"Eu": 0.116,
"Gd": 0.212,
"Tb": 0.13131,
"Dy": 0.015,
"Ho": 0.338,
"Er": 0.312,
"Tm": 1.029,
"Yb": -0.02,
"Lu": 0.2388,
"Hf": 0.1780,
"Ta": 0.328859,
"W": 0.81626,
"Re": 0.060396,
"Os": 1.077661,
"Ir": 1.564057,
"Pt": 2.12510,
"Au": 2.308610,
"Hg": -0.5,
"Tl": 0.320053,
"Pb": 0.356721,
"Bi": 0.942362,
"Po": 1.4,
"At": 2.41578,
"Rn": -0.7,
"Fr": 0.486,
"Ra": 0.1,
"Ac": 0.35
}
20 changes: 14 additions & 6 deletions moha/test/test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ def test_hub2():
"""
hubbard = HamHub([("C1", "C2", 1)],
alpha=0, beta=-1, u_onsite=np.array([1, 1]), sym=1)

ecore = hubbard.generate_zero_body_integral()
h = hubbard.generate_one_body_integral(basis='spatial basis', dense=True)
v = hubbard.generate_two_body_integral(sym=1,
Expand Down Expand Up @@ -41,10 +42,17 @@ def test_hub4():
nsites = np.linspace(2, 8, 4).astype(int)
for nsite in nsites:
nelec = nsite // 2
hubbard = HamPPP([(f"C{i}", f"C{i + 1}", 1) for i in range(1, nsite)] +
[(f"C{nsite}", f"C{1}", 1)],
alpha=0, beta=-1,
u_onsite=np.array([1 for i in range(nsite + 1)]))
hubbard = HamPPP([(f"C{i}",
f"C{i + 1}",
1) for i in range(1,
nsite)] + [(f"C{nsite}",
f"C{1}",
1)],
alpha=0,
beta=-1,
u_onsite=np.array([1 for i in range(nsite)]),
gamma=np.zeros((nsite,
nsite)))
ecore = hubbard.generate_zero_body_integral()
h = hubbard.generate_one_body_integral(basis='spatial basis',
dense=True)
Expand Down Expand Up @@ -73,8 +81,8 @@ def test_ethylene():
"""
a = -11.26
b = -1.45
hubbard = HamPPP([("C1", "C2", 1)], alpha=a, beta=b,
gamma=None, charges=None, sym=None)
hubbard = HamPPP([("C1", "C2", 1)], alpha=a, beta=b, gamma=np.zeros(
(2, 2)), charges=None, sym=None, u_onsite=[0, 0])
ecore = hubbard.generate_zero_body_integral()
h = hubbard.generate_one_body_integral(basis='spinorbital basis',
dense=True)
Expand Down

0 comments on commit 5ae0eb3

Please sign in to comment.