diff --git a/.gitignore b/.gitignore new file mode 100755 index 0000000..6748628 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.pyc +.ipynb_* diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..140d2b6 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,8 @@ +include COPYING +include COPYING.LESSER +include MANIFEST.in +include README.rst +include requirements.txt +include setup.py + +recursive-include fermilibpluginpsi4 *.py _psi4_template* diff --git a/README.rst b/README.rst index f20f7cb..5519597 100644 --- a/README.rst +++ b/README.rst @@ -1,7 +1,7 @@ -FermiLib Plugin to interface with Psi4 +FermiLib plugin to interface with Psi4 ====================================== -This repository contains a plugin which allows `FermiLib `__ to interface with the open-source quantum chemistry package `Psi4 `__. +This repository contains a plugin which allows `FermiLib `__ (licensed under Apache 2) to interface with the open-source quantum chemistry package `Psi4 `__ (licensed under GNU Lesser General Public License version 3). License ------- diff --git a/examples/psi4_demo.ipynb b/examples/psi4_demo.ipynb new file mode 100755 index 0000000..491c52e --- /dev/null +++ b/examples/psi4_demo.ipynb @@ -0,0 +1,208 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Running psi4 to populate the MolecularData class\n", + "\n", + "The module run_psi4.py provides a user-friendly way of running psi4 calculations in FermiLib. The basic idea is that once one generates a MolecularData instance, one can then call psi4 with a specification of certain options (for instance, how much memory to use and what calculations to do) in order to compute things about the molecule, update the MolecularData object, and save results of the calculation. The most common calculations users will want in FermiLib would probably be self-consistent field (aka Hartree-Fock calculations). Our code uses these \"SCF\" calculations compute orbitals, integrals, Hartree-Fock energy, and more. Other common calculations are CISD and FCI calculations which also compute the 1-RDM and 2-RDM associated with their answers, CCSD calculations which also compute the CCSD amplitudes (useful for UCC) and MP2 calculations which can also be used to UCCSD initial guesses. \n", + "Note that the \"delete_input\" and \"delete_output\" options indicate whether to save automatically generated psi4 input files or not.\n", + "\n", + "To use this plugin, you will need to personally download psi4. Python code in psi4_template instructs psi4 (not python) to load the MolecularData object, populate it with information from calculations and then save the MolecularData object as an HDF5. The module run_psi4 uses subprocess to actually run_psi4 and then loads the pickled MolecularData. Let us look at a simple example where we compute the energy of H$_2$ at various bond lengths.\n", + "\n", + "Warnings: electronic structure calculations are finicky. They sometimes fail for surprising reasons. If a particular calculation is not converging it is probably necessary to change some of the SCF options set in psi4_template. See the Psi4 documentation for more information or consult and electronic structure theory expert.\n", + "\n", + "Note finally that there is also a script to use pyscf in this repo. Integration with pyscf is not complete so use at your own risk." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "At bond length of 0.2 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of 0.1641558779392942 Hartree.\n", + "MP2 energy of 0.15902385897340335 Hartree.\n", + "FCI energy of 0.1574821241424691 Hartree.\n", + "Nuclear repulsion energy between protons is 2.64588604295 Hartree.\n", + "Spatial orbital 0 has energy of -0.85884365079039 Hartree.\n", + "Spatial orbital 1 has energy of 1.572449495506951 Hartree.\n", + "\n", + "At bond length of 0.4 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -0.9043953123281971 Hartree.\n", + "MP2 energy of -0.911467754694884 Hartree.\n", + "FCI energy of -0.9141497082141274 Hartree.\n", + "Nuclear repulsion energy between protons is 1.322943021475 Hartree.\n", + "Spatial orbital 0 has energy of -0.7452464479052874 Hartree.\n", + "Spatial orbital 1 has energy of 1.167776648029189 Hartree.\n", + "\n", + "At bond length of 0.6000000000000001 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -1.1011793445775693 Hartree.\n", + "MP2 energy of -1.1113798775246235 Hartree.\n", + "FCI energy of -1.1162860078265324 Hartree.\n", + "Nuclear repulsion energy between protons is 0.8819620143166668 Hartree.\n", + "Spatial orbital 0 has energy of -0.640927365839849 Hartree.\n", + "Spatial orbital 1 has energy of 0.8382934397898524 Hartree.\n", + "\n", + "At bond length of 0.8 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -1.110918875412235 Hartree.\n", + "MP2 energy of -1.125517302745064 Hartree.\n", + "FCI energy of -1.1341476663578607 Hartree.\n", + "Nuclear repulsion energy between protons is 0.6614715107375 Hartree.\n", + "Spatial orbital 0 has energy of -0.5545643581907187 Hartree.\n", + "Spatial orbital 1 has energy of 0.612776410591467 Hartree.\n", + "\n", + "At bond length of 1.0 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -1.0661936970704187 Hartree.\n", + "MP2 energy of -1.086742793881516 Hartree.\n", + "FCI energy of -1.1011503293035885 Hartree.\n", + "Nuclear repulsion energy between protons is 0.52917720859 Hartree.\n", + "Spatial orbital 0 has energy of -0.4845267279502007 Hartree.\n", + "Spatial orbital 1 has energy of 0.45764406664496715 Hartree.\n", + "\n", + "At bond length of 1.2000000000000002 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -1.005201854187592 Hartree.\n", + "MP2 energy of -1.0337457589181585 Hartree.\n", + "FCI energy of -1.0567407451325903 Hartree.\n", + "Nuclear repulsion energy between protons is 0.4409810071583334 Hartree.\n", + "Spatial orbital 0 has energy of -0.4265977900276036 Hartree.\n", + "Spatial orbital 1 has energy of 0.34424909086732186 Hartree.\n", + "\n", + "At bond length of 1.4000000000000001 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -0.9415807727061606 Hartree.\n", + "MP2 energy of -0.980650948280493 Hartree.\n", + "FCI energy of -1.0154682481426915 Hartree.\n", + "Nuclear repulsion energy between protons is 0.3779837204214286 Hartree.\n", + "Spatial orbital 0 has energy of -0.3774229438915012 Hartree.\n", + "Spatial orbital 1 has energy of 0.25900471560908955 Hartree.\n", + "\n", + "At bond length of 1.6 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -0.8818353849537497 Hartree.\n", + "MP2 energy of -0.9343173778082237 Hartree.\n", + "FCI energy of -0.9834727280932716 Hartree.\n", + "Nuclear repulsion energy between protons is 0.33073575536875 Hartree.\n", + "Spatial orbital 0 has energy of -0.3353992875908116 Hartree.\n", + "Spatial orbital 1 has energy of 0.19468565861959788 Hartree.\n", + "\n", + "At bond length of 1.8 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -0.8289518481024154 Hartree.\n", + "MP2 energy of -0.8979454626073924 Hartree.\n", + "FCI energy of -0.961816952119916 Hartree.\n", + "Nuclear repulsion energy between protons is 0.29398733810555555 Hartree.\n", + "Spatial orbital 0 has energy of -0.2996669223412511 Hartree.\n", + "Spatial orbital 1 has energy of 0.14603876496061327 Hartree.\n", + "\n", + "At bond length of 2.0 Bohr, molecular hydrogen has:\n", + "Hartree-Fock energy of -0.7838924536817647 Hartree.\n", + "MP2 energy of -0.8725561380421281 Hartree.\n", + "FCI energy of -0.9486411117424073 Hartree.\n", + "Nuclear repulsion energy between protons is 0.264588604295 Hartree.\n", + "Spatial orbital 0 has energy of -0.2695590236598946 Hartree.\n", + "Spatial orbital 1 has energy of 0.10906868098620753 Hartree.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8XNWZ//HPo96LJVmWLMuW5YJl44YwYEjoCSHJmnRa\nlg1svCQQyGZTSFnCK9lN2OSX3U1IWEKADSwlIZTgAAlgwKQY27hhXHDHtmzZlmRZ3arP7497ZY2k\nkTQjaebOSM/79ZrXzNy5M/exGPTVueeec0RVMcYYYwIV43UBxhhjoosFhzHGmKBYcBhjjAmKBYcx\nxpigWHAYY4wJigWHMcaYoFhwGGOMCYoFhzHGmKBYcBhjjAlKnNcFhEJubq5OmzbN6zKMMSZqbNiw\noVpV8wLZ19PgEJErgJ8CscADqnp3n9evA74BCNAAfEFV3x7qc6dNm8b69etDULExxoxNInIg0H09\nO1UlIrHAL4APAWXANSJS1me3/cCFqnom8H3g/vBWaYwxpi8v+ziWAHtUdZ+qtgG/AZb57qCqq1W1\n1n26BigKc43GGGP68DI4JgOHfJ5XuNsGchPwx5BWZIwxZkhR0TkuIhfjBMcFg+yzHFgOUFxcHKbK\njDFm/PGyxXEYmOLzvMjd1ouIzAceAJapas1AH6aq96tquaqW5+UFdGGAMcaYYfAyON4CZopIiYgk\nAFcDK3x3EJFi4Bngs6q6K6TVbHkS/mse3JXl3G95MqSHM8aYaOVZcKhqB3Ar8BKwA3hSVbeJyM0i\ncrO7251ADnCviGwWkZBcY7vyt/fQ+dxtUHcIUKg7ROdzt7Hyt/eE4nDGGBPVZCwuHVteXq7BjOM4\n9eM5JDUd6b89tZCkr+0YzdKMMSYiicgGVS0PZF+bcgRIaqoMarsxxoxnFhwAmQMMDxlouzHGjGMW\nHACX3klnbHKvTZ2xyXDpnR4VZIwxkcuCA1idegn/2vV5GmKzAGhLyuVfuz7P6tRLPK7MGGMijwUH\nsKWijo9cfzsbr3wBgGNn/hMfuf52tlTUeVyZMcZEnqgYOR5qN19YCsDx+jQOaw5tBzew9MNfZ2lp\nrseVGWNM5LEWh4+JGUnsiplJxoktXpdijDERy4Kjj9rseeS1H4HmE16XYowxEcmCo4+YorMAaNr/\nlseVGGNMZLLg6GPi7HMBqN75pseVGGNMZLLg6GNuyRT2dhXQdXij16UYY0xEsuDoIzMlnr3xs8g+\nudXrUowxJiJZcPhRnzOfrM4aqO8/8aExxox3Fhx+JEw9G4CTe9Z4XIkxxkQeCw4/Js9eQrvGUrvL\ngsMYY/qy4PCjrHgiO3UKMZWbvC7FGGMijgWHH8kJsRxMmk1u/XYYgwtdGWPMSFhwDKA5dwGp2ojW\n7PW6FGOMiSgWHANILlkCQM0uGwhojDG+LDgGMPWMxbRoAvV713pdijHGRBQLjgHMLsxmOyXEH9vs\ndSnGGBNRLDgGEB8bw+GUOUxs3AWd7V6XY4wxEcPT4BCRK0Rkp4jsEZE7/LwuIvIz9/UtIrI4nPW1\nTlxIIq10HtsezsMaY0xE8yw4RCQW+AXwIaAMuEZEyvrs9iFgpntbDvxPOGtML3U6yG2mXGOM6eFl\ni2MJsEdV96lqG/AbYFmffZYBj6hjDZAlIgXhKnDG7DM5qam2NocxxvjwMjgmA4d8nle424LdBwAR\nWS4i60VkfVVV1agUOD0vne2Uklz19qh8njHGjAVjpnNcVe9X1XJVLc/LyxuVz4yJEY6llzGxZS+0\nt4zKZxpjTLTzMjgOA1N8nhe524LdJ6Q6CxYTSxftFXZZrjHGgLfB8RYwU0RKRCQBuBpY0WefFcDf\nu1dXnQvUqWplOIvMmnEOAMetg9wYYwCI8+rAqtohIrcCLwGxwEOquk1EbnZfvw94EbgS2AM0A58L\nd51nzJrF0RezaT2wPtyHNsaYiORZcACo6os44eC77T6fxwrcEu66fE3OSuYNmUFZzRYvyzDGmIgx\nZjrHQ0VEqMmay8S2Q9By0utyjDHGcxYcgZhcDkCLna4yxhgLjkDkzXY6yG0EuTHGeNzHES3Kpk9l\nf1c+VGz0uhRjjPGctTgCkJuWyO64WWTVvuN1KcYY4zkLjgDVTTiT7I4qaDjqdSnGGOMpC44AxU1x\nOsgb963zuBJjjPGWBUeACmYvoUNjqNm1xutSjDHGUxYcASqbNondWgRHrIPcGDO+WXAEKCMpnn2J\ns8ip2wqqXpdjjDGeseAIQlPOAtK6GqB2v9elGGOMZyw4gpA49WwATu5Z63ElxhjjHQuOIBSdcRan\nNN6CwxgzrllwBGFuUQ47dBqxlZu8LsUYYzxjwRGEpPhYDiafwcTGHdDZ4XU5xhjjCQuOIJ2auJBE\nbUWr3vW6FGOM8YQFR5DSSpwO8mobCGiMGacsOIJUMnsB9ZpsU48YY8YtC44gzZqUwTamk3j8ba9L\nMcYYT1hwBCkuNobK1LlMbN4N7ae8LscYY8LOgmMY2ictJI5OOiq3eF2KMcaEnSfBISITROQVEdnt\n3mf72WeKiLwuIttFZJuI3O5Frf5kzeheStY6yI0x449XLY47gFdVdSbwqvu8rw7gX1S1DDgXuEVE\nysJY44BmzphNlWbS8t5bXpdijDFh51VwLAMedh8/DFzVdwdVrVTVje7jBmAHMDlsFQ5iWm4a26SU\n1Go7VWWMGX+8Co58Va10Hx8F8gfbWUSmAYuAiJgkKiZGqMqYR27rAThV73U5xhgTVgEHh4ikBPPB\nIrJSRLb6uS3z3U9VFRhwgQsRSQOeBr6sqgP+lhaR5SKyXkTWV1VVBVPqsGjhYmJQ2g7Zwk7GmPFl\nyOAQkaUish14132+QETuHep9qnqZqs7zc3sOOCYiBe7nFQDHBzh2PE5oPKaqzwxxvPtVtVxVy/Py\n8oYqb8RyZjod5FU73wz5sYwxJpIE0uL4L+CDQA2Aqr4NvH+Ex10B3OA+vgF4ru8OIiLAg8AOVf3P\nER5v1M0pLeFgVx7thzZ4XYoxxoRVQKeqVPVQn02dIzzu3cDlIrIbuMx9jogUisiL7j7nA58FLhGR\nze7tyhEed9QUZCbxbuwsMk5YB7kxZnyJC2CfQyKyFFD31NHtOFc4DZuq1gCX+tl+BLjSffxXQEZy\nnFASEU5mz2PCib9BYxWkhf70mDHGRIJAWhw3A7fgXAp7GFjoPh/3YorOAqDlgI3nMMaMH0MGh6pW\nq+p1qpqvqhNV9Xq3xTDu5c9aQqcK1dZBbowZRwK5qmqWiLwqIlvd5/NF5DuhLy3yzS2ZzB6dTFeF\nXZJrjBk/AjlV9Svgm0A7gKpuAa4OZVHRYkJqAnviZzGh7h3QAYeiGGPMmBJIcKSoat9Vi2zBbVdD\nznzSO+vg5EGvSzHGmLAIJDiqRaQUd3S3iHwSqBz8LeNHQnE5AA22IqAxZpwIJDhuAX4JnCEih4Ev\n41xpZYDJs8tp1Thqd9sU68aY8WHQcRwiEgOUq+plIpIKxLgz1RrX3OI83tViJlRu8roUY4wJi0Fb\nHKraBXzdfdxkodFfWmIc7yWeQV79duga6YB6Y4yJfIGcqlopIl91V+Sb0H0LeWVRpDlvAUnaglbv\n8roUY4wJuUCmHPmMe+87WlyB6aNfTnRKKTkbjsDJPWvJnjjH63KMMSakAmlxzFHVEt8bEBFLuEaK\nqbMW0qhJ1O+JiHWmjDEmpAIJjtUBbhu35kzOYqtOJ/7YZq9LMcaYkBvwVJWITMKZ2DBZRBbRM1Nt\nBhDUaoBjXWJcLIdT5nBW0wroaIO4BK9LMsaYkBmsj+ODwD8ARcBP6AmOeuBboS0r+rTlLyD+wNN0\nHd1KTNFir8sxxpiQGTA4VPVhEfk/4BpVfSyMNUWl9NJz4ADU7H6TPAsOY8wYFsg4jn8OUy1RbebM\nMmo0neb9tjaHMWZss3Eco2RGfjpbmUHycesgN8aMbTaOY5TExgjH08q4oOkJaG2ExDSvSzLGmJAI\nZAXAEj83Cw0/OgsWEUsXHYdt3ipjzNgVSIsDEZmHM+gvqXubqj4SqqKiVdasc2EPVO1cQ8H093ld\njjHGhEQgS8d+F7jHvV0M/Aj4u5Ec1O0neUVEdrv32YPsGysim0Tk+ZEcMxzmlE6nQnNpO2Ad5MaY\nsSuQzvFPApcCR1X1c8ACIHOEx70DeFVVZwKvus8HcjuwY4THC4viCSnskBmk1WzxuhRjjAmZQIKj\nxb0st0NEMoDjwJQRHncZ8LD7+GHgKn87iUgR8GHggREeLyxEhJrMeeS0V0JTjdflGGNMSAQSHOtF\nJAv4FbAB2Ai8OcLj5qtq9/KzR4H8Afb7b5z1QLpGeLywkcnO4L/Wg+s9rsQYY0JjyM5xVf2i+/A+\nEfkTkKGqQ56LEZGVwCQ/L327z+eriKif938EOK6qG0TkogCOtxxYDlBcXDzU7iGTO+scurYJ1bve\nZPKcD3pWhzHGhMpgkxwOOG+GiCxW1Y2DfbCqXjbI+4+JSIGqVopIAc7pr77OB/5ORK7EuZorQ0Qe\nVdXrBzje/cD9AOXl5f2CKFzmTS9irxaSeGiDVyUYY0xIDdbi+InP47NwTlN1U+CSERx3BXADcLd7\n/1zfHVT1m8A3AdwWx1cHCo1Ikp+RxPq4mbyv9m1QBZGh32SMMVFksEkOL+5+LCKbfJ+PgruBJ0Xk\nJuAA8Gn3OIXAA6p65SgeK+zqss8ko2YV1B+GzCKvyzHGmFEV0ABAnBbGqFHVGpxLfPtuPwL0Cw1V\nXQWsGs0aQil+yllQA0373yJ1oQWHMWZsCeSqKhOkSbPPpl1jqd29xutSjDFm1A3WOX4PPS2NIhH5\nme/rqnpbKAuLZmdOncgOLWbCYesgN8aMPYOdqvIdiGC/AYOQlZLAGwmzmVn/F+jqghhr2Bljxo5B\nVwAMZyFjTVPufJKP/glO7IXcmV6XY4wxo8b+FA6RxKnlANTtXetxJcYYM7osOEKkePYimjWROusg\nN8aMMRYcITK3aAJbtYS4o7aokzFmbBlyHIeI5AGfB6b57q+qN4aurOiXkhDHoeQzWNj0InS2Q2y8\n1yUZY8yoCKTF8RzO+hsrgRd8bmYIp/IWkKBt6LFtXpdijDGjJpCR4ymq+o2QVzIGpU5fAhVQu2ct\nEwoXel2OMcaMikBaHM+7M9SaIJXOnEetptFoV1YZY8aQQILjdpzwaBGRehFpEJH6UBc2FswuyOAd\nLSXx+Ntel2KMMaNmyOBQ1XRVjVHVZFXNcJ9nhKO4aJcQF0Nl2hxyW/ZBW7PX5RhjzKgYMDhE5Az3\nfrG/W/hKjG4dkxYSSxedR6zVYYwZGwbrHP8KzlKsP/Hz2kgXcho3MkvPhX1wYtca8qad53U5xhgz\nYoPNVbXcvR/NBZzGndkzZnDk5Ql0HFjndSnGGDMqbOR4iE3PS2MbpaRWb/G6FGOMGRUWHCEWGyNU\nZ8wlp7UCWmq9LscYY0bMgiMcCp1rCdoPbfS4EGOMGbkhg0NEnhGRD4uIhcwwTZh5LgA1O9/0uBJj\njBm5QMLgXuBaYLeI3C0is0Nc05hTNn0Ke7sKaD+0fuidjTEmwgUyAHClql4HLAbeA1aKyGoR+ZyI\n2JSvASjKTmZn7AwyTrzjdSnGGDNiAZ1+EpEc4B+AfwQ2AT/FCZJXhnNQEZkgIq+IyG73PnuA/bJE\n5CkReVdEdohIVA6EEBFqs+aR2VEN9ZVel2OMMSMSSB/Hs8BfgBTgo6r6d6r6W1X9EpA2zOPeAbyq\nqjOBV93n/vwU+JOqngEsAHYM83ieiyk6C4DWg295XIkxxoxMIC2On6lqmar+UFV7/bmsquXDPO4y\n4GH38cPAVX13EJFM4P3Ag+6x2lT15DCP57lJs5bQoTHWQW6MiXqBrMeRLSIf77OtDnhHVY8P87j5\nPiF0FMj3s08JUAX8r4gsADYAt6tq0zCP6al50yaxU6eQedguyTXGRLdAWhw3AQ8A17m3XwHfAP4m\nIp8d6E0islJEtvq5LfPdT1UVZ+6rvuJw+lH+R1UXAU0MfEoLEVkuIutFZH1VVVUA/6zwyktPZE/c\nLLJPbgP19881xpjoEEiLIx6Yo6rHAEQkH3gEOAf4M/B//t6kqpcN9IEickxEClS1UkQKAH8tlwqg\nQlW7V0F6ikGCQ1XvB+4HKC8vj8jfzA05Z5Ja9Qqc2Ac5pV6XY4wxwxJIi6OoOzRcx4EpqnoCaB/m\ncVcAN7iPb8BZ17wXVT0KHPIZN3IpsH2Yx4sI8cVOl1DTfpvw0BgTvQIJjlUi8ryI3CAi3b/kV4lI\nKjDczuq7gctFZDdwmfscESkUkRd99vsS8JiIbAEWAj8Y5vEiwpTZi2nRBGp3r/G6FGOMGbZATlXd\nAnwcuMB9/gjwtNs3Mawp11W1BqcF0Xf7EeBKn+ebgeFeuRVx5k7JZZtOo7Byk9elGGPMsA0aHCIS\nC6x01+R4OjwljV2ZyfEcSJzNmQ0vQ2cHxAaS28YYE1kGPVWlqp1AlzumwoyC5tz5JGorVL3rdSnG\nGDMsgfzJ2wi8IyKv4FwSC4Cq3hayqsawlJIlUAn1e9eQMWme1+UYY0zQAgmOZ9ybGQXTZs2n7m8p\n1O9dR8b5/+h1OcYYE7Qhg0NVHxaRZKBYVXeGoaYxbe7kTNbrdGYe2+x1KcYYMyyBTHL4UWAz8Cf3\n+UIRWRHqwsaqpPhYDqfMIadpD7Sf8rocY4wJWiDjOO4CluCO2XAvkZ0ewprGvLb8hcTRiR7d4nUp\nxhgTtECCo11V6/ps6wpFMeNFeuk5AJzYZTPlGmOiTyDBsU1ErgViRWSmiNwDrA5xXWPazBmzOKZZ\nNO+3tTmMMdEnkOD4EjAXaAWeAOqBL4eyqLFuVn46W7WU5Oq3vS7FGGOCFshVVc3At92bGQXxsTEc\nSy8jt+n/4FQdJNn4SmNM9AjkqqpZInK/iLwsIq9138JR3FjWWbDIua+whZ2MMdElkAGAvwPuw1nM\nqTO05Ywf2TPPgT1Qs3sNE2cMa65IY4zxRCDB0aGq/xPySsaZOdOn8V5XPnEH1ntdijHGBCWQzvE/\niMgXRaRARCZ030Je2RhXkpPKdiklrcbGchhjoksgLY7ulfq+5rNNsUGAIxITI5zInEdW/WpoOAbp\n+V6XZIwxAQnkqqqScBQyLhWdBduh/dB64ss+7HU1xhgTkAFPVYnI130ef6rPa1G9hGukyJ91Np0q\n1OyypWSNMdFjsD6Oq30ef7PPa1eEoJZxZ+60QnZpEZ0VG7wuxRhjAjZYcMgAj/09N8NQkJnErtiZ\nZNVuBVWvyzHGmIAMFhw6wGN/z80wiAh12WeS2lkHJw94XY4xJlpteRL+ax7cleXcb3kypIcbLDgW\niEi9iDQA893H3c/PHMlB3Ut6XxGR3e599gD7/bOIbBORrSLyhIgkjeS4kSiu+GwATr1nEx4aY4Zh\ny5Pwh9ug7hCgzv0fbgtpeAwYHKoaq6oZqpquqnHu4+7n8SM87h3Aq6o6E3jVfd6LiEwGbgPKVXUe\nEEvvfpcxoXD2Ylo1nhO7rYPcGBOg1kY48Ca8eS/84XZob+n9ensLvPq9kB0+kHEcobAMuMh9/DCw\nCviGn/3igGQRaQdSgCPhKC6c5hfnsV2nkn/Y5qwyxvjR1gzHtsKRTT236l2gQyyLVFcRspK8Co58\nVa10Hx8F+o1+U9XDIvL/gINAC/Cyqr4cxhrDYkJqAq/Fz6Ks/jXo6oSYWK9LMsZ4paO1T0hshuM7\nQN1pAlPzoHAxlF0FhYugcCE8cJl7mqqPzKKQlRmy4BCRlcAkPy/1mp5dVVVE+nW2u/0ey4ASnGVr\nfyci16vqowMcbzmwHKC4uHiE1YdXU+58Eo+9CFU7Ib/M63KMMeHQ0QZVO3q3JI5th6525/XkCTB5\nMcy6wg2JRZBRCNLnotZL73T6NHxPV8UnO9tDJGTBoaqXDfSaiBwTkQJVrRSRAuC4n90uA/arapX7\nnmeApYDf4FDV+4H7AcrLy6Pqqq+kqWfDMWjct440Cw5joseWJ52+hLoK5y/8S++E+Z/uv19nB1S9\n2ycktkFnq/N6UqYTDEtv7QmJzCn9Q8Kf7uMFUsco8epU1QqcObDudu+f87PPQeBcEUnBOVV1KTAm\np5ItnrWA+rXJNOxdS9p5/+B1OcaYQHRfzdT9l3731UzaBQULe4fE0Xegw90vMQMKFsA5y3tCIrsk\nsJAYyPxPhzQo+vIqOO4GnhSRm4ADwKcBRKQQeEBVr1TVtSLyFLAR6AA24bYoxpozp2SzRUsoPbrZ\n61KMMYF69Xv+r2Z69mZOD3WLT3VCovzGnpCYMB1iApmYPHJ5EhyqWoPTgui7/Qhwpc/z7wLfDWNp\nnkhLjONg0hksaVzhdI7FJXpdkjHGn64up1/iwGr/HdIAKFx1nxMSuTPH5AUvXrU4TB+nJi4kruIZ\n9Og7SFG51+UYY8DpwK7c7ATFwTfh4Bo4ddJ5TWL8XxKbOQUWXhPeOsPMgiNCpJUsgQqo27uWLAsO\nY7zR2ggV65zBdQffhIr1PX0TOTNhzkdh6lIoPg8q3gr71UyRwoIjQpTOmE3VnzNo27eOrAtv8boc\nY8aHpmonIA68CQdXQ+UWZ8yExMCk+VD+OSckis+DtLze753gLlUUxquZIoUFR4SYU5jJm1rKmcff\n9roUY8YmVWcy0e6QOPAm1Ox2XotLgsnl8L6vOCExZQkkpg/9mWG+milSWHBEiKT4WI6klnFhy+PQ\n2hDYl9YYMzDfjuzuVkWDO2tRUiZMORcWXQfFS50R2HZRSsAsOCJIx6SFxOx/jK7Dm4iZ/n6vyzEm\nsvUdfHfxtyBnhv+O7PQCpyXR3T8xsSzqL4n1kgVHBMmccS7sh9rda8mx4DBmYFuehBW39XRc1x2C\n33+h5/W+HdnZ00Y2wM70YsERQc4oncahl/KQA7Y2hzH9dHXB0bdh7+uw6u6e6Tp8peTAF9f278g2\no8qCI4LMyEtjKxmceeQVZyWvcXSVhjF+1R6Afatg3+uw7w1oOTH4/s0nLDTCwIIjgqx66hdcGPMe\nsbiDiuoO0fncbby+4xiXfeZL3hZnTDi0nIT3/uK0Kva9Dif2OdvTC5xZYksvhpIL4YFLwz6VuOlh\nwRFBLjh4L/F09toW29nCBQfvBSw4zBjU0QaH1/cExeENzmjs+FSYdgEsWQ7TL4a82b37KDyYStz0\nsOCIIElNlUFtNybqqDrrzux73QmLA3+DtkZnwN3ks+B9X3VaFZPLIS5h4M/xYCpx08OCI5JkFlnz\n24w9jcedforuVkWD+4fQhOkw/zNOUEx7HyRnBfe543TwXSSw4Igkl95J53O3EdvZZ6rmzCnQ2Q6x\n8d7UZUww2pqdsRTdrYrj25ztyRNg+oXOqafpF0H2VC+rNCNgwRFBVqdewvNdn+fO1KdIbKrkiOaw\npauEDx1cDY9/Bj79sI0oN97rO/Duku8404d3tyoOrYXONohNgOJz4dLvOq2KSQts0N0YIapRtcpq\nQMrLy3X9+uhbLPC+N/YyvyiTpaW5ALx96CTXP7iWZV2v8v3YB5BJ8+Da30F6vseVmnGr76p3feWf\nCaUXOa2K4vMgISWs5Y1HfX9vAKzeW82WijpuvrA04M8RkQ2qGtDU3Bb/EeTmC0t7/cdfMCWL399y\nPq+nXMEXOr9G5/Fd8OBlUL3bwyrNuNTRCntfg+f/2X9oJE+Ar+6GL/wVPvBvMONSC40wmV+Uya2P\nb+L1d49R3djK6r3V3Pr4JuYXZYbsmNbiiAJH605xw0PrSK3ZwhOp/0liTBdc8xvnNIAxoVJ7APa8\nArtXwv43oL15kJ0F7joZttIiwWj9pa+qtHZ0UX+qnfqWDhpOtdNwqsO9OY/r+9w3dO/b6m5vaadL\nITUhlsT4WH5+7aJedQUimBaH9XFEgUmZSTz5T+fxj4/EcfmB7/B89n+R8cgy+MQDznw8xoyGjlan\nU3vPStj9ClTvdLZnTYWF18HMy50WR/3h/u8dh1f+df+l/7OrFzGnIJ03dlVx14pt3HrxDP74TuXp\nX/T1p3zDoCcUfEOgvXPwP+BFID0xjvSkeNKT4shIiqcgM4lZSWlkJDvbNh08yeq9Ndx0QUnQoREs\na3FEkVPtndz6+CY27NjNC3m/oKBhK/KhH8E5y70uzUSrkwedkNiz0pnSo73J6dSeej7M/IATFjkz\negbf+evjiE+Gj/5szF0a29rRSVVDK8cbWjle30pVw6nTj4+7jytqW6hraR/ys9IS407/wk9PinNv\nbggkx59+npHku1/PvqkJccTEDDxJY/fpqevPKebRtQetxWF6JMXHct/1i/nWs/Fcsv4rPDPxIcr+\n+DWor4BL77IrVszQOtqcKcd3v+yERdW7zvasYlhwtRMUJe+HhFT/74+QgXcjOU3U1NrhBoAbBA1O\nEFTV9zw+3tDKyeb+gRAjkJOWyMR05zavMJN91Y289V4tHyjL59PlU3pCINm5T0uMI3aQX/oj1R0a\n3WFxbmlOr+ehYC2OKKSq/Oilnfxy1W7+d+KTXFi/As78FCz7hS1GY/o7eah3X0VbI8TEw7TzYcbl\nTljkzoqqacf7/rJcvaeaLz6+kW9+aA6FWUlOC6Gxd+ugyg2LprbOfp+XEBtDXnoieW4gTMxIZGJ6\nUr/HE1ITiIuN6VfHSP7SHykvrqryJDhE5FPAXcAcYImq+v0tLyJXAD8FYoEHVPXuQD5/rAdHtwf/\nup/vP7+NH+St5NqG/3VG3179mLO6mRm/Otrg0BrnFNTuV5xV8MAZSDrzcicsSt4PiWne1hmkupZ2\n3qtu4r2aJvZXN/HW/hOs3X+C5PhYGlo7/L4nNSGWiRlJPYGQnuQGQc/jvLREslLikSCDs1949Xke\nbaLhVNVW4OPALwfaQURigV8AlwMVwFsiskJVt4enxMh30wUl5KYl8C9PxnAkK5t/Ofgz5KEPwXW/\ng8zJXpdnRlvfgXe+p4jqDrutilecvoq2BqdVMfU8Z3nUGZf3nygwAjWcaue96mb21zQ5IeEGxXs1\nzZxoaju9nwgUZiYzKTOJitoWFhdn8eH5hadPIU3McFoIqYmh+xW3paKuV0gsLc3l59cuYktFXVQG\nRzA8PVVmBEDhAAATHElEQVQlIquAr/prcYjIecBdqvpB9/k3AVT1h0N97nhpcXRbtfM4X3h0I1ek\nvMtPun5MTHImXP80TJzjdWlmtPjrlI5NdKbuqKvomdYjowhmXuZ0bJe8PyJnGmhs7egJhOom9lc3\nc6DGeV7d2NZr34LMJKblpDItN5WS3JTTj4snpLDxYK3np4nGkmhocQRiMuA7418FcM5AO4vIcmA5\nQHFxcWgrizAXzZ7I458/hxt/HcP18l0e7vgx8Q9+EK553Jma2kS/V7/Xf+BdZyvsfsk5RXn595yw\nyDsjpK2KQM+nN7d18F518+nTSgdqmk63JKoaeq/cl5+RyLScVC6bk8/UHDcgclOZOiGV5IRYv3V4\n0SFseoSsxSEiK4FJfl76tqo+5+6zioFbHJ8ErlDVf3SffxY4R1VvHerY463F0W3P8UZueGgdyS2V\nrMj6T1IaD8LH7oN5n/C6NDMcrY3OtON7X4O19w2wU3gH3vn+wl5cnM1zmw/z/ed38NEFBajCfrcl\ncay+dzjkpSdSkpPK1JwUt/WQ6rYeUkhJCP7v19HqEDY9IqLFoaqXjfAjDgNTfJ4XudvMAGZMTOOp\nL5zHDQ+t4/3Vd/BywX1MeOpGqK+EpUPmrfFaVxcc3eIExd7X4OAa6GqHuGSIS4KOU/3fE4aBd11d\nysETzeyorGd7ZT3TclK4/oG1dPn8zfnEukPkpiUwNSeVC2bknW41dJ9aShvlvgZ/4bC0NNdaG2ES\nyaeq3gJmikgJTmBcDVzrbUmRryAzmSf/6Txueng95x28lZenPs7Ul7/tjPb9wL/bWI9IU1/pTj/+\nmjOzbHO1sz3/TDj3C1B6iTNZ4I4VYVnxrqWtk53HGpyQOFLPjkrn1n0Ja2yMMD03lRkT09h1rJEP\nzp3ErRfPYGpuChlJNu3/eOFJcIjIx4B7gDzgBRHZrKofFJFCnMtur1TVDhG5FXgJ53Lch1R1mxf1\nRpuslAQevekcbn18Ixe9+/c8PT2HxWvudcLjY/dDfJLXJY5f7S3OtB7drYrj7kWCqXnOxICllzgd\n3ul9zvKGYODd8YZTbjg0sL2ynu1H6thf3XS6JZGWGMecgnQ+eVYRZYUZzCnIYFZ++ulO6dsumcGj\naw/S0NpuoTHO2ADAMay9s4s7nn6HpzdWcO/0N7nyyD1QvNQZ65EywevyxgdVJxy6g+LAaueUU2yC\n05IovcS55c8LWWuwo7OL/dVNbjg4p5t2VNb3uoJpclby6XAoc29F2cn9prkYa2MXTI+IHwAYahYc\nPVSVu//0Lr98Yx/fmfYuN1XdjWSXwPVPOdNMmNHXWOUuauSGReNRZ3vu7J5WxdSlA0/r4UegncH1\np9p5t9LnVNPRenYebaC1owtwRkjPzE+jrMANicIM5kzKIDMlsBaDdUqPXRYcFhz9PPCXffzbCzu4\nsegw/9rwb0h8sjNQsGC+16VFv+6R2t1BUfm2sz0521nQqPQSZwW8EXRk+51i47GN3HhBCaqwvbKO\nHZUNHDzRM/V5dko8ZYUZvUKiNC+N+Fjr5zL9WXBYcPj1zMYKvv7UFi7PPcHP9QfEttbDZx5xfrGZ\ngfkbsV24CPa86gTFe391ZpWNiYOiJTDDPf1UsBBi/I9DCEZnl7K/upFnNh7mob/uJyctgSMnT9H9\nf64IlOSkng6H7qDIz0gMehoNM35ZcFhwDOj1ncf54qMbKUtr5ImUn5BQu8uZHHHB1V6XFpn8LpUq\n0P1re8J0KHVPP027AJIyRnS47pB453AdWyrq2Hq4jm1H6mk+fVUTdHbBvMkZXH12MWWFGczOTw/p\n1BpmfIiIcRwmMl08eyKPff4cbvz1W3yg/Q6eL/glac/+k3PF1QVfifi5jMKiqRqObIbKTfDnn0BH\n36VS1TkN9fnXYULJsA/THRJbKup453D/kEiKj6GsIINPl09h3uRMOruUu/+4g8+eO5VH1x5kel4q\ni4uzR/APNWZ4LDjGocXF2Tx183l89sF1vK/iFl6enkveq99zJsq78sejcnolajSfgCOboHKzc3/k\nbag7OPT7Wk4GFRqdXcq+KqclEUhInDk5k9K81NNTeHf3cfziusU2xYbxnJ2qGseOnGzh7x9ax6ET\njfxx7utM3/krmP1hZ0nahBSvyxt9LbVOS8I3KE76hMSE6U6/ROEiKFwIBQvgf86HukP9PytzCvzz\nVr+HCTQk5hdl+Q0Jf+xqJhNq1sdhwRGw2qY2bnz4Ld4+dJLfLtrK2dt/CNnTnAn06is9W+FtxFpO\nuuHgExS17/W8nj3NJyQWOSGRnNXvY1b+9h4u3vXvxHb2nK7qjE3m9Vnf5rLPfCkkIWGMFyw4LDiC\n0tzWwRcf28iqnVU8M/0PLDryBL16OrxaU3qw9Sd8napzLoH1DYkT+3pezyp2w2FhT0gEOABy9d5q\nnn/0p9yZ8hSJTZU0J0/iruZP0DDjY9Q0tVlImDHDgsOCI2jtnV1846ktfGXbJyiKqe73ekdMAnGz\nPgAJac5prPgUZwBb973v437bUiA+FeISAi/I39VM8clwxX9ATqnbH+EGxYm9PftkFkPhgt5BEeQo\n+abWDvZXN52+rdt/gjf31SBAhzsfh4WEGWssOCw4hqWrS5HvZSP0/04oIBPLoK0J2pt77oMRE98T\nIgndIZPaO4i6w2bDr6G1fvDPyyhy+iIKF/YERWpgHcVtHV0cqm1mf5UTDvuqm9hf3cj+6v5Tghdm\nJhETI1TUtnDZnHy+9sHZFhJmzLHLcc2wxMSIc0rIT2fwqZRCtl75AhNSE8hJTSAjKZ4Y1LlUta0Z\n2hrdQGl2BsO1Nfk87g4af9uanSub2ivc97jb/E0h3u26p5yQSMsb9N/T1aUcazjF/qom9lY3uSHh\nhMOh2hY6feYFz06JpyTXmRJ8ep6zXkT3mhGbDvWe1K+mqZXZkyJvZT1jwsWCw/R26Z39ThE1awJ3\nnLyKFfe9eXpbbIyQnRJPdkqCEyZpCWSnJJCTmkF2aq4bMIlk58Y796nxJMYFfplvww9nk956tP/2\nxALSZ17ea9vJ5janxVDVc3ppn7tedUt75+n9kuJjKMlNY25hJh+ZX+iEQ14qJTmpZKf6P41mK80Z\n05+dqjL97HrlQdL+9gMKqKGSHCrLv07ioqupaWqjtqmNmqY2TjS1cqKpnRNNrdQ2tVPT1Eptczu1\nzW0M9JVKS4wjOzWeCamJ5KS6QXM6cBLITk043aJpeOtxZqz9Fsn0zODaQgKvz/wO+ws/zD6f1kNt\nc/vpfWJjhCnZyW6LIY2SvFSm56YyPS+V/PSkfrO9DsUugzXjhfVxWHAM20inze7sUk42t1Hb3EZN\no3vf1MaJxjZONLdxoqnn1h1C3TO39nVV7F/5atyTFFLDEc3hRx2fZkWXs4Z6fkbi6XCY7p5WKslL\nZUp2Cglx1vdgTLAsOCw4hi3cf2GrKs1tnT2B0uyETHfgvLGziu2V9ZxTMoHrz53q9DuEYClSY8Y7\n6xw3wxbutZxFhNTEOFIT45gyofdo9dV7q/ntW4dOd0rnpCUwb3JmSOowxgTO2vQmIvmeIvvKB2bz\n82sXcevjm1i9t/8YE2NMeFlwmIi0paKuV7/K0tJcfn7tIrZU1HlcmTHG+jiMMcYE1cdhLQ5jjDFB\n8SQ4RORTIrJNRLpExG/CicgUEXldRLa7+94e7jqNMcb051WLYyvwceDPg+zTAfyLqpYB5wK3iEhZ\nOIozxhgzME8ux1XVHeBcijnIPpVApfu4QUR2AJOB7eGo0RhjjH9R0cchItOARcBabysxxhgTshaH\niKwEJvl56duq+lwQn5MGPA18WVUHnGdbRJYDy92njSKyM5h6g5ALRMtggmip1eocXdFSJ0RPreOh\nzqmB7ujp5bgisgr4qqr6vXZWROKB54GXVPU/w1nbQERkfaCXrHktWmq1OkdXtNQJ0VOr1dlbxJ6q\nEqcD5EFgR6SEhjHGGO8ux/2YiFQA5wEviMhL7vZCEXnR3e184LPAJSKy2b1d6UW9xhhjenh1VdWz\nwLN+th8BrnQf/xUIbvGE8Ljf6wKCEC21Wp2jK1rqhOip1er0MSanHDHGGBM6EdvHYYwxJjJZcPgQ\nkStEZKeI7BGRO/y8fp2IbBGRd0RktYgs8HntPXf7ZhEJ6QyLAdR5kYjU+fQN3Rnoe8Nc59d8atwq\nIp0iMsF9LZw/z4dE5LiIbB3gdRGRn7n/ji0istjntXD+PIeqMyK+nwHWGinf0aHqjJTv6JBTMIX1\ne6qqdnNO18UCe4HpQALwNlDWZ5+lQLb7+EPAWp/X3gNyI6TOi4Dnh/PecNbZZ/+PAq+F++fpHuv9\nwGJg6wCvXwn8EafP7dzu/+7h/HkGWKfn388gavX8OxpInRH0HS0AFruP04Fdfv6/D9v31FocPZYA\ne1R1n6q2Ab8BlvnuoKqrVbXWfboGKApzjRBAnSF6b6jrvAZ4IkS1DEpV/wycGGSXZcAj6lgDZIlI\nAeH9eQ5ZZ4R8P7trGepnOpCI+pn24eV3tFJVN7qPG4DuKZh8he17asHRYzJwyOd5Bf3/w/i6CSfd\nuymwUkQ2iDOKPVQCrXOp21z9o4jMDfK9oyHgY4lICnAFzgwB3cL18wzEQP+WcP48g+XV9zMYXn9H\nAxZJ31EZeAqmsH1Pbc3xYRCRi3H+x7zAZ/MFqnpYRCYCr4jIu+5fM17YCBSraqM4Y19+D8z0qJZA\nfBT4m6r6/uUXST/PqBIF30+w7+iwSIBTMIWatTh6HAam+Dwvcrf1IiLzgQeAZapa071dVQ+798dx\nxqgs8apOVa1X1Ub38YtAvIjkBvLecNbp42r6nAII488zEAP9W8L58wxIBHw/AxIh39FgeP4dFWcK\npqeBx1T1GT+7hO97Go6OnWi44bS+9gEl9HQgze2zTzGwB1jaZ3sqkO7zeDVwhYd1TqJnjM4S4CBO\nh9mQ7w1nne5+mTjnmFO9+Hn6HHMaA3fkfpjenY7rgvk3hrFOz7+fQdTq+Xc0kDoj5Tvq/mweAf57\nkH3C9j21U1UuVe0QkVuBl3CuQnhIVbeJyM3u6/cBdwI5wL3irCXSoc6EYvnAs+62OOBxVf2Th3V+\nEviCiHQALcDV6nyD/L7XwzoBPga8rKpNPm8P288TQESewLnKJ1ecqXC+C8T71PkizhUre4Bm4HOD\n/Rs9rNPz72cQtXr+HQ2wToiA7yg9UzC9IyKb3W3fwvljIezfUxs5bowxJijWx2GMMSYoFhzGGGOC\nYsFhjDEmKBYcxhhjgmLBYYwxJigWHCZquDOTbhaRt0Vko4gsHaXPvUhEng90+ygc7yoRKfN5vkpE\nBl0nWpzVMZ8K8jir3BlRN4vIjkCmxRCRxmCOYcYnCw4TTVpUdaGqLgC+CfzQ64KG6SqgbMi9fKjq\nEVX95DCOdZ2qLsQZB/AfIpIwjM/oRURs/Nc4Z8FholUGUAun1yH4sbtewjsi8hl3+0XuX91Pici7\nIvKYuCO23PUJ3hWRjcDHhzqYiKSKs3bDOhHZJCLL3O3/ICLPiMifRGS3iPzI5z03icgu9z2/EpGf\nu62kvwN+7LYESt3dP+Xut0tE3ufn+NPEXTNisGMOIg1oAjrdz7jG/VltFZH/6HOsf3dbdWtEJN/d\n9msRuU9E1gKBHM+MYfaXg4kmye6o2SSc9Qkucbd/HFgILABygbdEpHuyuUXAXOAI8DfgfHEW3fmV\n+/49wG8DOPa3cdZiuFFEsoB1IrLSfW2he5xWYKeI3IPzC/pfcdZ6aABeA95W1dUisgJnLYqnALpH\nH6vqEnfCv+8Clw1RT79jquohP/s9JiKtOBMIfllVO0WkEPgP4Cyc8H1ZRK5S1d/jTJ+xRlW/7QbS\n54F/cz+rCGc6k84Afl5mDLMWh4km3aeqzsCZ4voRtwVxAfCEqnaq6jHgDeBs9z3rVLVCVbuAzTjz\nEp0B7FfV3e40F48GcOwPAHe4wbUKJ7yK3ddeVdU6VT0FbAem4sy/9IaqnlDVduB3Q3x+96R1G9wa\nh+LvmP5cp6rz3Vq/KiJTcX42q1S1SlU7gMdwFjQCaAO6+3X61vI7Cw0D1uIwUUpV33RnU80bYtdW\nn8edDP87L8AnVHVnr40i54zSMbo/I9D3B3VMVa1yT8v1rbevdu2Zh6jv5zb52d+MQ9biMFFJRM7A\nmbCtBvgL8BkRiRWRPJy/ntcN8vZ3gWk+/QvXBHDIl4Av+fSRLBpi/7eAC0Uk2+1M/oTPaw04y3+G\njTgLES3CWUJ0nVtbrojE4vz73whnPSa6WYvDRJPuPg5wWgA3uOfsnwXOw5kuWoGvq+pRN1z6UdVT\n7qWpL4hIM07wDPWL/PvAfwNbRCQG2A98ZKCd1Vng5wc4v6RP4IRVnfvyb4BfichtOLPEhtJjItIC\nJAK/VtUNACJyB/A6zs/xBVV9LsR1mDHEZsc1JkREJE2dFe7icBb6eUhVn/W6LmNGyk5VGRM6d7kt\npK04LZTfe1yPMaPCWhzGGGOCYi0OY4wxQbHgMMYYExQLDmOMMUGx4DDGGBMUCw5jjDFBseAwxhgT\nlP8PNn8Mh+m1n+QAAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from fermilib.utils import MolecularData\n", + "from fermilibpluginpsi4.run_psi4 import run_psi4\n", + "\n", + "# Set molecule parameters.\n", + "basis = 'sto-3g'\n", + "multiplicity = 1\n", + "bond_length_interval = 0.2\n", + "n_points = 10\n", + "\n", + "# Set calculation parameters.\n", + "run_scf = 1\n", + "run_mp2 = 1\n", + "run_cisd = 0\n", + "run_ccsd = 0\n", + "run_fci = 1\n", + "delete_input = True\n", + "delete_output = True\n", + "\n", + "# Generate molecule at different bond lengths.\n", + "hf_energies = []\n", + "fci_energies = []\n", + "bond_lengths = []\n", + "for point in range(1, n_points + 1):\n", + " bond_length = bond_length_interval * float(point)\n", + " bond_lengths += [bond_length]\n", + " geometry = [('H', (0., 0., 0.)), ('H', (0., 0., bond_length))]\n", + " molecule = MolecularData(\n", + " geometry, basis, multiplicity,\n", + " description=str(round(bond_length, 2)))\n", + " \n", + " # Run Psi4.\n", + " molecule = run_psi4(molecule,\n", + " run_scf=run_scf,\n", + " run_mp2=run_mp2,\n", + " run_cisd=run_cisd,\n", + " run_ccsd=run_ccsd,\n", + " run_fci=run_fci)\n", + "\n", + " # Print out some results of calculation.\n", + " print('\\nAt bond length of {} Bohr, molecular hydrogen has:'.format(\n", + " bond_length))\n", + " print('Hartree-Fock energy of {} Hartree.'.format(molecule.hf_energy))\n", + " print('MP2 energy of {} Hartree.'.format(molecule.mp2_energy))\n", + " print('FCI energy of {} Hartree.'.format(molecule.fci_energy))\n", + " print('Nuclear repulsion energy between protons is {} Hartree.'.format(\n", + " molecule.nuclear_repulsion))\n", + " for orbital in range(molecule.n_orbitals):\n", + " print('Spatial orbital {} has energy of {} Hartree.'.format(\n", + " orbital, molecule.orbital_energies[orbital]))\n", + " hf_energies += [molecule.hf_energy]\n", + " fci_energies += [molecule.fci_energy]\n", + "\n", + "# Plot.\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "\n", + "plt.figure(0)\n", + "plt.plot(bond_lengths, fci_energies, 'x-')\n", + "plt.plot(bond_lengths, hf_energies, 'o-')\n", + "plt.ylabel('Energy in Hartree')\n", + "plt.xlabel('Bond length in Bohr')\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/fermilibpluginpsi4/__init__.py b/fermilibpluginpsi4/__init__.py new file mode 100644 index 0000000..2403ef8 --- /dev/null +++ b/fermilibpluginpsi4/__init__.py @@ -0,0 +1,21 @@ +# FermiLib plugin to interface with Psi4 +# Copyright (C) 2017 FermiLib Developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +""" +FermiLib plugin to interface with Psi4 +""" + +from ._version import __version__ diff --git a/fermilibpluginpsi4/_psi4_conversion_functions.py b/fermilibpluginpsi4/_psi4_conversion_functions.py new file mode 100755 index 0000000..9059e30 --- /dev/null +++ b/fermilibpluginpsi4/_psi4_conversion_functions.py @@ -0,0 +1,260 @@ +# FermiLib plugin to interface with Psi4 +# Copyright (C) 2017 FermiLib Developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Helper functions for parsing data files of different types.""" +from __future__ import absolute_import + +import numpy + +from fermilib.ops import InteractionOperator + + +def unpack_spatial_rdm(one_rdm_a, + one_rdm_b, + two_rdm_aa, + two_rdm_ab, + two_rdm_bb): + """ + Covert from spin compact spatial format to spin-orbital format for RDM. + + Note: the compact 2-RDM is stored as follows where A/B are spin up/down: + RDM[pqrs] = <| a_{p, A}^\dagger a_{r, A}^\dagger a_{q, A} a_{s, A} |> + for 'AA'/'BB' spins. + RDM[pqrs] = <| a_{p, A}^\dagger a_{r, B}^\dagger a_{q, B} a_{s, A} |> + for 'AB' spins. + + Args: + one_rdm_a: 2-index numpy array storing alpha spin + sector of 1-electron reduced density matrix. + one_rdm_b: 2-index numpy array storing beta spin + sector of 1-electron reduced density matrix. + two_rdm_aa: 4-index numpy array storing alpha-alpha spin + sector of 2-electron reduced density matrix. + two_rdm_ab: 4-index numpy array storing alpha-beta spin + sector of 2-electron reduced density matrix. + two_rdm_bb: 4-index numpy array storing beta-beta spin + sector of 2-electron reduced density matrix. + + Returns: + one_rdm: 2-index numpy array storing 1-electron density matrix + in full spin-orbital space. + two_rdm: 4-index numpy array storing 2-electron density matrix + in full spin-orbital space. + """ + # Initialize RDMs. + n_orbitals = one_rdm_a.shape[0] + n_qubits = 2 * n_orbitals + one_rdm = numpy.zeros((n_qubits, n_qubits)) + two_rdm = numpy.zeros((n_qubits, n_qubits, + n_qubits, n_qubits)) + + # Unpack compact representation. + for p in range(n_orbitals): + for q in range(n_orbitals): + + # Populate 1-RDM. + one_rdm[2 * p, 2 * q] = one_rdm_a[p, q] + one_rdm[2 * p + 1, 2 * q + 1] = one_rdm_b[p, q] + + # Continue looping to unpack 2-RDM. + for r in range(n_orbitals): + for s in range(n_orbitals): + + # Handle case of same spin. + two_rdm[2 * p, 2 * q, 2 * r, 2 * s] = ( + two_rdm_aa[p, r, q, s]) + two_rdm[2 * p + 1, 2 * q + 1, 2 * r + 1, 2 * s + 1] = ( + two_rdm_bb[p, r, q, s]) + + # Handle case of mixed spin. + two_rdm[2 * p, 2 * q + 1, 2 * r, 2 * s + 1] = ( + two_rdm_ab[p, r, q, s]) + two_rdm[2 * p, 2 * q + 1, 2 * r + 1, 2 * s] = ( + -1. * two_rdm_ab[p, s, q, r]) + two_rdm[2 * p + 1, 2 * q, 2 * r + 1, 2 * s] = ( + two_rdm_ab[q, s, p, r]) + two_rdm[2 * p + 1, 2 * q, 2 * r, 2 * s + 1] = ( + -1. * two_rdm_ab[q, r, p, s]) + + # Map to physicist notation and return. + two_rdm = numpy.einsum('pqsr', two_rdm) + return one_rdm, two_rdm + + +def parse_psi4_ccsd_amplitudes(number_orbitals, + n_alpha_electrons, n_beta_electrons, + psi_filename): + """Parse coupled cluster singles and doubles amplitudes from psi4 file. + + Args: + number_orbitals(int): Number of total spin orbitals in the system + n_alpha_electrons(int): Number of alpha electrons in the system + n_beta_electrons(int): Number of beta electrons in the system + psi_filename(str): Filename of psi4 output file + + Returns: + molecule(InteractionOperator): Molecular Operator instance holding ccsd + amplitudes + + """ + output_buffer = [line for line in open(psi_filename)] + + T1IA_index = None + T1ia_index = None + T2IJAB_index = None + T2ijab_index = None + T2IjAb_index = None + + # Find Start Indices + for i, line in enumerate(output_buffer): + if ('Largest TIA Amplitudes:' in line): + T1IA_index = i + + elif ('Largest Tia Amplitudes:' in line): + T1ia_index = i + + elif ('Largest TIJAB Amplitudes:' in line): + T2IJAB_index = i + + elif ('Largest Tijab Amplitudes:' in line): + T2ijab_index = i + + elif ('Largest TIjAb Amplitudes:' in line): + T2IjAb_index = i + + T1IA_Amps = [] + T1ia_Amps = [] + + T2IJAB_Amps = [] + T2ijab_Amps = [] + T2IjAb_Amps = [] + + # Read T1's + if (T1IA_index is not None): + for line in output_buffer[T1IA_index + 1:]: + ivals = line.split() + if not ivals: + break + T1IA_Amps.append([int(ivals[0]), int(ivals[1]), float(ivals[2])]) + + if (T1ia_index is not None): + for line in output_buffer[T1ia_index + 1:]: + ivals = line.split() + if not ivals: + break + T1ia_Amps.append([int(ivals[0]), int(ivals[1]), float(ivals[2])]) + + # Read T2's + if (T2IJAB_index is not None): + for line in output_buffer[T2IJAB_index + 1:]: + ivals = line.split() + if not ivals: + break + T2IJAB_Amps.append([int(ivals[0]), int(ivals[1]), + int(ivals[2]), int(ivals[3]), + float(ivals[4])]) + + if (T2ijab_index is not None): + for line in output_buffer[T2ijab_index + 1:]: + ivals = line.split() + if not ivals: + break + T2ijab_Amps.append([int(ivals[0]), int(ivals[1]), + int(ivals[2]), int(ivals[3]), + float(ivals[4])]) + + if (T2IjAb_index is not None): + for line in output_buffer[T2IjAb_index + 1:]: + ivals = line.split() + if not ivals: + break + T2IjAb_Amps.append([int(ivals[0]), int(ivals[1]), + int(ivals[2]), int(ivals[3]), + float(ivals[4])]) + + # Determine if calculation is restricted / closed shell or otherwise + restricted = T1ia_index is None and T2ijab_index is None + + # Store amplitudes with spin-orbital indexing, including appropriate + # symmetry + single_amplitudes = numpy.zeros((number_orbitals, ) * 2) + double_amplitudes = numpy.zeros((number_orbitals, ) * 4) + + # Define local helper routines for clear indexing of orbitals + def alpha_occupied(i): + return 2 * i + + def alpha_unoccupied(i): + return 2 * (i + n_alpha_electrons) + + def beta_occupied(i): + return 2 * i + 1 + + def beta_unoccupied(i): + return 2 * (i + n_beta_electrons) + 1 + + # Store singles + for entry in T1IA_Amps: + i, a, value = entry + single_amplitudes[alpha_occupied(i), + alpha_unoccupied(a)] = value + if (restricted): + single_amplitudes[beta_occupied(i), + beta_unoccupied(a)] = value + + for entry in T1ia_Amps: + i, a, value = entry + single_amplitudes[beta_occupied(i), + beta_unoccupied(a)] = value + + # Store doubles, include factor of 1/2 for convention + for entry in T2IJAB_Amps: + i, j, a, b, value = entry + double_amplitudes[alpha_occupied(i), + alpha_unoccupied(a), + alpha_occupied(j), + alpha_unoccupied(b)] = -value / 2. + if (restricted): + double_amplitudes[beta_occupied(i), + beta_unoccupied(a), + beta_occupied(j), + beta_unoccupied(b)] = -value / 2. + + for entry in T2ijab_Amps: + i, j, a, b, value = entry + double_amplitudes[beta_occupied(i), + beta_unoccupied(a), + beta_occupied(j), + beta_unoccupied(b)] = -value / 2. + + for entry in T2IjAb_Amps: + i, j, a, b, value = entry + double_amplitudes[alpha_occupied(i), + alpha_unoccupied(a), + beta_occupied(j), + beta_unoccupied(b)] = -value / 2. + + if (restricted): + double_amplitudes[beta_occupied(i), + beta_unoccupied(a), + alpha_occupied(j), + alpha_unoccupied(b)] = -value / 2. + + # Package into InteractionOperator. + molecule = InteractionOperator(0.0, + single_amplitudes, + double_amplitudes) + return molecule diff --git a/fermilibpluginpsi4/_psi4_template b/fermilibpluginpsi4/_psi4_template new file mode 100755 index 0000000..55cb60a --- /dev/null +++ b/fermilibpluginpsi4/_psi4_template @@ -0,0 +1,235 @@ +"""This is a template for psi4 input format.""" +import numpy +import sys + +from fermilib.config import * +from fermilib.ops._interaction_tensor import (one_body_basis_change, + two_body_basis_change) +from fermilib.utils import MolecularData + +sys.path.append('&THIS_DIRECTORY') +from _psi4_conversion_functions import * + + +# Set memory that job can use in megabytes. +memory &memory mb + +# Initialize molecular data. +_description = '&description' +if _description == 'None': + _description = None +molecule = MolecularData(&geometry, + '&basis', + &multiplicity, + &charge, + _description) + +# Set molecular geometry and symmetry. +molecule mol { +&geo_string +symmetry c1 +} +mol.set_multiplicity(&multiplicity) +mol.set_molecular_charge(&charge) + +# Set reference and guess. +if molecule.multiplicity == 1: + set reference rhf + set guess sad +else: + set reference rohf + set guess gwh + +# Set global parameters of calculation. +set globals { + basis &basis + freeze_core false + fail_on_maxiter true + df_scf_guess true + opdm true + tpdm true + soscf false + scf_type df + maxiter 1e6 + num_amps_print 1e6 + r_convergence 1e-6 + d_convergence 1e-6 + e_convergence 1e-6 + ints_tolerance EQUALITY_TOLERANCE + damping_percentage 0 +} + +# Run self-consistent field (SCF) calculation. +if &run_scf: + try: + hf_energy, hf_wavefunction = energy('scf', return_wfn=True) + if &verbose: + print('Hartree-Fock energy for {} ({} electrons) is {}.'.format( + molecule.name, molecule.n_electrons, hf_energy)) + except: + if &tolerate_error: + print('WARNING: SCF calculation failed.') + else: + raise + else: + # Get orbitals and Fock matrix. + molecule.hf_energy = hf_energy + molecule.nuclear_repulsion = mol.nuclear_repulsion_energy() + molecule.canonical_orbitals = numpy.asarray(hf_wavefunction.Ca()) + molecule.n_orbitals = molecule.canonical_orbitals.shape[0] + molecule.n_qubits = 2 * molecule.n_orbitals + molecule.orbital_energies = numpy.asarray(hf_wavefunction.epsilon_a()) + molecule.fock_matrix = numpy.asarray(hf_wavefunction.Fa()) + + # Get integrals using MintsHelper. + mints = MintsHelper(hf_wavefunction.basisset()) + molecule.one_body_integrals = one_body_basis_change( + numpy.asarray(mints.ao_kinetic()), molecule.canonical_orbitals) + molecule.one_body_integrals += one_body_basis_change( + numpy.asarray(mints.ao_potential()), molecule.canonical_orbitals) + two_body_integrals = numpy.asarray(mints.ao_eri()) + two_body_integrals.reshape((molecule.n_orbitals, molecule.n_orbitals, + molecule.n_orbitals, molecule.n_orbitals)) + two_body_integrals = numpy.einsum('psqr', two_body_integrals) + two_body_integrals = two_body_basis_change( + two_body_integrals, molecule.canonical_orbitals) + integrals_name = molecule.filename + '_eri' + numpy.save(integrals_name, two_body_integrals) + molecule.save() + + +# Perform MP2 energy calculation if there are at least two electrons. +if &run_mp2: + try: + assert molecule.n_electrons > 1 + mp2_energy = energy('mp2') + if &verbose: + print('MP2 energy for {} ({} electrons) is {}.'.format( + molecule.name, molecule.n_electrons, mp2_energy)) + except: + if &tolerate_error: + print('WARNING: MP2 calculation failed.') + else: + raise + else: + molecule.mp2_energy = mp2_energy + molecule.save() + + +# Perform configuration interaction singles and doubles (CISD) calculation. +if &run_cisd: + set qc_module detci + try: + cisd_energy, cisd_wavefunction = energy('cisd', return_wfn=True) + if &verbose: + print('CISD energy for {} ({} electrons) is {}.'.format( + molecule.name, molecule.n_electrons, cisd_energy)) + except: + if &tolerate_error: + print('WARNING: CISD calculation failed.') + else: + raise + else: + # For the functions below, "a" and "b" refer to "up and "down" spins. + molecule.cisd_energy = cisd_energy + + # Get 1-RDM from CISD calculation. + cisd_one_rdm_a = numpy.array(cisd_wavefunction.get_opdm( + 0, 0, 'A', True)).reshape(molecule.n_orbitals, molecule.n_orbitals) + cisd_one_rdm_b = numpy.array(cisd_wavefunction.get_opdm( + 0, 0, 'B', True)).reshape(molecule.n_orbitals, molecule.n_orbitals) + + # Get 2-RDM from CISD calculation. + cisd_two_rdm_aa = numpy.array(cisd_wavefunction.get_tpdm( + 'AA', False)).reshape(molecule.n_orbitals, molecule.n_orbitals, + molecule.n_orbitals, molecule.n_orbitals) + cisd_two_rdm_ab = numpy.array(cisd_wavefunction.get_tpdm( + 'AB', False)).reshape(molecule.n_orbitals, molecule.n_orbitals, + molecule.n_orbitals, molecule.n_orbitals) + cisd_two_rdm_bb = numpy.array(cisd_wavefunction.get_tpdm( + 'BB', False)).reshape(molecule.n_orbitals, molecule.n_orbitals, + molecule.n_orbitals, molecule.n_orbitals) + + # Get overall RDMs. + cisd_one_rdm, cisd_two_rdm = unpack_spatial_rdm( + cisd_one_rdm_a, cisd_one_rdm_b, cisd_two_rdm_aa, + cisd_two_rdm_ab, cisd_two_rdm_bb) + + # Store 1-RDM in molecule file, 2-RDM separately in other file. + molecule.cisd_one_rdm = cisd_one_rdm + cisd_rdm_name = molecule.filename + '_cisd_rdm' + numpy.save(cisd_rdm_name, cisd_two_rdm) + molecule.save() + + +# Perform exact diagonalization. +if &run_fci: + set qc_module detci + try: + fci_energy, fci_wavefunction = energy('fci', return_wfn=True) + if &verbose: + print('FCI energy for {} ({} electrons) is {}.'.format( + molecule.name, molecule.n_electrons, fci_energy)) + except: + if &tolerate_error: + print('WARNING: FCI calculation failed.') + else: + raise + else: + # For the functions below, "a" and "b" refer to "up and "down" spins. + molecule.fci_energy = fci_energy + + # Get 1-RDM from FCI calculation. + fci_one_rdm_a = numpy.array(fci_wavefunction.get_opdm( + 0, 0, 'A', True)).reshape(molecule.n_orbitals, molecule.n_orbitals) + fci_one_rdm_b = numpy.array(fci_wavefunction.get_opdm( + 0, 0, 'B', True)).reshape(molecule.n_orbitals, molecule.n_orbitals) + + # Get 2-RDM from FCI calculation. + fci_two_rdm_aa = numpy.array(fci_wavefunction.get_tpdm( + 'AA', False)).reshape(molecule.n_orbitals, molecule.n_orbitals, + molecule.n_orbitals, molecule.n_orbitals) + fci_two_rdm_ab = numpy.array(fci_wavefunction.get_tpdm( + 'AB', False)).reshape(molecule.n_orbitals, molecule.n_orbitals, + molecule.n_orbitals, molecule.n_orbitals) + fci_two_rdm_bb = numpy.array(fci_wavefunction.get_tpdm( + 'BB', False)).reshape(molecule.n_orbitals, molecule.n_orbitals, + molecule.n_orbitals, molecule.n_orbitals) + + # Get overall RDMs. + fci_one_rdm, fci_two_rdm = unpack_spatial_rdm( + fci_one_rdm_a, fci_one_rdm_b, + fci_two_rdm_aa, fci_two_rdm_ab, fci_two_rdm_bb) + + # Store 1-RDM in molecule file, 2-RDM separately in other file. + molecule.fci_one_rdm = fci_one_rdm + fci_rdm_name = molecule.filename + '_fci_rdm' + numpy.save(fci_rdm_name, fci_two_rdm) + molecule.save() + + +# Perform coupled cluster singles and doubles (CCSD) calculation. +if &run_ccsd: + set qc_module ccenergy + try: + ccsd_energy = energy('ccsd') + if &verbose: + print('CCSD energy for {} ({} electrons) is {}.'.format( + molecule.name, molecule.n_electrons, ccsd_energy)) + except: + if &tolerate_error: + print('WARNING: CCSD calculation failed.') + else: + raise + else: + molecule.ccsd_energy = ccsd_energy + + # Merge CC amplitudes into molecule by parsing + psi_filename = outfile_name() + ccsd_amplitudes = parse_psi4_ccsd_amplitudes( + 2 * molecule.n_orbitals, + molecule.get_n_alpha_electrons(), + molecule.get_n_beta_electrons(), + psi_filename) + molecule.ccsd_amplitudes = ccsd_amplitudes + molecule.save() diff --git a/fermilibpluginpsi4/_version.py b/fermilibpluginpsi4/_version.py new file mode 100644 index 0000000..1ccb4f7 --- /dev/null +++ b/fermilibpluginpsi4/_version.py @@ -0,0 +1,18 @@ +# FermiLib plugin to interface with Psi4 +# Copyright (C) 2017 FermiLib Developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Define version number here and read it from setup.py automatically""" +__version__ = "0.1a0" diff --git a/fermilibpluginpsi4/run_psi4.py b/fermilibpluginpsi4/run_psi4.py new file mode 100755 index 0000000..a5963a3 --- /dev/null +++ b/fermilibpluginpsi4/run_psi4.py @@ -0,0 +1,210 @@ +# FermiLib plugin to interface with Psi4 +# Copyright (C) 2017 FermiLib Developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Functions to prepare psi4 input and run calculations.""" +from __future__ import absolute_import + +import os +import re +import subprocess + + +def create_geometry_string(geometry): + """This function converts MolecularData geometry to psi4 geometry. + + Args: + geometry: A list of tuples giving the coordinates of each atom. + example is [('H', (0, 0, 0)), ('H', (0, 0, 0.7414))]. Distances in + atomic units. Use atomic symbols to specify atoms. + + Returns: + geo_string: A string giving the geometry for each atom on a line, e.g.: + H 0. 0. 0. + H 0. 0. 0.7414 + """ + geo_string = '' + for item in geometry: + atom = item[0] + coordinates = item[1] + line = '{} {} {} {}'.format(atom, + coordinates[0], + coordinates[1], + coordinates[2]) + if len(geo_string) > 0: + geo_string += '\n' + geo_string += line + return geo_string + + +def generate_psi4_input(molecule, + run_scf, + run_mp2, + run_cisd, + run_ccsd, + run_fci, + verbose, + tolerate_error, + memory): + """This function creates and saves a psi4 input file. + + Args: + molecule: An instance of the MolecularData class. + run_scf: Boolean to run SCF calculation. + run_mp2: Boolean to run MP2 calculation. + run_cisd: Boolean to run CISD calculation. + run_ccsd: Boolean to run CCSD calculation. + run_fci: Boolean to FCI calculation. + verbose: Boolean whether to print calculation results to screen. + tolerate_error: Whether to fail or merely warn when Psi4 fails. + memory: Int giving amount of memory to allocate in MB. + + Returns: + input_file: A string giving the name of the saved input file. + """ + # Create Psi4 geometry string. + geo_string = create_geometry_string(molecule.geometry) + + # Find the psi4_directory. + psi4_directory = os.path.dirname(os.path.realpath(__file__)) + + # Parse input template. + template_file = psi4_directory + '/_psi4_template' + input_template = [] + with open(template_file, 'r') as stream: + for line in stream: + input_template += [line] + + # Populate contents of input file based on automatic parameters. + input_content = [re.sub('&THIS_DIRECTORY', + psi4_directory, line) + for line in input_template] + + # Populate contents of input file based on MolecularData parameters. + input_content = [re.sub('&geometry', str(molecule.geometry), line) + for line in input_content] + input_content = [re.sub('&basis', molecule.basis, line) + for line in input_content] + input_content = [re.sub('&charge', str(molecule.charge), line) + for line in input_content] + input_content = [re.sub('&multiplicity', str(molecule.multiplicity), line) + for line in input_content] + input_content = [re.sub('&description', str(molecule.description), line) + for line in input_content] + input_content = [re.sub('&geo_string', geo_string, line) + for line in input_content] + + # Populate contents of input file based on provided calculation parameters. + input_content = [re.sub('&run_scf', str(run_scf), line) + for line in input_content] + input_content = [re.sub('&run_mp2', str(run_mp2), line) + for line in input_content] + input_content = [re.sub('&run_cisd', str(run_cisd), line) + for line in input_content] + input_content = [re.sub('&run_ccsd', str(run_ccsd), line) + for line in input_content] + input_content = [re.sub('&run_fci', str(run_fci), line) + for line in input_content] + input_content = [re.sub('&tolerate_error', str(tolerate_error), line) + for line in input_content] + input_content = [re.sub('&verbose', str(verbose), line) + for line in input_content] + input_content = [re.sub('&memory', str(memory), line) + for line in input_content] + + # Write input file and return handle. + input_file = molecule.filename + '.inp' + with open(input_file, 'w') as stream: + stream.write(''.join(input_content)) + return input_file + + +def clean_up(molecule, delete_input=True, delete_output=False): + input_file = molecule.filename + '.inp' + output_file = molecule.filename + '.out' + run_directory = os.getcwd() + for local_file in os.listdir(run_directory): + if local_file.endswith('.clean'): + os.remove(run_directory + '/' + local_file) + try: + os.remove('timer.dat') + except: + pass + if delete_input: + os.remove(input_file) + if delete_output: + os.remove(output_file) + + +def run_psi4(molecule, + run_scf=True, + run_mp2=False, + run_cisd=False, + run_ccsd=False, + run_fci=False, + verbose=False, + tolerate_error=False, + delete_input=True, + delete_output=False, + memory=8000): + """This function runs a Psi4 calculation. + + Args: + molecule: An instance of the MolecularData class. + run_scf: Optional boolean to run SCF calculation. + run_mp2: Optional boolean to run MP2 calculation. + run_cisd: Optional boolean to run CISD calculation. + run_ccsd: Optional boolean to run CCSD calculation. + run_fci: Optional boolean to FCI calculation. + verbose: Boolean whether to print calculation results to screen. + tolerate_error: Optional boolean to warn or raise when Psi4 fails. + delete_input: Optional boolean to delete psi4 input file. + delete_output: Optional boolean to delete psi4 output file. + memory: Optional int giving amount of memory to allocate in MB. + + Returns: + molecule: The updated MolecularData object. + + Raises: + psi4 errors: An error from psi4. + """ + # Prepare input. + input_file = generate_psi4_input(molecule, + run_scf, + run_mp2, + run_cisd, + run_ccsd, + run_fci, + verbose, + tolerate_error, + memory) + + # Run psi4. + output_file = molecule.filename + '.out' + try: + process = subprocess.Popen(['psi4', input_file, output_file]) + process.wait() + except: + print('Psi4 calculation for {} has failed.'.format(molecule.name)) + process.kill() + clean_up(molecule, delete_input, delete_output) + if not tolerate_error: + raise + else: + clean_up(molecule, delete_input, delete_output) + + # Return updated molecule instance. + molecule.load() + return molecule diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..aeeebeb --- /dev/null +++ b/requirements.txt @@ -0,0 +1 @@ +fermilib diff --git a/setup.py b/setup.py new file mode 100755 index 0000000..017170a --- /dev/null +++ b/setup.py @@ -0,0 +1,32 @@ +import os + +from setuptools import setup, find_packages + +# This reads the __version__ variable from projectq/_version.py +exec(open('fermilibpluginpsi4/_version.py').read()) +# Readme file as long_description: +long_description = open('README.rst').read() +# Read in requirements.txt +requirements = open('requirements.txt').readlines() +requirements = [r.strip() for r in requirements] + + +setup( + name='fermilibpluginpsi4', + version=__version__, + author='Ryan Babbush, Jarrod McClean, Damian Steiger, Ian Kivlichan, ' + 'Thomas Haener, Vojtech Havlicek, Matthew Neeley, Wei Sun', + author_email='ryanbabbush@gmail.com, jarrod.mcc@gmail.com, ' + 'fermilib@projectq.ch', + url='http://www.projectq.ch', + description='A plugin allowing FermiLib to interface with Psi4.', + long_description=long_description, + install_requires=requirements, + license='GNU Lesser General Public License (LGPL), Version 3 or any later' + ' version', + packages=find_packages(), + include_package_data=True, + package_data={ + '': [os.path.join('fermilibpluginpsi4', '_psi4_template')] + } +)