diff --git a/README.md b/README.md index 8e979d6..9b744f9 100644 --- a/README.md +++ b/README.md @@ -33,6 +33,9 @@ export C_INCLUDE_PATH+=:$PETSC_DIR/include/:$SLEPC_DIR/include:$PETSC_DIR/arch-l make ex1 ``` +[![asciicast](https://asciinema.org/a/Ng3tSNDoWBkV5C9ZYvbxCW43B.png)](https://asciinema.org/a/Ng3tSNDoWBkV5C9ZYvbxCW43B) + + _Using DEHam_ --------------- @@ -53,16 +56,42 @@ true # Restrict the hole to the 1'st (i.e. half of natom) Family of states. .1430,-0.20,0.0000 # The three types of links this line gives J, K .1430,-0.20,0.0000 # -1.00,0.0,0.00 # This line gives t -1 # Currently unused (Perhaps can be used for potential energy per site in the future.) -1 # The total number of roots +0.,0.,0.,0.,0.,0.,0.,0.,0. # Energy of each orbital + one extra term +2 # The total number of roots +1 # I The position of the first +1 # I SBox +1 # I +1 # I +1 # II The positions of the second +1 # II SBox +1 # II +1 # II +1 # III +1 # III The positions of the third +1 # III SBox +1 # III +1 # positio of the hole +0 # fix the position of the first hole during the CI +0 # fix the position of the second hole during the CI +0 # Print the wavefunction. It is stored in the FIL666 file after the run ``` 2. running DEHam ```shell -mpiexec -n [nprocs] ./ex1 inpfile +mpiexec -n [nprocs] ./bin/ex1 input_file.inp ``` +_Sample Application_ +-------------------- + +A 2D t-J model Hamiltonian description and setup for using DEHam to solve for few low lying states +is provided in the notbooks folder. Please have a look about the details of using DEHam to study +t-J Hamiltonians. + +![](https://raw.githubusercontent.com/v1j4y/DEHam/master/notebooks/graph.png) + + _Publications using this code_ ------------------------------- diff --git a/examples/eight.inp b/examples/eight.inp new file mode 100644 index 0000000..9c11848 --- /dev/null +++ b/examples/eight.inp @@ -0,0 +1,30 @@ +16 +77 +4 +4 +0 +true +1,2,3,4,5,6,7, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,12,13,14,15 +2,3,4,5,6,7,8,16,15,14,13,12,11,10, 9,10, 11, 12,13,14,15,16 +1,1,1,1,1,1,1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3 +.1430,-0.20,0.0000 +.1430,-0.20,0.0000 +-1.00,0.0,0.00 +0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 +2 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +0 +0 +0 diff --git a/examples/seven.inp b/examples/seven.inp new file mode 100644 index 0000000..8a7cf8d --- /dev/null +++ b/examples/seven.inp @@ -0,0 +1,30 @@ +14 +15 +2 +3 +0 +true +1,2,3,4,5,6, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13 +2,3,4,5,6,7,14,13,12,11,10, 9, 8, 9,10,11,12,13,14 +1,1,1,1,1,1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3 +.1430,-0.20,0.0000 +.1430,-0.20,0.0000 +-1.00,0.0,0.00 +0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 +2 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +0 +0 +0 diff --git a/examples/small.inp b/examples/small.inp index 8e1a91e..7c40809 100644 --- a/examples/small.inp +++ b/examples/small.inp @@ -1,6 +1,6 @@ 8 -140 -1 +35 +4 1 0 true @@ -10,5 +10,21 @@ true .1430,-0.20,0.0000 .1430,-0.20,0.0000 -1.00,0.0,0.00 -1 +0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 2 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +0 +0 +0 diff --git a/examples/tJ_2x2.inp b/examples/tJ_2x2.inp new file mode 100644 index 0000000..d4baba3 --- /dev/null +++ b/examples/tJ_2x2.inp @@ -0,0 +1,30 @@ +4 +1 +1 +1 +0 +false +1,1,2,3 +2,3,4,4 +1,1,1,1 +0.1430,0.0,0.00 +0.1430,0.0,0.00 +-1.000,0.0,0.00 +0.0,0.0,0.0,0.0,0.0 +2 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +0 +0 +0 diff --git a/examples/tJ_3x3.inp b/examples/tJ_3x3.inp new file mode 100644 index 0000000..14ea37c --- /dev/null +++ b/examples/tJ_3x3.inp @@ -0,0 +1,30 @@ +9 +7 +1 +1 +0 +false +1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8,9 +2,3,1,5,6,4,8,9,7,4,5,6,7,8,9,1,2,3 +1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 +0.1430,0.0,0.00 +0.1430,0.0,0.00 +-1.000,0.0,0.00 +0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 +2 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +0 +0 +0 diff --git a/examples/tJ_4x4.inp b/examples/tJ_4x4.inp new file mode 100644 index 0000000..a193afb --- /dev/null +++ b/examples/tJ_4x4.inp @@ -0,0 +1,30 @@ +16 +33 +4 +1 +0 +false +1,2,3,4,5,6,7,8, 9,10,11,12,13,14,15,16, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16 +2,3,4,1,6,7,8,5,10,11,12, 9,14,15,16,13, 5, 6, 7, 8, 9,10,11,12,13,14,15,16, 1, 2, 3, 4 +1,1,1,1,1,1,1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 +0.1430,0.0,0.00 +0.1430,0.0,0.00 +-1.000,0.0,0.00 +0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 +2 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +0 +0 +0 diff --git a/examples/three.inp b/examples/three.inp index 8706cf2..0ac607f 100644 --- a/examples/three.inp +++ b/examples/three.inp @@ -10,5 +10,21 @@ true .1430,-0.20,0.0000 .1430,-0.20,0.0000 -1.00,0.0,0.00 +0.0,0.0,0.0,0.0,0.0,0.0 1 1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +0 +0 +0 diff --git a/examples/twenty.inp b/examples/twenty.inp new file mode 100644 index 0000000..b15a1ff --- /dev/null +++ b/examples/twenty.inp @@ -0,0 +1,30 @@ +20 +77 +4 +3 +0 +true +1,2,3,4,5,6,7,8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,12,13,14,15,16,17,18,19 +2,3,4,5,6,7,8,9,10,20,19,18,17,16,15,14,13,12, 11, 12,13,14,15,16,17,18,19,20 +1,1,1,1,1,1,1,1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3 +.1430,-0.20,0.0000 +.1430,-0.20,0.0000 +-1.00,0.0,0.00 +0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 +2 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +0 +0 +0 diff --git a/notebooks/graph.png b/notebooks/graph.png new file mode 100644 index 0000000..d7e911a Binary files /dev/null and b/notebooks/graph.png differ diff --git a/notebooks/t_J_Model.ipynb b/notebooks/t_J_Model.ipynb new file mode 100644 index 0000000..d10d2f1 --- /dev/null +++ b/notebooks/t_J_Model.ipynb @@ -0,0 +1,501 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 1. t-J Model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "-----" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. t-J Hamiltonian" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### a. Basis" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The $t-J$ Hamiltonian with a S=1/2 spin-model is constructed in a basis of $\\mathcal{R}^3$ given by the vector B as follows:\n", + "\n", + "$$\n", + "B=\n", + " \\begin{bmatrix}\n", + " 0 \\\\\n", + " \\uparrow \\\\\n", + " \\downarrow\n", + " \\end{bmatrix}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### b.) Operators" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Hamiltonian is constructed using two types of operators.\n", + "\n", + "1. The particle creation and annihilation operators $\\hat{C}^\\dagger,\\hat{C}$\n", + "2. The spin ladder operators $\\hat{S}^+, \\hat{S}^-, \\hat{S}^z$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The matrix form of the five operators in the above basis is given as follows:\n", + "\n", + "$$\n", + "\\hat{C}^{\\dagger}_{\\uparrow}= \\left( \\hat{C}_{\\uparrow} \\right)^\\dagger=\n", + " \\begin{bmatrix}\n", + " 0 & 0 & 0 \\\\\n", + " 1 & 0 & 0 \\\\\n", + " 0 & 0 & 0\n", + " \\end{bmatrix}\n", + "$$\n", + "\n", + "$$\n", + "\\hat{C}^{\\dagger}_{\\downarrow}= \\left( \\hat{C}_{\\downarrow} \\right)^\\dagger=\n", + " \\begin{bmatrix}\n", + " 0 & 0 & 0 \\\\\n", + " 0 & 0 & 0 \\\\\n", + " 1 & 0 & 0\n", + " \\end{bmatrix}\n", + "$$\n", + "\n", + "$$\n", + "\\hat{S}^{z}\\ \\ \\ \\ =\\ \\ \\ \\ \\ \n", + " \\begin{bmatrix}\n", + " 0 & 0 & 0 \\\\\n", + " 0 & \\frac{1}{2} & 0 \\\\\n", + " 0 & 0 & -\\frac{1}{2}\n", + " \\end{bmatrix}\n", + "$$\n", + "\n", + "$$\n", + "\\hat{S}^{+}_{\\downarrow}= \\left( \\hat{S}^- \\right)^\\dagger=\n", + " \\begin{bmatrix}\n", + " 0 & 0 & 0 \\\\\n", + " 0 & 0 & 0 \\\\\n", + " 0 & 1 & 0\n", + " \\end{bmatrix}\n", + "$$\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### c.) Hamiltonian " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Nearest Neighbor coupling**\n", + "\n", + "The $t-J$ Hamiltonian is written in the above basis with the above five operators and the two parameters $t$ which governs the kinetic energy of the particles and the magnetic-coupling parameter $J$ which governs the coupling of the localized spins on neighboring sites.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\hat{H_{t-J}} = \\sum^{N-1}_i \\left\\{ t_{i,i+1} \\sum_{\\sigma} \\left(\\hat{C}^{\\dagger}_{\\sigma,i}\\cdot\\hat{C}_{\\sigma,i+1} + h.c. \\right) + J_{i,i+1} \\left( \\frac{\\hat{S}^+_{i}\\cdot\\hat{S}^-_{i+1} + \\hat{S}^-_{i}\\cdot\\hat{S}^+_{i+1}}{2} + \\hat{S}^z_i\\cdot\\hat{S}^z_{i+1} \\right)\\right\\}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$h.c.$ stands for Hermitian conjugate." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 2. 2D Lattice" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "----" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1.) Labeling\n", + "Here a 4x4 lattice is shown with periodic boundary conditions where each site and bond is labelled." + ] + }, + { + "cell_type": "code", + "execution_count": 150, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVwU9f8H8NeyKLAsopSComKaeCCHJ9nhkeKRllmUmlZmh6mVlaX2zbOsTDPFNE1Ns9O0TMsDw7zKr+ZtWWpGoKIoKvd97Of3xycQZJdjf+zM7Hdez8djH7EzQ/v2zey8Zz6fz3zGIIQQICIi0gkXtQMgIiJSEgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgufk1i1apXaITillStXqh2CU2Le7LNixQq1Q6AqMAghhNpBUOUMBgP4p6o+5s0+zJt9mDfnwCs+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFeOMGTNmqB0EVa5+/fro0qWL2mE4nZSUFPTv31/tMJyOEAI9e/ZUOwynw7w5B97OQP/TcnNz4e7urnYYRKQhLHxERKQr7OMjIiJdYeEjIiJdYeEjIiJdYeEjIiJdYeHToMzMTEybNg1BQUHw9vZG/fr1cdttt+GTTz5ROzRNKywsxEcffYR+/fohJCQEISEh6N+/P5YuXYqCggK1w3NKzzzzjNohaFZRURE++ugjTJ06FXv37i2zbtasWSpFRVXBUZ0aNGjQIAwePBi9e/fG2rVrkZWVhaFDh2LWrFnw9/fH22+/rXaImjRs2DDUrVsXjz/+OBo3bgwASEhIwOrVq5GcnIyvv/5a5Qi1KTk52epyIQRCQ0ORkJCgcETO4amnnkJ2dja6dOmCzz77DN27d8f7778PAOjQoQOOHDmicoRkCwufBoWGhuL48eMl7zt37oyDBw/CYrGgbdu2OHXqlIrRaVdgYCD++uuvaq/TO6PRiICAgDKP0yl+vM6FCxeQn5+vYnTaFRISgt9++w2AbG0YO3Ysrl69iq+++gq33XYbjh49qnKEZAubOjXI09MTv/zyCwDg+++/h4+PDwDAxcWFz/qqgI+PD9atWweLxVKyzGKx4Ouvv0a9evVUjEzbmjdvjl27diEuLq7k9c8//yAuLg6+vr5qh6dZpU8IXF1dsWzZMoSGhuLuu+9GZmamipFRZVj4NGjJkiV4+eWXUa9ePcyZMwcLFy4EAFy5cgXjxo1TOTrtWrNmDb755hv4+fkhMDAQgYGB8PPzw/r167FmzRq1w9OsF198ESkpKVbXTZw4UeFonEenTp0QHR1dZtn06dPxxBNPID4+Xp2gqErY1KlRJ0+exIULF3DbbbfBbDaXLI+Ojka/fv1UjEzbfv31VxgMBrRo0QKnTp3Cvn370LZtW9xzzz1qh6ZpBw4cgMFgQOfOnfHnn38iOjoarVu3Zt4qwbw5JxY+DVq4cCE+/PBDtG7dGseOHUNUVBQGDRoEgJ3mFZk5cya2bt2KwsJCRERE4MCBA+jRowdiYmLQt29fvP7662qHqEk35u3XX39Fz549mbdKMG9OTJDmtGvXTmRkZAghhIiLixMdO3YUCxYsEEIIERYWpmZomtauXTtRWFgosrKyhJeXl0hLSxNCCJGdnS2Cg4NVjk67mDf7MG/Oy1XtwkvlWSyWkubNZs2aYdeuXYiMjMTZs2c5uKUCrq6uMBqNMJlMaNGiBerUqQMA8PDwgIsLu7NtYd7sw7w5L/51NMjX1xfHjh0reW82m7Fp0yZcvXoVv//+u4qRaVvt2rWRnZ0NADh8+HDJ8rS0NB6IKsC82Yd5c17s49OghIQEuLq6ws/Pr9y6vXv34o477lAhKu3Ly8uDm5tbueVXr15FYmIigoODVYhK+5g3+zBvzouFj4iIdIXX404iNTVV7RCcUm5urtohOKWMjAy1Q3BKzJtzYOFzEj179lQ7BKfUsGFDtUNwSiEhIWqH4JSYN+fApk4nUTx3IlUP82Yf5s0+zJtz4BUfERHpCgsfERHpCgsfERHpCgsfERHpCgsfERHpCgufRqWkAPPmAS1aAHLazjzUrQv06AFs3gwUFakcoFZdvgy8+SYQEAB4eiIXAHx8gH79gB07AI64syohAXjtNcDfHzCZACAXPj7AAw8A//0v02ZTbCwwfjzQsCHg4YEcALj5ZmD4cIBPUdEs3s6gMenpwNixwLffAi4uwL9TAZZhNgMeHsCsWcAzzygfoyZduSKTsXUrYDAA1m5cN5sBb2/gvfeAoUOVj1GDzp0DnnoK2LNHvs/LK7veYJCF0NcXWLQI6N9f+Rg16fRpmbhDh+RZaEFB2fVGI+DmBjRrBixdCtx1lyphknUsfBqSmCi/HwkJ5Q9A1phMwKhRwMKF8gClW//8A9x5J3D1avkDkDUmEzBpEjBtmuNj07Dff5ctCGlpVWtB8PCQ5wxjxzo8NG3btw/o2xfIzKzapbCHB7ByJU+2NISFTyMyM4GOHeUxvLCw6r/n6Qm8/DLwxhuOi03Trl4FQkOBS5cAi6Xqv2cyAXPmAOPGOS42DTt3DggLk03q1eHhAaxaBQwZ4pi4NO/kSSA8HKju1GQeHsDGjUBEhGPiomphH59GTJ0KnD0LjB4NHDwoW+pWrbq+vk0buTw5Wb5iYuSyrCx5Fq7bpxU9/7xs5hwzxnriSps6VZ6h9+ol25BffVVeXuvQyJGyWX3cOOtpCwiQqcrIuP6aMgXIyQGeeEJeJerSkCHyLNVW4gBZ5BYvlvtlaiqwe7dM3EMPVa0phxyOhU8DcnOBFSvkd+LiRdl3t3Jl2W0uXgQiI+U4jZtvBr7/HlizRq4rKADmz1cm1lGjRqFBgwZo165dybKpU6ciJCQEYWFh6NOnDy5evKhMMMnJwIYNMgG2EleseXN54Ckdm8UCfPihIqFay1uxefPmwWAw4OrVq4rEEh8vW+uKiipPW926gJeXfM2aJZcZDMDq1YqEajVvM2bMgL+/P8LCwhAWFoYtW7YoE8zRo3IwixAVJ27ZMvlFbdNG/vell+TyoiJg/XpFQrW1v33wwQdo3bo1goKCMHHiREVi0STFn/lO5Xz6qRBmsxDyGyVfb74pxKpVZZcVv4xGIcaOFSIr6/oyDw8h0tIcH+vu3bvF4cOHRVBQUMmytFIfHBUVJUaPHu34QIQQYu5cIUymqiVu61Yh+vcXIi5OiF69ri/39hYiP9/hoVrLmxBCnDt3TvTp00c0bdpUXLlyxeFxCCHEK68IUbt2xWkLCJDbGo3W98HGjYWwWBwfq7W8TZ8+XcydO9fxH36jESPKJ+TGxLVqJb+IXl7WExcaqkio1vK2Y8cO0atXL5GbmyuEEOLy5cuKxKJFvOLTgGXLZOtJVaSkyCvEDz4A3n77+nJXV0CJE99u3brBx8enzLI6deqU/JyVlQWDUiNtVqywPuz1RpGR8nJ669by64SQTVEOZi1vAPDSSy9hzpw5yuUM8motP79q2549C5w/Ly9sbrrp+vLUVOD4ccfEV5qtvClOCGDduspHAXXpIpM2c6Zs6vztN3lPSLHTpxVpXreWtyVLlmDy5MklD89t0KCBw+PQKhY+DUhMrPq29erJEfnPPSdbXorl58tb2NTy+uuvo0mTJvjiiy/whlIjbarSNGg2yzOE8eOtrxdCtcRt3LgR/v7+CA0NVfRzq/Jox6tXgU6dZF9fx46yqfOLL66vNxrV3d8WLVqEkJAQjBo1CinVHaFjj5ycqo06a9wYCA6WnaCNGskv6urVQOvWcn3t2kBSkmNjteGvv/7Czz//jPDwcHTv3h0HDx5UJQ4t0GThK+7z0MsrNja+WvnJzpa3Bn36KVC/vlyWl5ePF198tUbjmjt3bpVjeuutt3D+/HkMHz4cixYtshJzNvz8/Go0vtRr1yoPbMYM4LPP5Fm4FRkZGXhixIgajWtNcedrBbKzs/H2229XepJgsVhgNBprNL6CgsrvXcjKAg4flhc4SUny+N23b/FkCkBaWir69bu3RuM6XsVLyDFjxiA2NhbHjh1Dw4YNMWHCBKvb1WRs3p6eyKvKPR85OfIsdNYs2fe8Zw+wcyfQpw8AIDU9Hbd17FijsZ0/f75KeSssLERycjL279+PuXPn4uGHHy73CKWMjIwaja1FixaafDivJgvfhAkTIITQzSskpFm1c+TiIkfk+/vL9yZTbSxfPrdG43r11VerHdfw4cPx7bfflltuMplw6dKlGo2vbvE/viK9egEvvCAvqxMTgSZNgLVrgX879r3q1MGq776r0biGVuF+rdjYWMTFxSE0NBTNmjVDQkICOnTogEuXLpXZzsXFBUVFRTUan9lsrNofsxQhiuOR//X2rouff/6hRuOq6pWvr68vjEYjXFxc8PTTT+PAgQM2Yq652NKKiuDmUoXD5W+/2U4egLpmM/afOlWjsTVp0qRKeWvcuDEeeOABGAwGdOnSBS4uLuUGVHl5edVobLGxsfDy8qpSfErSZOHTm759ZQsIcH3CB6Ox7M+9e8v7rlxcZLPT++/L/r6TJ+XvCQHccYc68Z85c6bk540bN6J1cbOOo919t0wOYDtxvXoB7drJ5IWFydF4o0fL4eaAPDsPD1cm3lKCg4ORlJSE+Ph4xMfHo3Hjxjhy5Aj8/Pwc/tldu17/2VbaunQBAgPlCE4fHzlJws6d8hYIQKZN4RbaEoml+ga+++47qyNla5yLC1D66eq2Erdnj7xJ8rXX5Pvbbwd69gS2bZO/5+oqRxir4P7778fOnTsByGbP/Px83HzzzarEojpBqouPF8LdXQhAiOnTy6+fPl2IyEghTp4UIiNDiKQkITZtEiI4WJQMFuvUSZlYhw4dKvz8/ISrq6vw9/cXK1asEA888IAICgoSwcHBYuDAgSIhIUGZYI4evT6q01bibhxVV3pUp8EgxIABioRqLW+lBQQEKDaqMybm+ihiW2kbOlSIf/4RIjNTiIsXhVi9WghfX1EyqnjUKEVCtZq3ESNGiHbt2ong4GBx7733iosXLyoTzFdfXR+tWdH+1ratEP/9r0zeH38Icf/9crmbmxD/+Y8ioVrLW15enhg+fLgICgoS7du3Fz/99JMisWgRZ27RiF695BzK9jCb5T20kZE1G5NTCA4GTpyw73fNZmDTJqB795qNSeMsFtnia+/tlh4ewIED8kJaV/LzZad68WVvdbm7A2fOyAEwpCo2dWrEnDnFs+JXT61ach7cQYNqPCTnsGCBPBJXl5sb0L490K1bzcekcS4uQFSUfWnz8JATVeuu6AGyP2L2bPu+qCYT8PjjLHoaYZwxY8YMtYMgOfI5KEjOyFLVuTpr15az5u/dC5S6lU5fmjeXSdi+veqJc3MDbrkF2LVL/qxDbdvKArhvX9Xm9QZk0QsOlhfJrq6OjU+zOncGrl2Tg1iqmjiTSc4+/8UX10cHkar4V9CQ+++Xhc9slpNP22IwyG3atZP38un4PlTpqafkvR0eHhWfjbu4yMR27Srb6jQ42kxJr70mL5jd3Su++jMaZVr79pX3+uv0XOG6efPkkz3c3StOhqurTOzQofJswVj90bTkGOzj06DMTODzz2XzZ1KS/P6kpaXB29sbubmyP/DVV2XXlEHPjyO6UWqq7Ox87z3ZD2M0luQNeXnAwIHAK6/IIYtMXIkrV4Dly+V8r3l5MjXp6TJveXlyitOXX5aDYqmUCxeAJUvkCOF/nwySlp4u97eCAuDRR4EXX7x+8zppBgufhgkhp4VKSADuvfchbNu2DsHB8mHPVAEh5N3XiYmIvO8+fBMTI4/aeh26XUVFRfKBA0lJwKBBkdi+/Rt06iRnCqIKFBQAv/4KXL2KBwYPxvodO+TJVUXNNqQqFj4nYTAYwD9V9TFv9mHe7MO8OQf28RERka6w8BERka6w8BERka6w8BERka6w8BERka6w8DmJf/75R+0QnBLzZh/mzT7Hjh1TOwSqAt7OQEREusIrPiIi0hUWPiIi0hUWPiIi0hUWPiIi0hUWPvqfdu3aNbVDICKNYeHToOjo6JKf09LS8OSTTyIkJASPPPIILl++rGJk2jZ58mRcvXoVAHDo0CE0b94c4eHhCAgIwO7du1WOTrs6dOiAWbNmITY2Vu1QnMqhQ4fQs2dPjBgxAufPn0dERAS8vb3RuXNnHD16VO3wqAIsfBr0n//8p+TnCRMmoGHDhvjhhx/QuXNnjB49WsXItG3z5s24+d9HD7366qv4+uuv8ffffyMmJgYTJkxQOTrtSklJQWpqKnr27IkuXbpg/vz5uHjxotphad7YsWMxceJEDBgwALfffjtGjx6NtLQ0zJ49G2PHjlU7PKoAC5/GHTp0CLNmzUJAQABeeuklxMfHqx2SZhUWFqKwsBAAkJOTg86dOwMAAgMDkZeXp2ZomlavXj289957OHfuHObNm4czZ86gQ4cO6NmzJ5YtW6Z2eJpVUFCA/v37Y9iwYTAYDIiMjAQA9OrVC7m5uSpHRxVh4dOgpKQkvP/++5g3bx7S09PLPN/L8u+Tnqm8sWPH4p577sGOHTvQr18/jB8/Hrt378b06dMRxseHV8ldd92FDz/8EBcuXMCkSZOwb98+tUPSLHd3d/z4449Yt24dDAYDNmzYAADYvXs3jEajytFRRVzVDoDKe/rpp5GRkQEAePzxx3H16lXUr18fly5d4gG8As8//zyCg4OxZMkS/PXXXygsLMSZM2dw//33Y8qUKWqHp1mBgYHllhmNRvTr1w/9+vVTISLnsHTpUkycOBEuLi7Ytm0blixZgpEjR8Lf3x/Lly9XOzyqAKcs06hTp07hwoULCA8Ph9lsLlkeHR3Ng1EFmDf7MG/2OXnyJC5evMi8ORtBmrNw4UIRGBgoBg0aJAICAsSGDRtK1rVv317FyLQtKiqKebMD9zf7REVFiVatWjFvTohNnRq0bNkyHD58GGazGfHx8YiMjER8fDzGjx9fpr+Pylq+fDnzZgfub/ZZvnw5Dh06xLw5IRY+DbJYLCXNJs2aNcOuXbsQGRmJs2fP8gtVAebNPsybfZg358VRnRrk6+tb5rleZrMZmzZtwtWrV/H777+rGJm2MW/2Yd7sw7w5Lw5u0aCEhAS4urrCz8+v3Lq9e/fijjvuUCEq7WPe7MO82Yd5c14sfEREpCts6nQS2dnZaofglJg3+zBv9rly5YraIVAV8IrPSRgMBnaY24F5sw/zZh/mzTnwio+IiHSFhY+IiHSFhY+IiHSFhY+IiHSFhY+IiHSFhU/rrlwBjh9HOACcPAnwAZdVcukSICfV6ILTp4H8fLUjchIXLgBHj6IzAJw5A/z7YF+yTQjg3Dng8GEA6ITYWICPzdQ2Fj4tKioCtmwBevQAmjQBunVDNACEhwM33ww8/7w8KFEZhYXA+vUyTc2aAd27A8A2dOoE1K8PTJokD1B0g7w84MsvgdBQ4NZbgR498CMAtG8P+PoCM2bIMwkqIzsb+PhjoFUroHVr4O67ASAGISGAvz8wZw5w7ZraUZI1vI9Pa44fB+65B0hPBzIzrW9TqxZgNAL9+skDloeHsjFq0N69wP33y2P4v8/wLad2bcDFBRgyBFi+XKZR97ZuBYYOlZctthLn7i7XjxkDzJsnk6hza9YATz0FGAy2v6YeHjJt//kPMGWK3Ja0gYVPS/buBfr2BbKyqra9u7s83fzlF6DUQzD1ZvNm4OGH5Rl4VZhMQKdOQEyMLIa69fnnwDPPADk5VdveZAIiIoBvv5UnXjq1YIEsZlVNm6enPLdYvpzFTytY+LTizBmgY0fbZ922uLsDXbsC27fr8kz8yBHgrruqXvSKeXgAAwcCa9c6Ji7N27kTGDCg6kfvYiYT8MQTwKJFjolL49avB0aMqH7aPD1lU/vUqY6Ji6pHf0dKrZo0yfqVXuvWwE8/Aampsjjef3/Z9bm5wMGD8vJFBVFRUWjXrh2CgoKwYMECxT//uedk0Rs3TqYhNxdYtarsNnffLccFZWUBO3YATZvKA9fmzcUDEhxv1KhRaNCgAdq1a1eyLDk5GREREWjZsiUiIiKQkpKiTDBCXL/Ss5W4WrWAdeuAuDi5vewwvd6xFRenSKjW8rZu3ToEBQXBxcUFhw4dUiQOQHa9P/tsxWkLDwd+/FH27SUlyRMrPz+57739tnJ9ftby9uqrr6J169YICQnB4MGDkZqaqkwwGsTCpwVJSbKv5cahYEYjsHEjsGkT4OMjD1affw60bFl2u8xM4N13lYv3XydOnMDy5ctx4MABHD9+HJs2bcLff/+t2Of/9Rdw9Kj8+eJFYNYsYOXKstvcdJM8S586Vabw0CHg66/lutxc2WWlhJEjRyI6OrrMstmzZ6NXr144c+YMevXqhdmzZysTzK+/AomJ8mdbiQNkE/qIEde3LWaxAAsXOj5OWM9bu3btsH79enTr1k2RGIpt3Xp9ULWttNWrByxbJgdXBQTIBpziwmgwACtWKBOrtbxFRETgxIkT+O233xAYGIh33nlHmWC0SJD63nxTCHd3IeS59fVXUJAQGRlll23bJsQbb5Tf1t1diLg4RcNeu3atGDVqVMn7N954Q7z77ruKff64cULUqlU2DW++KcSqVdffP/20EHv3Xn9vMgmRnS1Eq1bX05acrEy8cXFxIigoqOR9YGCguHjxohBCiIsXL4rAwEBlAnnwQSFcXCpOXOnX+fNCdO9edpnZLEROjiLh3pi3Yt27dxcHDx5UJAYhhLjrrvKpqShtgBDt2wuRnn79fYMGQhQVKROvrbwJIcT69evFI488okwgGsQrPi1Yv77q9+cZDECp5osSRqNsx1NQu3bt8PPPP+PatWvIzs7Gli1bcP78ecU+/4cfgIKCircJCpIDZYtlZwOxsXI5IAe3/Pyz42KsyOXLl9GwYUMAgJ+fHy5fvqzMB8fE/P9vNHNxkR2sOiGEHHtWXd26AX/8cf19Vpbc/9S2cuVK9O/fX+0wVMPCpwW2+nZOn5bNoK++Cri6yhF13bvLAQY3ys+3/f9xkDZt2mDSpEno06cP+vXrh7CwMBgVHO2Xnl75NmYzkJZWdllaGuDlJX+2WBRPm1UGgwEGpYb8VXXUcEUMBm0kTiEZGdUfyBocDEybJr++xVxd1U/bW2+9BVdXVwwfPlzdQFSkycI3e/bskgOBHl6x8fHWE1FYKAezDBggbyCeMEH2licklNs0r6AAL73yisNifPbZZ60+nPTJJ5/E4cOHsWfPHtSrVw+BgYEl6woKChyat9TU5Er3pcxMoE6dssvq1Lk+eDYzMx0jRz5ao3GtWbOm0rgAwNfXF4n/9p8lJiaiQYMG5baxWCwwGo01Gl9BUVGV4qtIaloa+g8cWKNxHS99aV4DajI2b28zCiprXiilRQvZJzh+vOwqLZaWlobw8E41Glt1Wlk++eQTbNq0CV988QUMhvInWhkZGTUaW7NmzZBp60ZHFWmy8E2ePBlCCN28WnTpYjsZv/8uZ3C5+WZ5w3rz5sCBA+U2czObMf/TTx0W49KlS2GycqWZlJQEADh37hzWr1+PRx55pGRdrVq1HJq35s19Kt2X/vhDTkhSzGSSB6Xi5qc6depgy5bPajSuoUOHVhoXANx3331YvXo1AGD16tUYNGhQuW1cXFxQVFRUo/HVqlevSvFVpK6XF7YeOFCjcYWW/kPVgJqMzWLJhKtr1WY8aNpU3l305ptyLFppHh7eiIs7VKOxNWnSpEpxRUdHY86cOfj++++tfpcBwMvLq0Zji4+Ph1mD9xhrsvDpzmOPyRt9rAkOBtzc5I1nEyYADRsCn3xSfruCAlkYFfbggw+ibdu2uPfee7F48WLUrVtXsc9+5BF5GyMgm6Hc3OR/S//83XeyS/SBB+SyadOA336TrciAHKJePFLfkYYNG4auXbvi9OnTaNy4MT7++GNMnjwZMTExaNmyJbZv347Jkyc7PhAAiIyUbW6A7cQBsgPUza38z4D8uUMHh4dqLW/fffcdGjdujH379mHAgAHo27evw+MwGOSESsUXSbbS1qiR7GpftAj46KPy/5+GDeVoT0ezlrfnnnsOGRkZiIiIQFhYGJ599lnHB6JVgtSXliaEh4ewOixszhw57DAjQ4gtW4Ro0aL8Ni4uQjzwgNr/CsVduCCEm5tMwfTp5ddPny7X9eolxMmTcjTnzp1CBATI5bVqyZGhuvPHH9f3t4oSZ22UcECA/N033nBsjBq0Z48Qnp4Vp23aNPlzRkbZFyB/d+lSJSMmWzhzi1aMGgV89pl9s+GbTLJtpWvXmo9L4+65B9i2zb5Biu7ucsRnqW5J/ejUyf67993dgfh4OYG1jgghexpsdclXxmQCLl/W9eyCmsGmTq2YPVvebW2lw7lCJpOcqPK22xwTl8YtXlx+8EpVmEzACy/otOgB8q5qW83rFTGZ5BQkOit6gPxqfvGFfXPCm0zA0qUselrBwqcVDRoAe/bI4lfc/1IZkwno1UvXs9/ecou82PX2rvpUpSaTfEKDUhOlaFJwsJwRqDrFz2SSc3W99JLj4tK4228HvvqqesXPZAJmzgQefdRxcVH1sPBpSWCgfHpqeLj8Ztl6bo6np1z/wgvAhg1VL5T/ozp2lFORhYbKtNi638pslqmbNk1ON6nTc4XrevSQd2W3aiUTY+vMwWyWZxbvvy8fMqdzgwbJOQACAmRqbO1HXl7yPHbVKuCVV5SNkSqhdicj2XDqlBBjxsgecaNR5ABCGAxCNGsmxIcfynmQqJzjx4UYOVIOenF1FQLIEQaDEK1by6mlsrPVjlCjDhwQ4uGHhahdWwhXV5FbvL+FhQmxZo0QeXlqR6g5FosQu3cLMXCg3Nfk9Hm5wsVFiK5dhdiwQYiCArWjJGs4uMUZ5OTgJpMJ1woKdH91V1VCyFn0PT19UFiYrOfHxyHxVtwAAB50SURBVFWPEEBWFup5eSGlsFDXz92rDiGKJ0uoi6KiVD0+IcypsPA5CYPBAP6pqo95sw/zZh/mzTnwvISIiHSFhY+IiHSFhY+IiHSFhY+IiHSFhc9J+PhU/iQCKo95s48LhyXahXlzDhzV6SRycnLgYc9cSTrHvNnHYrHwIG4H5s05sPAREZGu8NSEiIh0hYWPiIh0hYWPiIh0hYWPiIh0hYWPiKiGJCUlqR0CVQELn5Pp37+/2iFoVnp6Ol577TU8+uij+PLLL8usGzt2rEpRad+lS5cwZswYjBs3DteuXcOMGTMQHByMhx9+GImJiWqHp1nJycllXteuXUOXLl2QkpKC5ORktcOjCvB2Bg06cuSI1eVCCAwcOJAHIxsefPBBtGzZErfddhtWrlyJWrVq4csvv4Sbmxs6dOhgM696169fPwwYMABZWVn48ssvMXz4cDzyyCPYsGEDtm/fjo0bN6odoia5uLggICCgzLKEhAQ0btwYBoMB//zzj0qRUWVY+DTIaDSie/fuVh9vsn//fuTk5KgQlfaFhYXh2LFjJe/feustbNmyBd9//z0iIiJY+Gxo3749jh49CgBo2rQpzp07V7LuxpzSdfPmzUNMTAzmzp2L4OBgAMAtt9yCuLg4lSOjyvCpphrUpk0bfPTRR2jZsmW5dU2aNFEhIueQl5dXZuaM119/Hf7+/ujWrRsyMzNVjk67LBZLyc+PPfZYmXVFRUVKh+M0JkyYgCFDhuCll15CkyZNMHPmTBgMBrXDoipgH58GzZgxo8zBqLQPPvhA4Wicx7333osdO3aUWTZy5EjMmzcPtWvXVikq7Rs0aFDJicGsWbNKlv/9999o1aqVWmE5hcaNG2PdunXo0aMHIiIikJ2drXZIVAVs6tSoU6dO4cKFCwgPD4fZbC5ZHh0djX79+qkYmbbZytvWrVs5MKgC3N/sUzpvRqMRsbGxaNeuHfOmdYI0JyoqSgQGBopBgwaJgIAAsWHDhpJ17du3VzEybVu4cCHzZgfmzT78njov9vFp0PLly3H48GGYzWbEx8cjMjIS8fHxGD9+vNUBLyQtW7aMebMD82Yffk+dFwufBlkslpLmpmbNmmHXrl2IjIzE2bNn+YWqAPNmH+bNPsyb8+LgFg3y9fUtM4TcbDZj06ZNuHr1Kn7//XcVI9M25s0+zJt9mDfnxcEtGpSQkABXV1f4+fmVW7d3717ccccdKkSlfcybfZg3+zBvzouFj4iIdIVNnU4iJiZG7RCcEvNmnzVr1qgdglNi3pwDr/ichMFgYIe5HZg3+zBv9mHenAOv+IiISFdY+IiISFdY+IiISFdY+IiISFdY+IiISFc4ZZmG5eQAmzcD588DwGtYvBjo2BEIDwf42K8KZGQAmzYBiYmYDABLlgBduwJhYWpHpmkpKTJtSUkAMAkffQR06wa0aaN2ZBqXlARs2QJcuYKJALB8OdCrF9C8udqRkQ28nUGDYmOBhQuBlStlgcvLA/LzC+Hu7gqjEfDzAyZOBIYPBzw91Y5WQ/78E3j/feDLLwFXVyAvDwX5+ajl4SET2by5TNxDDwHu7mpHqxlHjgDz5gHr15ekDQUFBTCZakEIICgImDQJGDQIqFVL7Wg1Qgjgv/8F3nsPiI4GjEYgPx8FBQWoZTIBFgvQubPc3/r3l+tJM1j4NGb5cmD8eKCwECgosL2dpydQpw6waxcQGKhYeNokBDB7NvDmm0B+PlDRU8PNZsDXVyaucWPFQtQii0Uel5csAXJz5XtbzGZ53vDTT8DNNysXoyYVFgJPPQWsWyebZSo6hJrNQPv2sunGy0u5GKlCLHwaMn8+MGUKUNWHOBsMsvgdOKDz4jd5MvDBB1VPnNEI+PgAR48C/v6OjU2jhACefBL4+uuqp612baBhQ3mF6OPj2Pg0y2IB7r9fngFUNXFubsCttwK//somGo3g4BaN+PFH4PXXq/5dAuTBKz0d6NEDyMpyWGja9tVX1St6gLwiTEkB7r5bnr3rUFRU9YoeIC+mExOBe+6p+CLnf9rrr1ev6AGy7Tg2FhgyxHFxUbWw8GnEa6/JVpMb7dwpl2dkyNepU2XXFxc/taYIbNasGYKDgxEWFoZOnTop++FCyKs9WwehIUNkv19mJvD338Cdd15fV1gIXLwoByWoIDU1FZGRkWjdujXatGmDffv2KfbZBQXAzJnl01a8jxW/CgtlX3Np+fnAiROAguGWOH36NMLCwkpederUwYIFC5QLIDNTnjFY298CAmRzZnKyPDv44IOy/Xq5ubJg3vgFVsj8+fMRFBSEdu3aYdiwYcjNzVUlDs1Q7mHvZMuJE0J4eAghj+RlXzt3CvHkk9bXlX7deqsQFovysQcEBIgrV64o/8FCCLFrlxBms/WE9O4tRHy8EOHhQhgMQjRqJF83bnf77aqE/thjj4nly5cLIYTIy8sTKSkpin32N9/YTlvxy9NTiIwMIe66q/w6g0GIwYMVC9eqwsJC4evrK+Lj45X70CVLZGKsJWzzZiFWrRLCzU0IX18hfvtNiOefL7uNq6sQzz6rXLz/SkhIEM2aNRPZ2dlCCCEeeughsWrVKsXj0BJe8WnAggXyTPr/IzEROHiwZuJxGu+9Z7uNd+ZM4I03ZL+KEPLq7uLF8tsdOSKboRSUlpaGPXv24MknnwQA1K5dG3Xr1lXs8+fMkRcvFXnwQTlK/+efy68TAti6FbhyxTHxVcVPP/2EFi1aICAgQLkPnTvX9v52yy3A2rWyWfPyZTnSMyio7DaFhcCnn1avmbSGFBYWIicnB4WFhcjOzkajRo0Uj0FLWPg0YN++igcivvOOPMj88gvQvbv1bSwWOVZDaQaDAX369EHHjh2xbNkyZT/88GHrnU0uLkCnTkD9+sCZM/JGyA8+sH4LQ+3aQKmnaCshLi4O9evXxxNPPIH27dvjqaeeQpaCnbQnTlS+zeOPy2O0LW5ushVZLWvWrMGwYcOU+8CiIiAuzvb6BQuAoUMBDw+gUSN5C0N0dPntjEYgPt5hYVrj7++PV155BU2bNkXDhg3h7e2NPn36KBqD1rDwaUB6uu11kybJYeT+/sCyZcAPP1i/L7agAEhLc1yMtvzyyy84cuQItm7disWLF2PPnj3KfbitYuHrKwtaZCRw113yxvX27eWQ2RsVFSmeuMLCQhw5cgRjxozB0aNH4enpidmzZyv2+db6kktr2lSeYK1eXfF2auxvAJCfn4/vv/8eDz30kHIfmp5e8U2Me/bIK7z0dODCBeDQIWDDhvLbubgonriUlBRs3LgRcXFxuHjxIrKysvD5558rGoPmqN3WeiOLxSJmzpwpAOjodarSPrzi19atQjz3nLV12QIY57AYR48eLbKysir8202fPl3MnTu35H1+fr5D83bZVpLq1pUBPPbY9WUPPCDEkSPltk0DxBAHxjhlyhRRVFRUJk+JiYkiICCg5P2ePXvEPffcU/L+2rVrDt7fcircx15/XXafVrwvpgigl8NiXLJkic39bMOGDSIiIqLc8tjYWIfF4waIQlvJMBhkf/J//iNE7dpC+PgIsWGDEO++W27bVECEOPBv+/HHH5fLy9q1a8WoUaNK3q9evVqMGTOm5H10dLTD4mnUqJHIzMy0fdBQieYKnx5FRFSt6AFCbNlSvs8cEMLLS4j165WNOzMzU6Snp5f83LVrV7F161blAggNtZ2oc+eEePTR6+8HD7Za+ISXlxA//6xczP+68847xalTp4QQ8oThlVdeUeyz/f0r3sdOnxbiiScq3sZkEuLkScVCLmPIkCFi5cqVyn+wrRFBN90k19epc33ZoEFC/P57+W3d3IS4fFnRsPfv3y/atm0rsrKyhMViEY899phYuHChojFoDQufBmzcKI+/N35HvL2F6NNHfleMRiEeeUSIzEwhWra0fvzOzVU27tjYWBESEiJCQkJE27ZtxaxZs5QNYMUK2wejmTOFOHBAiPr15RXgnj1CvPFG+e38/IS44YpMCUePHhUdO3YUwcHBYtCgQSI5OVmxz37rLdujiLt2lftYZaM+27ZVLNwyMjMzhY+Pj0hNTVX+w194QYhatawnJDZWiEmT5BfV21uehX7xRfntundXPm4hxLRp00SrVq1EUFCQGDFihMhV+mChMSx8GlBYKE8ab/yO3HyzPHanpwuRkiLEvn1ylP6N29WuLcSECWr/K1SQlSUvPawdiFxdhVi8WCYuMVGIqCh5BlF6Gw8PIUo1zepFUpIQ7u7W07Z0qRCfflpx0TObhVi9Wu1/hQr+/tt24kJD5b1HyclCXLkixNdfC9GgQfmz082b1f5XkBCCU5ZpxNtvA2+9Zd9IZw8POcKuWbMaD0v7nn8eWLFC3iBcXZ6ewLlzupx/6+GHgY0b7buNxtsbuHRJp/N89+gB7N1r34w/jRrJ/Y0TVquOozo1YuJEOQK/ugcTkwn48EOdFj1A3pTWsqUcxVkdJpOc7kaHRQ+QI4QbNZJPY6gOk0k+ukiXRQ+4vs+4VPPQ6ekpZwli0dMEFj6NcHWVMx6Fh8sruKrw8ADefRcYOdKhoWmbh4ec161Nm6ofjT085FXiwIGOjU3D6taVFy4BAfKevMoYDPLY/c03ZWd+0x0/P5k4P7+qPaPJxUXOJL9tGxAa6vj4qEpY+DTEbAa2bwemTpWPfrH2FBNXV3l879JFnnk/95zycWrOTTcB+/cDEybII7q1xNWuLRPXvTuwYweg5M3PGtWokZy4ZswYue+ZzeW3cXOTr3795OPn+vdXPk7NufVWOenByJHyJMraExc8PGTi7r9f3tN3xx2Kh0m2sY9Po4qK5LRQixcDZ88CJ0/+g5CQ5rj9dvm8vtat1Y5QowoK5I3DS5cCFy7gn9On0TwsTD6J4bnn5NRSVE5urny83IoVsv/ur7/+Qfv2zXHPPbIw6vTpTZXLzAS++EJOc3PlCv4+cwa3duoEDB4MPP20nD2INIeFz0kYDAbwT1V9zJt9mDf7MG/OgU2dRESkKyx8RESkKyx8RESkKyx8RESkKyx8RESkKyx8TiIyMlLtEJxS79691Q7BKXF/sw/3N+fA2xmcRFpaGry9vdUOw+mkpqaibt26aofhdLi/2Yf7m3Ng4SMiIl1hUycREekKCx8REekKCx8REekKCx8REekKCx8RkR3S0tIwefJktG7dGj4+PrjpppvQpk0bTJ48GampqWqHRxVg4dOgDh06YNasWYiNjVU7FKdy6NAh9OzZEyNGjMD58+cREREBb29vdO7cGUePHlU7PM3KzMzEtGnTEBQUBG9vb9SvXx+33XYbPvnkE7VD07SHH34Y9erVw65du5CcnIxr165h586dqFevHh5++GG1w6MK8HYGDbrlllvw4IMPYu3atfDz88OwYcMwZMgQNGrUSO3QNK1Lly6YOXMmUlNTMXHiRMyfPx+RkZH46aefMGXKFOzbt0/tEDVp0KBBGDx4MHr37o21a9ciKysLQ4cOxaxZs+Dv74+3335b7RA1qVWrVjh9+nS115H6WPg0qEOHDjhy5AgA4Oeff8ZXX32F9evXo02bNhg2bBieeeYZlSPUpvbt25dc2TVt2hTnzp2zuo7KCg0NxfHjx0ved+7cGQcPHoTFYkHbtm1x6tQpFaPTrj59+qB37954/PHH4evrCwC4fPkyPvnkE8TExGD79u0qR0i2sKlT4+666y58+OGHuHDhAiZNmsSrlgq4u7vjxx9/xLp162AwGLBhwwYAwO7du2E0GlWOTrs8PT3xyy+/AAC+//57+Pj4AABcXFz4UNUKfP3117h27Rq6d++OevXqwcfHBz169EBycjLWrl2rdnhUEUGaM2TIELVDcErHjh0Tffr0Ef369RMnT54UL7zwgvD29hZt27YVe/fuVTs8zTp+/Ljo3LmzqFu3rrjjjjvE6dOnhRBCJCUliaioKJWj07aTJ0+KmJgYkZGRUWb51q1bVYqIqoJNnU5m1apVeOKJJ9QOw+kwb/Zh3mxbuHAhFi9ejDZt2uDYsWOIiorCoEGDAJTtriDtYeFzMjf2XVHVMG/2Yd5sCw4Oxr59+2A2mxEfH4/IyEg8+uijGD9+PPuUNc5V7QCovJCQEKvLhRC4fPmywtE4D+bNPsybfSwWC8xmMwCgWbNm2LVrFyIjI3H27Fn2jWocC58GXb58Gdu2bUO9evXKLBdC4Pbbb1cpKu1j3uzDvNnH19cXx44dQ1hYGADAbDZj06ZNGDVqFH7//XeVo6OKsPBp0MCBA5GZmVnyhSqtR48eygfkJJg3+zBv9vn000/h6lr2EOrq6opPP/0Uo0ePVikqqgr28RERka7wPj4nkZOTo3YITik9PV3tEJwS9zf7cH9zDix8TuKVV15ROwSnNGbMGLVDcErc3+zD/c05sKnTSRgMBo4UswPzZh/mzT7Mm3PgFR8REekKCx8REekKCx8REekKCx8REekKCx8REekKC59WCQHs3w888QTQrRsOAEBEBPDGG8ClS2pHp1lCADt3AiNGAHfdBQC/om9fYO5c4No1taPTMIsF2LoVGDIEuPNO7AeA/v2BRYsA3ptmU2EhsH498OCDgJzdbR8GDgQ+/hjIzlY7OrKFtzNojRDAJ58Ab78NJCYCOTnyoFTM3V1uExEBzJwJdOigWqhaUlgILFkCzJkDpKYCWVkyTcVMJpnG++4DZswA2rRRLVRtyc0FoqKA99+XR+rMzLLrPT1l4oYMAaZNA265RZ04NSYjQ55MLVok972MjLLrzWaZtpEjgddfBxo1UiVMsoGFT0vy8oBhw4Bt2yo/XTQYAA8PYOVKeVDSscxM4N57gQMHKk+bi4tM27ffAn37KhOfZiUnyxOokyflCVZFjEZZBLduLb600a2LF4EePYDz5+V5Q0Vq1QLq1AF27ABsPASDVMDCpxUWC/DAA8CPP1Z+ECrNZAK+/BL49wGYelNQANx9N3DoUOUHodJMJnkM79bNcbFpWnY2EB4O/PUXkJ9f9d/z9AR++QWwMqG1HqSkAO3bAxcuyCu9qvL2lvvorbc6LjaqOvbxacWiRUBMTPmiN24ccPCgPKqvWlX+97KzgUcekaehOjR9OnDkSNmiV7s2sGIFEB8vu6eOHgX69Sv7e9nZ8ipRt91Xzz0H/P13+aL32WdyX0pLA06fBp58suz6rCygT5/qFcv/IY8+KnsgbBW9W2+VX+HPPiu7PCNDtjDwMkMbWPi0wGIBZs+23k538SIwa5Zs0qzo95cscVx8NowaNQoNGjRAu3btFP9sQLYML1pUPm2urrIZqnt3eaY9ZQqwdi0QEFB2u6Ki8gcoJZw/fx49e/ZE27ZtERQUhKioKGUDSE0FvvrK+iXyO+8AzZrJxN13n9z3buxHzs0FNmxQJNSyH5uLLl26IDQ0FEFBQZg+fbqin3/+PPDTTxXX/MWL5XnqjSwWIClJDrxSQ1FREdq3b4+BAweqE4DGsPBpwfbt5XvHi333HbBxY8VDEnNzZQUoKHBMfDaMHDkS0dHRin5mad98Y315drYc93P2rDzD3rwZiIsDOnYsu11WlhygoPRZuKurK+bNm4c///wT+/fvx+LFi/Hnn38qF8Ann8jOTmv+/PP6kV0I+WrRouw2GRnAu+86NERr3NzcsGPHDhw/fhzHjh1DdHQ09u/fr9jnL15c8b4yZIg8p/jpJ+vrMzPl/qaGqKgotOGIrhIsfFqwcGH50XTVVVQk+wcV1K1bN/j4+Cj6maUtWGD7fKG0Bg2AwEDgjz/Kr7t2TQ6KUVLDhg3R4d+rKC8vL7Rp0wYXLlxQLoCFCyseBbR4sTwrOH1atutt2VJ+m5MngdhYx8VohcFggNlsBgAUFBSgoKAABoNBsc//6CPZymCNl5e80+jllyv+f+zcqfxtNQkJCdi8eTOeeuopZT9Yw1j4tODvv////4+CAtmppSPnzlW+jasr8MUXwOrV8jh+I4NB3bTFx8fj6NGjCA8PV+5DK+sPHjdOHsnvvFPepGbtaF+7tiqJKyoqQlhYGBo0aICIiAjF8pafX3F/8Jtvynv3Kjt/cXOrfJua9uKLL2LOnDlwsXWVr0OazMT8+fNhMBh08zpn7YhcTQW5uZj43HMOi/HZZ59FdjXvyC0+I3fUKykprcLPNxhkH15+vhzLYU1GRiaGDh3lsBinTp0KS+n7MEvJzMzEgw8+iAULFqBOnToly5OTkx2aN4uty5bSLBZg716gcWPAyjPmUtPScF/v3g6LcebMmVbDMhqNOHbsGBISEnDgwAGcOHGiZN3Ro0cdFo+bWz1YLNbzFhoK9O4NzJ9feVrT01MRGtrVYXGuvGEswKZNm9CgQQN0vLGdv5Rt27Y5LJ6WLVsioyrNMkoTpL62bYt7U2y/3nxTiFWrbK/38BBiyRLFQ4+LixNBQUGKf64QQvj5VZyylSuF2LFDCHd329vUqSPEunXKx56fny/69Okj5s2bp/yHe3hUvr8Vv5YvF2LBgvLLvb2F2LlT+dhLmTlzppg7d64in1VQIITBYD1F48cLkZkpRGKifGVkCJGdLcThw+W39fIS4sQJRUIWQggxefJk4e/vLwICAoSvr6/w8PAQw4cPVy4AjdLkFZ/udOhge7CB0SjbR4zGsj9b265tW8fGqTEVDSZdskTOznLvvRXf31dYCAQF1XxsFRFC4Mknn0SbNm3wcmWdQo7QsqX15fXryxEanp5yf+zTR06oYG20Rl4e0KqVY+O8wZUrV5CamgoAyMnJQUxMDFq3bq3IZ7u6Av7+1tctWybH/4SFydfSpXJAlbUJEgoKyo8udqR33nkHCQkJiI+Px5o1a3D33Xfj888/Vy4AjWLh04IXX5RTkVkzZYo8cr/2mryJKDdXLrtRvXrFk1MqZtiwYejatStOnz6Nxo0b4+OPP1b08ydMkFND3ahpU+DZZ+VB6NIlOQAmI0Pe7nijVq2Un75s7969+Oyzz7Bjxw6EhYUhLCwMW6wNIHGUV1+1njghZLNmQoK8U/u99+S++cMP5bft3h1o2NDxsZaSmJiInj17IiQkBJ07d0ZERISiw/NfeklOfHCjnBzg8uXrr8xM+TW9erXsdkYj8NBD1lNPyuLMLVrRtq0cKWcPT095/9Xzz9dsTBpnscg5EC9ftu/3vbyA5ct1OONbbq4c6mpv34vZLO/j69WrZuPSuJQUub9VZ4ag0jw8gP/+V7eT3mgKr/i04q23rJ9OVoWbG/DYYzUbjxNwcZEzt9iTNhcXOYfi4ME1H5fmubvLy2V7EufqKi+p77675uPSuHr15MNS7EmbmxvQqROLnlaw8GnF4MGyWam63yqzWd4A7+3tmLg07tlnZRNmddJmMMirvV275Kh8XZo6VU5Q7eFR9d8xGuXRPyZGJlGHoqJkl7ytnglratUC/PzkPBSkDSx8WjJrFjB5sjyKV3bPjbu7PAjt3i1nzdUpg0HeWPz00zJtlR2PTSbZyvfrrzqfMNjFBVi3Tk6M7ulZ+fYmE9CkCXD4sK6fsVOrlpwnomfPqqXN01P2Ix88KL+upA0sfFpiMMgz8V275JMt3d3LXsoUX6rcfLN8yNfp03weH+QxfMECIDoaGDBANiuVPiM3GOSFccOGcnaNU6cUH5CoTbVqyRsd16+XTZfu7jJ5xVxc5JG7WTM50OXECVn8dM7DA9i0Cfj8c6BrV5m20i0HxU9wat1aziR48KAcMEvawcEtWnbtmpxd+dw5fDB7Np6fPl0WugEDrN/SQADkLFvr1skZMubMWYg33ngB4eHyJmNOXlGBs2dlEUxMxMK5c/HCm2/K0Zt33qnbps2qOHNGjvW5fBmYN28+3nrrJUREAJ07qx0Z2cLC5yQMBgP4p6o+5s0+zJt9mDfnwPNfIiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+IiLSFRY+J1FUVKR2CE6JebNPYWGh2iE4JebNORiEEELtIIiIiJTCKz4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItIVFj4iItKV/wOWKFdG5mwvRQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib as mplt\n", + "import networkx as nx\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "Lx = 4\n", + "Ly = 4\n", + "#fig = plt.figure()\n", + "#draw2DLattice(Lx,Ly,fig)\n", + "draw2DLattice(Lx,Ly)\n", + "#fig.savefig(\"/tmp/graph.png\",dpi=1000)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2.) Generating input file \n", + "\n", + "In this work, we will be using the [DEHam](https://github.com/v1j4y/DEHam) code to solve the Hamiltonian using Exact Diagonalization method. \n", + "\n", + "The input file for the code is as follows:\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + " 1. The DEHam program requires an input file which \n", + " has the topology of the Hamiltonian and the various parameters\n", + " as explained below in a sample inputfile:\n", + "\n", + "```python\n", + "8 \t\t\t\t\t\t\t\t# The number of orbitals (total)\n", + "1 \t\t\t\t\t\t\t\t# The largest number of non-zero elements per row (Multiple of Ndet)\n", + "1 \t\t\t\t\t\t\t\t# The total number of processors used in parallel (Multiple of Ndet)\n", + "1 \t\t\t\t\t\t\t\t# The number of holes\n", + "0 \t\t\t\t\t\t\t\t# The isz (ms-1/2) value\n", + "true \t# Restrict the hole to the 1'st (i.e. half of natom) Family of states. *false* for no restrictions\n", + "1,2,3,1,2,3,4,5,6,7\t\t\t\t\t# (L1) The topology of the system is specified here\n", + "2,3,4,8,7,6,5,6,7,8\t\t\t\t\t# (L2) first and second line contain the two sites linked\n", + "1,1,1,2,2,2,2,3,3,3\t\t\t\t\t# (ltype) third line contains the type of link (1 for t or J, 2 for K and 3 for none)\n", + ".1430,-0.20,0.0000\t\t\t\t\t# The three types of links this line gives Jz, Kz\n", + ".1430,-0.20,0.0000\t\t\t\t\t# The three types of links this line gives Jx,y; Kx,y\n", + "-1.00,0.0,0.00\t\t\t\t\t\t# This line gives t(i,i+1), t(i,i+2), t(i,i+3)\n", + "0.,0.,0.,0.,0.,0.,0.,0.,0. \t\t # Energy of each orbital + additional energy\n", + "2 \t\t\t\t# The total number of roots\n", + "1 \t\t\t\t# I The position of the first\n", + "1 \t\t\t\t# I SBox\n", + "1 \t\t\t\t# I\n", + "1 \t\t\t\t# I\n", + "1 \t\t\t\t# II The positions of the second\n", + "1 \t\t\t\t# II SBox\n", + "1 \t\t\t\t# II\n", + "1 \t\t\t\t# II\n", + "1 \t\t\t\t# III\n", + "1 \t\t\t\t# III The positions of the third\n", + "1 \t\t\t\t# III SBox\n", + "1 \t\t\t\t# III\n", + "1 \t\t\t\t# positio of the hole\n", + "0 \t\t\t\t# fix the position of the first hole during the CI\n", + "0 \t\t\t\t# fix the position of the second hole during the CI\n", + "0 \t\t\t\t# Print the wavefunction. It is stored in the FIL666 file after the run\n", + "```\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### a.) Script for $t-J$ input file" + ] + }, + { + "cell_type": "code", + "execution_count": 139, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]\n", + "[2, 3, 4, 1, 6, 7, 8, 5, 10, 11, 12, 9, 14, 15, 16, 13, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1, 2, 3, 4]\n", + "[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n" + ] + } + ], + "source": [ + "l1, l2, ltype = getInput(Lx,Ly)\n", + "print(l1)\n", + "print(l2)\n", + "print(ltype)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "------" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 3.) Misc" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Draw Lattice" + ] + }, + { + "cell_type": "code", + "execution_count": 149, + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [], + "source": [ + "\"\"\"\n", + "Function that draws a 2D Lattice\n", + "\"\"\"\n", + "def draw2DLattice(Lx,Ly):\n", + " n=Lx+2\n", + " m=Ly+2\n", + "\n", + " G = nx.grid_2d_graph(n,m) # Lx, Ly grid\n", + "\n", + " pos = dict( (n, n) for n in G.nodes() );\n", + " #labels = Dict( ((i, j), i + (Ly-j) * (Lx+1) ) for (i,j) in G.nodes() );\n", + " #edges=Dict()\n", + " #for (i,edge) in enumerate(G.edges)\n", + " # edges[edge]=i\n", + " # print(i,\"\\t\",edge,\"\\n\")\n", + " #end\n", + "\n", + "\n", + " for j in range(0,m-1):\n", + " G.remove_edge((0,j),(0,j+1))\n", + " G.remove_edge((n-1,j),(n-1,j+1))\n", + "\n", + " for i in range(0,n-1):\n", + " G.remove_edge((i,0),(i+1,0))\n", + " G.remove_edge((i,m-1),(i+1,m-1))\n", + "\n", + "\n", + "\n", + " edgescenter = dict()\n", + " countedge = 1\n", + " for j in range(1,m-1):\n", + " for i in range(1,n-1):\n", + " edgescenter[((i,j),(i+1,j))] = countedge\n", + " countedge += 1\n", + "\n", + "\n", + " for j in range(1,m-1):\n", + " for i in range(1,n-1):\n", + " edgescenter[((i,j),(i,j+1))] = countedge\n", + " countedge += 1\n", + "\n", + "\n", + " labelscenter = dict()\n", + " countlabels = 1\n", + " for j in range(1,m-1):\n", + " for i in range(1,n-1):\n", + " labelscenter[(i,j)]=countlabels\n", + " countlabels+=1\n", + "\n", + "\n", + " color_map=[]\n", + " countcolor=0\n", + " for i in range(0,n):\n", + " for j in range(0,m):\n", + " if i==0 or j==0 or i==n-1 or j==m-1:\n", + " color_map.append(\"white\")\n", + " else:\n", + " if(countcolor%2 == 0):\n", + " color_map.append(\"blue\")\n", + " else:\n", + " color_map.append(\"red\")\n", + " countcolor+=1\n", + " countcolor+=1\n", + "\n", + "\n", + "# nx.draw(G, pos=pos, labels=labelscenter, node_color=color_map, font_color=\"white\", ax=fig.add_subplot(111));\n", + "# nx.draw_networkx_edge_labels(G,pos,edge_labels=edgescenter,font_color=\"black\", ax=fig.add_subplot(111));\n", + " nx.draw(G, pos=pos, labels=labelscenter, node_color=color_map, font_color=\"white\");\n", + " nx.draw_networkx_edge_labels(G,pos,edge_labels=edgescenter,font_color=\"black\");\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Generate Input" + ] + }, + { + "cell_type": "code", + "execution_count": 132, + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [], + "source": [ + "\"\"\"\n", + "Function that generates the input\n", + "\"\"\"\n", + "def getInput(Lx,Ly):\n", + " n=Lx+2\n", + " m=Ly+2\n", + "\n", + "\n", + " labelscenter = dict()\n", + " countlabels = 1\n", + " for j in range(1,m-1):\n", + " for i in range(1,n-1):\n", + " labelscenter[(i,j)]=countlabels\n", + " countlabels+=1\n", + "\n", + " edgescenter = dict()\n", + " countedge = 1\n", + " for j in range(1,m-1):\n", + " for i in range(1,n-1):\n", + " edgescenter[((i,j),(i+1,j))] = countedge\n", + " countedge += 1\n", + " for j in range(1,m-1):\n", + " for i in range(1,n-1):\n", + " edgescenter[((i,j),(i,j+1))] = countedge\n", + " countedge += 1\n", + "\n", + "\n", + " l1 = []\n", + " l2 = []\n", + " itype = [1 for x in range(0,len(edgescenter))]\n", + "\n", + " for x in edgescenter:\n", + " if Lx + 1 not in x[1]:\n", + " l1.append(labelscenter[x[0]])\n", + " l2.append(labelscenter[x[1]])\n", + " else:\n", + " l1.append(labelscenter[x[0]])\n", + " if x[1][0] > Lx:\n", + " l2.append(labelscenter[(x[1][0]-Lx,x[1][1])])\n", + " elif x[1][1] > Ly:\n", + " l2.append(labelscenter[(x[1][0],x[1][1]-Ly)])\n", + "\n", + "# print(l1,\"\\n\")\n", + "# print(l2,\"\\n\")\n", + "# print(itype,\"\\n\")\n", + " return l1,l2,itype" + ] + }, + { + "cell_type": "code", + "execution_count": 133, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "([1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9],\n", + " [2, 3, 1, 5, 6, 4, 8, 9, 7, 4, 5, 6, 7, 8, 9, 1, 2, 3],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])" + ] + }, + "execution_count": 133, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "getInput(3,3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "header=r\"\"\"\n", + "\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python3", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/src/elem_diag.irp.f b/src/elem_diag.irp.f index 85b244f..96fc01e 100644 --- a/src/elem_diag.irp.f +++ b/src/elem_diag.irp.f @@ -18,10 +18,10 @@ ! if(.not.redo)write(6,*)'vijayyves' do i=1,nlientot if(yalt(i))then -! xmatd=-(xj1+xeneJ(i)*xbJ)+xmatd - xmatd= -(xjz(i))+xmatd - if(yw)write(6,*)xmatd,'xmatd' - if(yw)write(6,*)'xjz',xjz(i) +! xmatd=-(xj1+xeneJ(i)*xbJ)+xmatd + xmatd= -(xjz(i))+xmatd + if(yw)write(6,*)xmatd,'xmatd' + if(yw)write(6,*)'xjz',xjz(i) endif ! if(yrep1(i))then ! xmat=xv1+xmat diff --git a/src/ex1.c b/src/ex1.c index c2a2659..98f89f0 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -121,6 +121,7 @@ int main(int argc,char **argv) &getdata.fix_trou1, &getdata.fix_trou2, &getdata.FAM1, + &mpiid, tcol, val); // if(i%getdata.npar == 0 && mpiid==0){ @@ -128,9 +129,11 @@ int main(int argc,char **argv) // } for(ll=0;ll0) { ierr = PetscPrintf(PETSC_COMM_WORLD, - " k ||Ax-kx||/||kx|| \n" - " ----------------- ----------------- ------------------\n");CHKERRQ(ierr); + " k ||Ax-kx||/||kx|| \n" + " ----------------- ----------------- ------------------\n");CHKERRQ(ierr); for(i=0;i i=",i +! write(6,*)(ideter2(iik),iik=1,natom) + ideter2=0 + Touch deter + call adr(deter,iaa) + call elem_diag(xmat) + countcol+=1 +! print *,"1->id=",iaa," val=",xmat*1.0d0 + col(countcol)=iaa + val(countcol)=xmat*1.0d0 + + tistart2+=1 + detfound2+=countcolfull(count2) +! if(i.ne.1)then + count2+=1 +! endif + endif + imat4=iaa+1 + jmat4=foundetadr(i) + dmat4=foundetdmat(i) + if(jmat4.le.(nt1*nt2) .and. dmat4 .ne. 0d0)then + countcol+=1 +! print *,"2->id=",jmat4," val=",dmat4 + col(countcol)=jmat4 + val(countcol)=dmat4 + endif enddo +! do i=0,200 +! print *,i," det=",det(i,1)," add=",det(i,2) +! enddo + return end diff --git a/src/foundet.irp.f b/src/foundet.irp.f index e95c062..de6c850 100644 --- a/src/foundet.irp.f +++ b/src/foundet.irp.f @@ -1,7 +1,7 @@ use iso_c_binding BEGIN_PROVIDER[integer, foundet,(natomax,maxlien)] &BEGIN_PROVIDER[integer(C_SIZE_T), foundetadr,(maxlien)] -&BEGIN_PROVIDER[real, foundetdmat,(maxlien)] +&BEGIN_PROVIDER[real*8, foundetdmat,(maxlien)] &BEGIN_PROVIDER[integer(C_SIZE_T), foundadd,(maxlien,3)] &BEGIN_PROVIDER[integer(C_SIZE_T), foundaddh,(maxlien,3)] &BEGIN_PROVIDER[integer, detfound] diff --git a/src/get_s2.c b/src/get_s2.c index 081213e..7a82498 100644 --- a/src/get_s2.c +++ b/src/get_s2.c @@ -675,7 +675,8 @@ void get_s2(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *n *xymat2=*xymat2+xmat2; *xymat3=*xymat3+xmat3; *xymat4=*xymat4+xmat4; -// if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3); +// if(mpiid==0)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3); +// printf(" %d) ii = %d norm = %18f xymat = %18f 3 = %18f 4 = %18f\n",mpiid, ii, *norm, *xymat, *norm3, *xymat2, *xymat3); } ierr = PetscTime(&tt2); diff --git a/src/getdet.irp.f b/src/getdet.irp.f index acba73a..ba30c33 100644 --- a/src/getdet.irp.f +++ b/src/getdet.irp.f @@ -24,12 +24,12 @@ subroutine getdet(add,ideter) detb = det(ib,1) deta = deth(ia,1) if(FAM1) then - if(fix_trou1 .eq. fix_trou2) then - deta = ISHFT(deta,-(natom/2)) - else - natom2 = natom - (fix_trou2 - fix_trou1) - deta = ISHFT(deta, -natom2) - endif + if(fix_trou1 .eq. fix_trou2) then + deta = ISHFT(deta,-(natom/2)) + else + natom2 = natom - (fix_trou2 - fix_trou1) + deta = ISHFT(deta, -natom2) + endif endif ! do while (i.le.(ib)) ! const=1 @@ -50,7 +50,6 @@ subroutine getdet(add,ideter) ! i+=1 ! write(6,14)deta,deta ! enddo - const=0 if(FAM1) then if(fix_trou1 .eq. fix_trou2) then natom2 = natom/2 @@ -66,6 +65,12 @@ subroutine getdet(add,ideter) ideter((natom2)-i)=3 endif enddo +! do i=ib-2,ib+2 +! print *,i," det=",det(i,1)," add=",det(i,2) +! enddo +! print *,"add=",ib," detb=",detb +! write(6,*)(ideter(i),i=1,natom) + const=0 do i=0,natom-1 if(ideter(natom-i).eq.1)then if(BTEST(detb,const))then diff --git a/src/natom.irp.f b/src/natom.irp.f index 58bf864..d588596 100644 --- a/src/natom.irp.f +++ b/src/natom.irp.f @@ -94,33 +94,41 @@ BEGIN_PROVIDER [integer, natom] ! FAM1=.TRUE. yham=.TRUE. - write(6,*)'HAMILTONIEN t-J' - write(6,*)'Le nombre de trou est : ',ntrou - write(6,*)'Famille 1 : ',FAM1 - if(FAM1) then - if(fix_trou1 .ne. fix_trou2) write(6,*)'Trou fixe entre :', fix_trou1, "et ", fix_trou2 + if(mpiid==0)then + write(6,*)'HAMILTONIEN t-J' + write(6,*)'Le nombre de trou est : ',ntrou + write(6,*)'Famille 1 : ',FAM1 + if(FAM1) then + if(fix_trou1 .ne. fix_trou2) write(6,*)'Trou fixe entre :', fix_trou1, "et ", fix_trou2 + endif endif !--------------------------------------------- - write(6,*)' ' - write(6,*)' ' - write(6,*)'LECTURE DES ATOMES, DES LIAISONS, DES INTEGRALES' - write(6,*)' ' - write(6,*)' ' + if(mpiid==0)then + write(6,*)' ' + write(6,*)' ' + write(6,*)'LECTURE DES ATOMES, DES LIAISONS, DES INTEGRALES' + write(6,*)' ' + write(6,*)' ' + endif !-----------Lecture 1ier voisin ! read(5,clust) jclust=jclust+1 if(yw)write(6,*)' ' - write(6,*)'================ CLUSTER',jclust,'==================' - write(6,*)' ' + if(mpiid==0)then + write(6,*)'================ CLUSTER',jclust,'==================' + write(6,*)' ' + endif !------------------------ do i=1,maxlien if(l1(i).eq.0)EXIT nlien=i enddo + if(mpiid==0)then if(yw)write(6,*)(l1(i),i=1,maxlien) write(6,*)'Liaisons entre les atomes',nlien + endif do i=1,natomax y(i)=.false. @@ -140,18 +148,22 @@ BEGIN_PROVIDER [integer, natom] y(l1(i))=.true. y(l2(i))=.true. ltyp(nlientot+i)=ktyp(i) + if(mpiid==0)then write(6,*)'Les atomes ',l1(i),' et ',l2(i),' forment la liaison', & ilien(l1(i),l2(i)),'qui est de type',ltyp(i) if(yw)write(6,*)'iliatom1',iliatom1(i) + endif enddo nlientot=nlientot+nlien do i=1,natomax if(y(i))natom=natom+1 enddo - write(6,*)'=============================================' - write(6,*)'Le nombre total d atomes est ',natom - write(6,*)'=============================================' + if(mpiid==0)then + write(6,*)'=============================================' + write(6,*)'Le nombre total d atomes est ',natom + write(6,*)'=============================================' + endif !---------------------------------------------------------------- !----------------Voisins--------------------- @@ -210,24 +222,28 @@ BEGIN_PROVIDER [integer, natom] !C xjjxy=(/.1000d0,-0.8d0,0.000d0/) !C xtt=(/-1.0000d0,0.d0,0.0d0/) - write(6,*)'Nombre de J differents',ityp do ikl=1,nlientot Ykl(ikl)=.false. enddo - do il=1,nlientot - write(6,*)'type de liaison',il,ltyp(il) - enddo - do iiki=1,ityp - write(6,*)'type de J',xjjz(iiki) - enddo + if(mpiid==0)then + write(6,*)'Nombre de J differents',ityp + do il=1,nlientot + write(6,*)'type de liaison',il,ltyp(il) + enddo + do iiki=1,ityp + write(6,*)'type de J',xjjz(iiki) + enddo + endif do il=1,nlientot if(.not.ykl(il))then xjz(il)=xjjz(ltyp(il)) - write(6,*)'Parametres : Jz',il,'=',xjz(il) xjxy(il)=xjjxy(ltyp(il)) - write(6,*)'Parametres : Jxy',il,'=',xjxy(il) xt(il)=xtt(ltyp(il)) - write(6,*)'Parametre : t',il,'=',xt(il) + if(mpiid==0)then + write(6,*)'Parametres : Jz',il,'=',xjz(il) + write(6,*)'Parametres : Jxy',il,'=',xjxy(il) + write(6,*)'Parametre : t',il,'=',xt(il) + endif ykl(il)=.true. endif enddo @@ -248,26 +264,28 @@ BEGIN_PROVIDER [integer, natom] xenediagT=0.000d0 xspar=-0.00d0 xsperp=-0.00d0 - write(6,*)'coucoudslect3' - write(6,*)'coucou' - write(6,*)'Parametres pour le t-J' - write(6,*)'xj1 = ',xj1 - write(6,*)'xj2 = ',xj2 - write(6,*)'xt1 = ',xt1 - write(6,*)'xt2 = ',xt2 - write(6,*)'xv1 = ',xv1 - write(6,*)'xv2 = ',xv2 - write(6,*)'xv3 = ',xv3 - write(6,*)'xbj = ',xbj - write(6,*)'xbt = ',xbt - write(6,*)'xeneparJ = ',xeneparJ - write(6,*)'xeneperpJ = ',xeneperpJ - write(6,*)'xeneparT = ',xeneparT - write(6,*)'xeneperpT = ',xeneperpT - write(6,*)'xenediagJ = ',xenediagJ - write(6,*)'xenediagT = ',xenediagT - write(6,*)'xspar = ',xspar - write(6,*)'xsperp = ',xsperp + if(mpiid==0)then + write(6,*)'coucoudslect3' + write(6,*)'coucou' + write(6,*)'Parametres pour le t-J' + write(6,*)'xj1 = ',xj1 + write(6,*)'xj2 = ',xj2 + write(6,*)'xt1 = ',xt1 + write(6,*)'xt2 = ',xt2 + write(6,*)'xv1 = ',xv1 + write(6,*)'xv2 = ',xv2 + write(6,*)'xv3 = ',xv3 + write(6,*)'xbj = ',xbj + write(6,*)'xbt = ',xbt + write(6,*)'xeneparJ = ',xeneparJ + write(6,*)'xeneperpJ = ',xeneperpJ + write(6,*)'xeneparT = ',xeneparT + write(6,*)'xeneperpT = ',xeneperpT + write(6,*)'xenediagJ = ',xenediagJ + write(6,*)'xenediagT = ',xenediagT + write(6,*)'xspar = ',xspar + write(6,*)'xsperp = ',xsperp + endif !=================================================================== !==================================================================== ! calcul des plaquettes pour un 2D t-J @@ -315,11 +333,13 @@ BEGIN_PROVIDER [integer, natom] enddo enddo enddo - write(6,*)'Le systeme comporte ',nplac,' plaquettes.' - do kko=1,nplac - write(6,*)'La plaquette ',kko,' est contituee des atomes',& - iplac(1,kko),' ',iplac(2,kko),' ',iplac(3,kko),' et ',iplac(4,kko) - enddo + if(mpiid==0)then + write(6,*)'Le systeme comporte ',nplac,' plaquettes.' + do kko=1,nplac + write(6,*)'La plaquette ',kko,' est contituee des atomes',& + iplac(1,kko),' ',iplac(2,kko),' ',iplac(3,kko),' et ',iplac(4,kko) + enddo + endif !=================================================================== ! isz=0 IPREC=8 @@ -332,6 +352,7 @@ BEGIN_PROVIDER [integer, natom] xseuil=1.0E-008 ysuiv=.FALSE. yec=.TRUE. + if(mpiid==0)then write(6,*)'Spin total',isz write(6,*)'Nombre de vecteurs demande',nvec write(6,*)'Nombre maximal d iterations de Davidson',niter @@ -344,12 +365,15 @@ BEGIN_PROVIDER [integer, natom] ! write(6,*)Emin ! write(6,*)Emax ! write(6,*)M0 + endif if(yham)then IAL0=(natom-ntrou)/2+mod(natom-ntrou,2)+isz else IAL0=(natom)/2+mod(natom,2)+isz endif - write(6,*)'=======nombre de centres de spin alpha=====',ial0 + if(mpiid==0)then + write(6,*)'=======nombre de centres de spin alpha=====',ial0 + endif natrest=natom-ntrou !C calculating nalpha and nbeta diff --git a/src/nt1.irp.f b/src/nt1.irp.f index 8bd0124..1832b17 100644 --- a/src/nt1.irp.f +++ b/src/nt1.irp.f @@ -16,5 +16,7 @@ BEGIN_PROVIDER [integer(C_SIZE_T), nt1] endif endif nt1= nint(gamma(1.0*(natom2+1))/(gamma(1.0*(natom2-ntrou+1))*gamma(1.0*(ntrou+1))),selected_int_kind(16)) - write(6,*)'nt1',nt1 + if(mpiid==0)then + write(6,*)'nt1',nt1 + endif END_PROVIDER diff --git a/src/nt2.irp.f b/src/nt2.irp.f index 0ff4c5d..51fee91 100644 --- a/src/nt2.irp.f +++ b/src/nt2.irp.f @@ -7,5 +7,7 @@ BEGIN_PROVIDER [integer(C_SIZE_T), nt2] ! call combin(idet2(1,nt2+1),natrest,ial0,nt2,32,jrangmax) nt2= nint(gamma(1.0*(natom-ntrou+1))/((gamma(1.0*(nalpha+1))*gamma(1.0*(nbeta+1)))),selected_int_kind(16)) - print *,"nt2=",nt2 + if(mpiid==0)then + print *,"nt2=",nt2 + endif END_PROVIDER diff --git a/src/providdet.irp.f b/src/providdet.irp.f index d651244..afac84c 100644 --- a/src/providdet.irp.f +++ b/src/providdet.irp.f @@ -8,7 +8,7 @@ BEGIN_PROVIDER[integer(C_SIZE_T),det,(nt2,2)] implicit none ! integer(kind=selected_int_kind(16))::dethsh integer(C_SIZE_t)::a - integer(C_SIZE_T)::i,count + integer(C_SIZE_T)::i,count,iik integer::const i=1 a=0 @@ -18,14 +18,7 @@ BEGIN_PROVIDER[integer(C_SIZE_T),det,(nt2,2)] If(ntrou.ge.1)then const=0 -! dethsh = ISHFT(deth,-natom/2) -! i=nt1 do while (i.le.(nt1)) -! if(a.eq.dethsh)then -! addh=i-1 -! EXIT -! endif - i+=1 a+=1 do while(popcnt(a).ne.ntrou) @@ -42,21 +35,16 @@ BEGIN_PROVIDER[integer(C_SIZE_T),det,(nt2,2)] deth(count, 1) = a endif deth(count,2)=i-1 -! write(6,16)ISHFT(a,natom/2),ISHFT(a,natom/2),i-1 enddo -! if(a.eq.dethsh )then -! count+=1 -! deth(1,1)=ISHFT(a,natom/2) -! deth(1,2)=nt1 -! endif - endif !C if det=0 then exit a=0 i=0 count=0 - print *,'nt2=',nt2,'nbeta=',nbeta + if(mpiid==0)then + print *,'nt2=',nt2,'nbeta=',nbeta + endif do while (i.lt.(nt2)) i+=1 @@ -79,6 +67,9 @@ BEGIN_PROVIDER[integer(C_SIZE_T),det,(nt2,2)] endif ! write(6,16)a,a,i enddo +! do i=0,nt2 +! print *,i," det=",det(i,1)," add=",det(i,2) +! enddo 10 FORMAT(B64,I8,F8.2) diff --git a/src/provide_mpiid.irp.f b/src/provide_mpiid.irp.f new file mode 100644 index 0000000..8786ff4 --- /dev/null +++ b/src/provide_mpiid.irp.f @@ -0,0 +1,7 @@ +use iso_c_binding +BEGIN_PROVIDER [integer(C_SIZE_T), mpiid] + BEGIN_DOC + ! simply store the mpiid for pinting + END_DOC +END_PROVIDER + diff --git a/src/stimsyr.h b/src/stimsyr.h index 5964517..798c85f 100644 --- a/src/stimsyr.h +++ b/src/stimsyr.h @@ -17,5 +17,6 @@ void unit_l1_( long int *, _Bool *, long int *, + long int *, double *); diff --git a/src/unit_FIL44.irp.f b/src/unit_FIL44.irp.f index 96da0c8..46a7ddd 100644 --- a/src/unit_FIL44.irp.f +++ b/src/unit_FIL44.irp.f @@ -88,6 +88,7 @@ enddo ! print *,"tistart =", tistart,"countcol =", countcol,"\n",(tcountcol(i),i=1,nrows) ! print *,"" -! print *,(tcol(i),i=1,maxlien) +! print *,(tval(i),i=1,maxdet) +! print *,(tcol(i),i=1,maxdet) deallocate(ideter2) end diff --git a/src/unit_l1.irp.f b/src/unit_l1.irp.f index fc0dc40..b897a01 100644 --- a/src/unit_l1.irp.f +++ b/src/unit_l1.irp.f @@ -13,10 +13,11 @@ tfix_trou1, & tfix_trou2, & tfam1, & + tmpiid, & tcol,tval) use iso_c_binding implicit none - integer,INTENT(INOUT)::tistart, tnrows + integer,INTENT(INOUT)::tistart, tnrows, tmpiid integer,INTENT(INOUT)::tntrou, tisz integer,INTENT(INOUT)::tfix_trou1, tfix_trou2 logical*1,INTENT(INOUT)::tfam1 @@ -45,7 +46,8 @@ fix_trou2 = tfix_trou2 tcol=0 tval=0d0 - provide l1 l2 ktyp xtt xjjxy xjjz ntrou + mpiid=tmpiid + provide l1 l2 ktyp xtt xjjxy xjjz ntrou mpiid !print *,"l1" !print *,l1 !print *,"xjjz"