diff --git a/Final/Final - Tight-binding propagation method.ipynb b/Final/Final - Tight-binding propagation method.ipynb index ca2e07b..1082744 100644 --- a/Final/Final - Tight-binding propagation method.ipynb +++ b/Final/Final - Tight-binding propagation method.ipynb @@ -606,7 +606,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 17, "metadata": { "deletable": false, "nbgrader": { @@ -625,7 +625,16 @@ "source": [ "def hopping(i, j, n):\n", " \"\"\"\n", + " Calculates hopping matrix elements t_ij for sigma = 0.25 in a 1D\n", + " chain of n atoms at distance a = 1 from eachother.\n", " \n", + " Args:\n", + " i: origin site index\n", + " j: destination site index\n", + " n: number of atoms in the chain\n", + "\n", + " Returns:\n", + " Hopping parameter t_ij.\n", " \"\"\"\n", " \n", " sigma = .25\n", @@ -633,10 +642,11 @@ " \n", " # This 'infinity' is large enough, as the Gaussians decay quite quickly\n", " # away from the atomic positions, which we already saw in the overlap\n", - " # above.\n", - " # TODO: Remember https://www.wolframalpha.com/input?i=integrate+1%2F%28.25*sqrt%282pi%29%29*exp%28-.5*%28x%2F.25%29%5E2%29+for+x%3D-1..1.\n", - " infty = 1e2\n", - " h = 1e-6\n", + " # above. In fact, 99.7% of all probability mass is under the integral\n", + " # for x radius of 3*sigma from the centers x_i.\n", + " # TODO: Make integration bounds relative to n and a, and maybe sigma.\n", + " infty = 5e1\n", + " h = 1e-5\n", " x = np.arange(-infty, infty, h)\n", " \n", " def V(x):\n", @@ -653,21 +663,12 @@ " #V = lambda x: np.sum( -1/( np.abs(np.subtract.outer(x, positions)) + 0.001 ), axis=1 )\n", " \n", " integrand = lambda x: atomic_basis(x, positions[i], sigma)*V(x)*atomic_basis(x, positions[j], sigma)\n", - " return integrate(integrand, x)\n", - "\n", - "# Yann had\n", - "# 1) a plot of the peaks as before\n", - "# 2) a weird line under it with peaks downward at these x\n", - "# and slowly rising energy values outside the chain\n", - "# 3)\n", - "# For i = 1 ...\n", - "# t_{i,i-1} = -0.1203\n", - "# t_{i,i+1} = ..." + " return integrate(integrand, x)" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 27, "metadata": { "deletable": false, "nbgrader": { @@ -687,15 +688,88 @@ "name": "stdout", "output_type": "stream", "text": [ - "-3.08808170797426e-06 -0.13849169656144708\n" + "For i = 0 ...\n", + "\tt_{i,i+1} = -0.13849171527765194\n", + "\tt_{i,i+2} = -3.0880816778755327e-06\n", + "\tt_{i,i+3} = -1.8833562169250255e-15\n", + "\n", + "For i = 1 ...\n", + "\tt_{i,i-1} = -0.13849171527765194\n", + "\tt_{i,i+1} = -0.14871538196685558\n", + "\tt_{i,i+2} = -3.1307050912510583e-06\n", + "\tt_{i,i+3} = -1.9456306931892115e-15\n", + "\n", + "For i = 2 ...\n", + "\tt_{i,i-1} = -0.14871538196685555\n", + "\tt_{i,i+1} = -0.1536327589565012\n", + "\tt_{i,i-2} = -3.0880816778755327e-06\n", + "\tt_{i,i+2} = -3.152251398227753e-06\n", + "\tt_{i,i+3} = -1.976347912754562e-15\n", + "\n", + "For i = 3 ...\n", + "\tt_{i,i-1} = -0.1536327589565012\n", + "\tt_{i,i+1} = -0.15605828154196247\n", + "\tt_{i,i-2} = -3.1307050912510583e-06\n", + "\tt_{i,i+2} = -3.1616580016478533e-06\n", + "\tt_{i,i-3} = -1.8833562169250255e-15\n", + "\tt_{i,i+3} = -1.98575211950378e-15\n", + "\n", + "For i = 4 ...\n", + "\tt_{i,i-1} = -0.15605828154196244\n", + "\tt_{i,i+1} = -0.15680086554505168\n", + "\tt_{i,i-2} = -3.152251398227753e-06\n", + "\tt_{i,i+2} = -3.161664327360904e-06\n", + "\tt_{i,i-3} = -1.9456306931892115e-15\n", + "\tt_{i,i+3} = -1.976348151505985e-15\n", + "\n", + "For i = 5 ...\n", + "\tt_{i,i-1} = -0.15680086554505168\n", + "\tt_{i,i+1} = -0.156058300394482\n", + "\tt_{i,i-2} = -3.1616580016478533e-06\n", + "\tt_{i,i+2} = -3.152251395817713e-06\n", + "\tt_{i,i-3} = -1.9763479127545625e-15\n", + "\tt_{i,i+3} = -1.945630454073864e-15\n", + "\n", + "For i = 6 ...\n", + "\tt_{i,i-1} = -0.15605830039448199\n", + "\tt_{i,i+1} = -0.153632740075246\n", + "\tt_{i,i-2} = -3.161664327360904e-06\n", + "\tt_{i,i+2} = -3.1306987607180875e-06\n", + "\tt_{i,i-3} = -1.9857521195037804e-15\n", + "\tt_{i,i+3} = -1.8833562165610968e-15\n", + "\n", + "For i = 7 ...\n", + "\tt_{i,i-1} = -0.153632740075246\n", + "\tt_{i,i+1} = -0.14871538193811829\n", + "\tt_{i,i-2} = -3.152251395817713e-06\n", + "\tt_{i,i+2} = -3.088087998768686e-06\n", + "\tt_{i,i-3} = -1.976348151505985e-15\n", + "\n", + "For i = 8 ...\n", + "\tt_{i,i-1} = -0.14871538193811829\n", + "\tt_{i,i+1} = -0.1384917341014353\n", + "\tt_{i,i-2} = -3.1306987607180875e-06\n", + "\tt_{i,i-3} = -1.945630454073864e-15\n", + "\n", + "For i = 9 ...\n", + "\tt_{i,i-1} = -0.1384917341014353\n", + "\tt_{i,i-2} = -3.088087998768686e-06\n", + "\tt_{i,i-3} = -1.8833562165610968e-15\n", + "\n" ] } ], "source": [ - "print(hopping(0, 2, 10), hopping(0, 1, 10))\n", + "n = 10\n", "\n", - "# YOUR CODE HERE\n", - "#raise NotImplementedError()" + "for i in range(n):\n", + " print(\"For i =\", i, \"...\")\n", + " for r in range(1, 4):\n", + " if i - r >= 0:\n", + " print(\"\\tt_{{i,i-{}}} = {}\".format(r, hopping(i, i - r, n)))\n", + " if i + r < n:\n", + " print(\"\\tt_{{i,i+{}}} = {}\".format(r, hopping(i, i + r, n)))\n", + " print()" ] }, { @@ -723,7 +797,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 54, "metadata": { "deletable": false, "nbgrader": { @@ -740,10 +814,30 @@ }, "outputs": [], "source": [ - "# YOUR CODE HERE\n", - "raise NotImplementedError()\n", - "\n", - "# Yann said we just have to be able to solve tri-diagonal matrix diagonalization, but more general is ok." + "def QREig(T, eps=1e-6, k_max=10000):\n", + " \"\"\"\n", + " Follows the method of the QR decomposition based diagonalization routine\n", + " for tridiagonal matrices. The matrix T is diagonalized, resulting in\n", + " all diagonal elements being an eigenvalue.\n", + " \n", + " Args:\n", + " T: a tridiagonaliz matrix.\n", + " eps: the desired accuracy.\n", + " k_max: maximum number of iterations after which to cut off\n", + " \n", + " Returns:\n", + " A one dimensional array with the eigenvalues of the matrix T.\n", + " \"\"\"\n", + " \n", + " e = eps + 1\n", + " k = 0\n", + " while e > eps and k < k_max:\n", + " k += 1\n", + " Q, R = np.linalg.qr(T)\n", + " T = np.matmul(R,Q)\n", + " e = np.sum(np.abs(np.diag(T, k=1)))\n", + " print(k)\n", + " return np.diag(T)" ] }, { @@ -769,7 +863,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 75, "metadata": { "deletable": false, "nbgrader": { @@ -784,10 +878,46 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10000\n", + "1000\n", + "[-3.35311057 0.14162732 3. 4.21148325]\n", + "10000\n" + ] + } + ], "source": [ - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "def test_QREig():\n", + " # Test case one\n", + " T = np.array([\n", + " [1,4,0,0],\n", + " [3,4,1,0],\n", + " [0,2,3,4],\n", + " [0,0,1,3]\n", + " ])\n", + " # Eigenvalues are roots of λ^4 - 11*λ^3 + 25*λ^2 + 31*λ - 46.\n", + " eigenvalues_of_T = np.array([-1.45350244, 1., 4.65531023, 6.79819221])\n", + " assert np.allclose(np.sort(QREig(T)), eigenvalues_of_T)\n", + " #print(QREig(T))\n", + " \n", + " # Test case two\n", + " T = np.array([\n", + " [1,4,0,0],\n", + " [3,0,1,0],\n", + " [0,2,0,4],\n", + " [0,0,0,3]\n", + " ])\n", + " #eigenvalues_of_T = np.array([-3.35311, .141627, 3. , 4.21148])\n", + " eigenvalues_of_T = np.sort(eig(T)[0])\n", + " #print(QREig(T+np.eye(len(T)), k_max=1000000)-1)\n", + " print(np.sort(QREig(T, k_max=1000)))\n", + " assert np.allclose(np.sort(QREig(T)), eigenvalues_of_T)\n", + "\n", + "test_QREig()" ] }, { @@ -815,7 +945,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 77, "metadata": { "deletable": false, "nbgrader": { @@ -833,8 +963,14 @@ "outputs": [], "source": [ "def TBHamiltonian(n):\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()" + " \"\"\"\n", + " \n", + " \"\"\"\n", + " \n", + " # TODO: Comment on the weird 20% differences in hopping parameters.\n", + " \n", + " t = hopping(i, j, n)\n", + " H = np.eye(n) + np.eye(n, -1) + np.eye(n, 1)" ] }, { @@ -856,6 +992,9 @@ }, "outputs": [], "source": [ + "for n in [10, 20, 40, 80, 100]:\n", + " pass\n", + "\n", "# Do your plotting here ...\n", "\n", "# YOUR CODE HERE\n",