Skip to content

Commit 69dee4a

Browse files
ilfreddyilfreddy
authored
Add files via upload
1 parent 4aa0194 commit 69dee4a

File tree

2 files changed

+383
-0
lines changed

2 files changed

+383
-0
lines changed

‎content/notebooks/H-atom-4c.ipynb

+299
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,299 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "90e862bb-783a-4f97-a227-8f6e0b4012b8",
6+
"metadata": {},
7+
"source": [
8+
"# The Dirac equation for the Hydrogen atom\n",
9+
"\n",
10+
"In this notebook we will illustrate how one can solve the Dirac equation for a Hydrogen atom using the Multiwavelet framework provided by VAMPyR\n",
11+
"\n",
12+
"The Dirac equation can be coincisely written as follows:\n",
13+
"\n",
14+
"\\begin{equation}\n",
15+
"(\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p} + V) \\phi = \\epsilon \\phi \n",
16+
"\\end{equation}\n",
17+
"\n",
18+
"where $\\phi = (\\phi^{L\\alpha}, \\phi^{L\\beta}, \\phi^{S\\alpha}, \\phi^{S\\beta})^t$ represents a 4-component spinor, $\\boldsymbol{\\alpha} = \n",
19+
"\\begin{pmatrix}\n",
20+
"0_2 & \\boldsymbol{\\sigma} \\\\\n",
21+
"\\boldsymbol{\\sigma} & 0_2 & \n",
22+
"\\end{pmatrix}\n",
23+
"$ and\n",
24+
"$\\beta = \n",
25+
"\\begin{pmatrix}\n",
26+
"1_2 & 0_2 \\\\\n",
27+
"0_2 & -1_2\n",
28+
"\\end{pmatrix}\n",
29+
"$ are the 4x4 Dirac matrices, $\\boldsymbol{\\sigma}$ is a cartesian vector collecting the three 2x2 Pauli matrices, $\\mathbf{p}$ is the momentum operator, $c$ is the speed of light, $m$ is the electron mass and $V$ is the nuclear potential.\n",
30+
"\n",
31+
"As for the non-relativistic case, the equation is first rewritten in its integral formulation:\n",
32+
"$$\\phi = \\frac{1}{2mc^2}(\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p} + \\epsilon) \\left[ G_\\mu \\star (V \\psi) \\right]$$\n",
33+
"\n",
34+
"where $G_\\mu(x) = \\frac{e^{-\\mu |x|}}{4 \\pi |x|}$ is the Helmholtz Green's kernel and $\\mu = \\sqrt{\\frac{m^2c^4-\\epsilon}{mc^2}}$. An initial guess can be obtained by taking a Slater orbital or a Gaussian function for the $\\psi^{L\\alpha}$ component and then applying the restricted kinetic balance:\n",
35+
"\n",
36+
"$$\n",
37+
"\\begin{pmatrix}\n",
38+
"\\phi^{S\\alpha} \\\\\n",
39+
"\\phi^{S\\beta}\n",
40+
"\\end{pmatrix}\n",
41+
"= \\frac{1}{2c}\\boldsymbol{\\sigma} \\cdot \\mathbf{p} \n",
42+
"\\begin{pmatrix}\n",
43+
"\\phi^{L\\alpha} \\\\\n",
44+
"0\n",
45+
"\\end{pmatrix}\n",
46+
"$$\n",
47+
"The guess obtained is then plugged into the integral form of the Dirac equation, which can then be iterated until convergence"
48+
]
49+
},
50+
{
51+
"cell_type": "markdown",
52+
"id": "54c75006-0358-42c5-9075-37a5b8a909ba",
53+
"metadata": {},
54+
"source": [
55+
"We start by loading the relevant packages: the 3d version of `vampyr`, `numpy`, the `complex_function` class which deals with complex funtions and the `orbital` class which deals with 4-component spinors. Each complex function is handled as a pair of `function_tree`s and each spinor is handled as a 4c vector of complex functions. The `nuclear_potential` package is self-explanatory"
56+
]
57+
},
58+
{
59+
"cell_type": "code",
60+
"execution_count": 13,
61+
"id": "37e0e6b2-886e-4415-8612-d2a652f24c4f",
62+
"metadata": {},
63+
"outputs": [],
64+
"source": [
65+
"from vampyr import vampyr3d as vp\n",
66+
"from orbital4c import orbital as orb\n",
67+
"from orbital4c import nuclear_potential as nucpot\n",
68+
"from orbital4c import complex_fcn as cf\n",
69+
"import numpy as np"
70+
]
71+
},
72+
{
73+
"cell_type": "markdown",
74+
"id": "7fbe4360-9936-47d4-bb6f-d51bb3b34b69",
75+
"metadata": {},
76+
"source": [
77+
"As a reference, the exact Dirac energy for the ground state Hydrogen atom is computed in the function below."
78+
]
79+
},
80+
{
81+
"cell_type": "code",
82+
"execution_count": 14,
83+
"id": "0ef3d827-810b-4411-8ba3-bd63271033e4",
84+
"metadata": {},
85+
"outputs": [
86+
{
87+
"name": "stdout",
88+
"output_type": "stream",
89+
"text": [
90+
"Exact Energy -0.5000066565989982\n"
91+
]
92+
}
93+
],
94+
"source": [
95+
"def analytic_1s(light_speed, n, k, Z):\n",
96+
" alpha = 1/light_speed\n",
97+
" gamma = orb.compute_gamma(k,Z,alpha)\n",
98+
" tmp1 = n - np.abs(k) + gamma\n",
99+
" tmp2 = Z * alpha / tmp1\n",
100+
" tmp3 = 1 + tmp2**2\n",
101+
" return light_speed**2 / np.sqrt(tmp3)\n",
102+
"\n",
103+
"light_speed = 137.03599913900001\n",
104+
"alpha = 1/light_speed\n",
105+
"k = -1\n",
106+
"l = 0\n",
107+
"n = 1\n",
108+
"m = 0.5\n",
109+
"Z = 1\n",
110+
"atom = \"H\"\n",
111+
"\n",
112+
"energy_1s = analytic_1s(light_speed, n, k, Z)\n",
113+
"print('Exact Energy',energy_1s - light_speed**2, flush = True)"
114+
]
115+
},
116+
{
117+
"cell_type": "markdown",
118+
"id": "6277d776-2236-496a-bdd5-41929d11ac3d",
119+
"metadata": {},
120+
"source": [
121+
"The `MultiResolutionAnalysis` object defining the simulation box is constructed and passed to the classes for complex functions and spinors"
122+
]
123+
},
124+
{
125+
"cell_type": "code",
126+
"execution_count": 15,
127+
"id": "5ac41574-6f43-411d-841f-f8fdf0749b66",
128+
"metadata": {},
129+
"outputs": [],
130+
"source": [
131+
"mra = vp.MultiResolutionAnalysis(box=[-30,30], order=6)\n",
132+
"prec = 1.0e-4\n",
133+
"origin = [0.1, 0.2, 0.3] # origin moved to avoid placing the nuclar charge on a node\n",
134+
"\n",
135+
"orb.orbital4c.light_speed = light_speed\n",
136+
"orb.orbital4c.mra = mra\n",
137+
"cf.complex_fcn.mra = mra"
138+
]
139+
},
140+
{
141+
"cell_type": "markdown",
142+
"id": "d5abdbe4-ba5b-49d3-ae7e-febab5ced4fe",
143+
"metadata": {},
144+
"source": [
145+
"We construct a starting guess by taking a simple Gaussian function and initialize the real part of the $\\phi^{L\\alpha}$ component of the spinor with it. Thereafter the restricted kinetic balance is employed. This is implemented in the `init_small_components` method of the `orbital` class"
146+
]
147+
},
148+
{
149+
"cell_type": "code",
150+
"execution_count": 16,
151+
"id": "95826623-b759-4318-afc5-e0c225cb5f1d",
152+
"metadata": {},
153+
"outputs": [],
154+
"source": [
155+
"a_coeff = 3.0\n",
156+
"b_coeff = np.sqrt(a_coeff/np.pi)**3\n",
157+
"gauss = vp.GaussFunc(b_coeff, a_coeff, origin)\n",
158+
"gauss_tree = vp.FunctionTree(mra)\n",
159+
"vp.advanced.build_grid(out=gauss_tree, inp=gauss)\n",
160+
"vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss)\n",
161+
"gauss_tree.normalize()\n",
162+
"\n",
163+
"spinor_H = orb.orbital4c()\n",
164+
"La_comp = cf.complex_fcn()\n",
165+
"La_comp.copy_fcns(real = gauss_tree)\n",
166+
"spinor_H.copy_components(La = La_comp)\n",
167+
"spinor_H.init_small_components(prec/10)\n",
168+
"spinor_H.normalize()"
169+
]
170+
},
171+
{
172+
"cell_type": "markdown",
173+
"id": "f6a40ddd-7933-46ee-8128-96caf89f3682",
174+
"metadata": {},
175+
"source": [
176+
"The nuclear potential is defined and projected onto the `V_tree`"
177+
]
178+
},
179+
{
180+
"cell_type": "code",
181+
"execution_count": 17,
182+
"id": "4f0cda6f-217c-4941-8ce5-fa52576c2d98",
183+
"metadata": {},
184+
"outputs": [],
185+
"source": [
186+
"Peps = vp.ScalingProjector(mra, prec)\n",
187+
"f = lambda x: nucpot.coulomb_HFYGB(x, origin, Z, prec)\n",
188+
"V_tree = Peps(f)\n",
189+
"\n",
190+
"default_der = 'BS'"
191+
]
192+
},
193+
{
194+
"cell_type": "markdown",
195+
"id": "1678bc97-d85d-4a48-9aea-13901631534d",
196+
"metadata": {},
197+
"source": [
198+
"The orbital is optimized by iterating the integral version of the Dirac equation as follows:\n",
199+
"1. Application of the Dirac Hamiltonian $f^n = \\hat{h}_D \\phi^n = (\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p}) \\phi^n$\n",
200+
"2. Application of the potnetial operator $g^n = \\hat{V} \\phi^n$\n",
201+
"3. Sum $h^n = f^n + g^n$\n",
202+
"4. Expectation value of the energy $\\left\\langle \\phi^n | h^n \\right\\rangle$\n",
203+
"5. Calculation of the Helmholtz parameter $\\mu$\n",
204+
"6. Convolution with the Helmholtz kernel $t^n = G_\\mu \\star g^n$\n",
205+
"7. application of the shifted Dirac Hamiltonian $\\tilde{\\phi}^{n+1} = (\\hat{h}_D + \\epsilon) t^n$\n",
206+
"8. normalization of the new iterate\n",
207+
"9. calculation of the change in the orbital $\\delta \\phi^n = \\phi^{n+1}-\\phi^n$\n",
208+
"\n",
209+
"Once the orbital error is below the requested threshold the iteration is interrupted and the final energy is computed."
210+
]
211+
},
212+
{
213+
"cell_type": "code",
214+
"execution_count": 18,
215+
"id": "40a443b6-33f4-4771-8b06-ccd12c9da35d",
216+
"metadata": {},
217+
"outputs": [
218+
{
219+
"name": "stdout",
220+
"output_type": "stream",
221+
"text": [
222+
"Energy -0.14180720307558659\n",
223+
"Error 0.3648353655184581\n",
224+
"Energy -0.49611934473068686\n",
225+
"Error 0.0024839139737050653\n",
226+
"Energy -0.4997319277244969\n",
227+
"Error 0.00023304055836475225\n",
228+
"Energy -0.49998214000152075\n",
229+
"Error 2.3987883028320866e-05\n",
230+
"Final Energy -0.5000042447172746\n",
231+
"Exact Energy -0.5000066565989982\n",
232+
"Difference -2.411881723674014e-06\n"
233+
]
234+
}
235+
],
236+
"source": [
237+
"orbital_error = 1\n",
238+
"while orbital_error > prec:\n",
239+
" # 1. Application of the Dirac Hamiltonian\n",
240+
" hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)\n",
241+
" # 2. Application of the potnetial operator\n",
242+
" v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)\n",
243+
" # 3. Sum\n",
244+
" add_psi = hd_psi + v_psi\n",
245+
" # 4. Expectation value of the energy\n",
246+
" energy = (spinor_H.dot(add_psi)).real\n",
247+
" print('Energy',energy-light_speed**2)\n",
248+
" # 5. Calculation of the Helmholtz parameter\n",
249+
" mu = orb.calc_dirac_mu(energy, light_speed)\n",
250+
" # 6. Convolution with the Helmholtz kernel\n",
251+
" tmp = orb.apply_helmholtz(v_psi, mu, prec)\n",
252+
" tmp.crop(prec/10)\n",
253+
" # 7. application of the shifted Dirac Hamiltonian\n",
254+
" new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = default_der) \n",
255+
" new_orbital.crop(prec/10)\n",
256+
" # 8. normalization of the new iterate\n",
257+
" new_orbital.normalize()\n",
258+
" delta_psi = new_orbital - spinor_H\n",
259+
" # 9. calculation of the change in the orbital\n",
260+
" orbital_error = (delta_psi.dot(delta_psi)).real\n",
261+
" print('Error',orbital_error, flush = True)\n",
262+
" spinor_H = new_orbital\n",
263+
"\n",
264+
"# Computing the final energy\n",
265+
"hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)\n",
266+
"v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)\n",
267+
"add_psi = hd_psi + v_psi\n",
268+
"energy = (spinor_H.dot(add_psi)).real\n",
269+
"print('Final Energy',energy - light_speed**2)\n",
270+
"\n",
271+
"energy_1s = analytic_1s(light_speed, n, k, Z)\n",
272+
"\n",
273+
"print('Exact Energy',energy_1s - light_speed**2)\n",
274+
"print('Difference',energy_1s - energy)"
275+
]
276+
}
277+
],
278+
"metadata": {
279+
"kernelspec": {
280+
"display_name": "Python 3 (ipykernel)",
281+
"language": "python",
282+
"name": "python3"
283+
},
284+
"language_info": {
285+
"codemirror_mode": {
286+
"name": "ipython",
287+
"version": 3
288+
},
289+
"file_extension": ".py",
290+
"mimetype": "text/x-python",
291+
"name": "python",
292+
"nbconvert_exporter": "python",
293+
"pygments_lexer": "ipython3",
294+
"version": "3.10.13"
295+
}
296+
},
297+
"nbformat": 4,
298+
"nbformat_minor": 5
299+
}

‎content/notebooks/H-test.py

+84
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
from vampyr import vampyr3d as vp
2+
from orbital4c import orbital as orb
3+
from orbital4c import nuclear_potential as nucpot
4+
from orbital4c import complex_fcn as cf
5+
import numpy as np
6+
7+
def analytic_1s(light_speed, n, k, Z):
8+
alpha = 1/light_speed
9+
gamma = orb.compute_gamma(k,Z,alpha)
10+
tmp1 = n - np.abs(k) + gamma
11+
tmp2 = Z * alpha / tmp1
12+
tmp3 = 1 + tmp2**2
13+
return light_speed**2 / np.sqrt(tmp3)
14+
15+
light_speed = 137.03599913900001
16+
alpha = 1/light_speed
17+
k = -1
18+
l = 0
19+
n = 1
20+
m = 0.5
21+
Z = 1
22+
atom = "H"
23+
24+
energy_1s = analytic_1s(light_speed, n, k, Z)
25+
print('Exact Energy',energy_1s - light_speed**2, flush = True)
26+
27+
mra = vp.MultiResolutionAnalysis(box=[-30,30], order=6)
28+
prec = 1.0e-4
29+
origin = [0.1, 0.2, 0.3] # origin moved to avoid placing the nuclar charge on a node
30+
31+
orb.orbital4c.light_speed = light_speed
32+
orb.orbital4c.mra = mra
33+
cf.complex_fcn.mra = mra
34+
35+
a_coeff = 3.0
36+
b_coeff = np.sqrt(a_coeff/np.pi)**3
37+
gauss = vp.GaussFunc(b_coeff, a_coeff, origin)
38+
gauss_tree = vp.FunctionTree(mra)
39+
vp.advanced.build_grid(out=gauss_tree, inp=gauss)
40+
vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss)
41+
gauss_tree.normalize()
42+
43+
spinor_H = orb.orbital4c()
44+
La_comp = cf.complex_fcn()
45+
La_comp.copy_fcns(real = gauss_tree)
46+
spinor_H.copy_components(La = La_comp)
47+
spinor_H.init_small_components(prec/10)
48+
spinor_H.normalize()
49+
50+
Peps = vp.ScalingProjector(mra, prec)
51+
f = lambda x: nucpot.coulomb_HFYGB(x, origin, Z, prec)
52+
V_tree = Peps(f)
53+
54+
print('V', V_tree)
55+
56+
default_der = 'BS'
57+
58+
orbital_error = 1
59+
while orbital_error > prec:
60+
hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)
61+
v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
62+
add_psi = hd_psi + v_psi
63+
energy = (spinor_H.dot(add_psi)).real
64+
print('Energy',energy-light_speed**2)
65+
mu = orb.calc_dirac_mu(energy, light_speed)
66+
tmp = orb.apply_helmholtz(v_psi, mu, prec)
67+
tmp.crop(prec/10)
68+
new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = default_der)
69+
new_orbital.crop(prec/10)
70+
new_orbital.normalize()
71+
delta_psi = new_orbital - spinor_H
72+
orbital_error = (delta_psi.dot(delta_psi)).real
73+
print('Error',orbital_error, flush = True)
74+
spinor_H = new_orbital
75+
76+
hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)
77+
v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
78+
add_psi = hd_psi + v_psi
79+
energy = (spinor_H.dot(add_psi)).real
80+
print('Final Energy',energy - light_speed**2)
81+
energy_1s = analytic_1s(light_speed, n, k, Z)
82+
83+
print('Exact Energy',energy_1s - light_speed**2)
84+
print('Difference',energy_1s - energy)

0 commit comments

Comments
 (0)