Merge pull request #3 from eimrek/release/1.5.0

Release/1.5.0

example.ipynb

@@ -84,7 +84,7 @@
     "index = 0\n",
     "print(\"Orbital index: %d, spin %d\" % (index, spin))\n",
     "print(\"Energy: %.6f eV\" % mfh_model.evals[spin][index])\n",
-    "mfh_model.plot_orbital(mo_index=index, spin=spin)\n",
+    "mfh_model.plot_mo_eigenvector(mo_index=index, spin=spin)\n",
     "\n",
     "# corresponding eigenvector (each element corresponds in order to atoms defined in geom/mfh_model.ase_geom)\n",
     "evec = mfh_model.evecs[spin][index]"
@@ -154,6 +154,59 @@
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Natural orbitals"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geom = ase.io.read(\"geom/clars_goblet.xyz\")\n",
+    "\n",
+    "# \"open shell\" case, normal MFH\n",
+    "mfh_model = tbmfh.MeanFieldHubbardModel(geom, [2.7, 0.1, 0.4], charge=0, multiplicity=1)\n",
+    "mfh_model.run_mfh(u = 3.0, print_iter=False, plot=False)\n",
+    "mfh_model.calculate_natural_orbitals()\n",
+    "\n",
+    "# \"closed shell\" case (just tight-binding)\n",
+    "tb_model = tbmfh.MeanFieldHubbardModel(geom, [2.7, 0.1, 0.4], charge=0, multiplicity=1)\n",
+    "tb_model.run_mfh(u = 0.0, print_iter=False, plot=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_orb = 2\n",
+    "h = 8.0\n",
+    "\n",
+    "for i_rel in np.arange(num_orb, -num_orb, -1):\n",
+    "    i_mo = int(np.around(0.5 * (mfh_model.num_spin_el[0] + mfh_model.num_spin_el[1]))) + i_rel - 1\n",
+    "    \n",
+    "    fig, axs = plt.subplots(nrows=1, ncols=7, figsize=(7 * mfh_model.figure_size[0], mfh_model.figure_size[1]))\n",
+    "    \n",
+    "    mfh_model.plot_no_eigenvector(i_mo, ax=axs[0])\n",
+    "    mfh_model.plot_orb_squared_map(axs[1], mfh_model.no_evecs[i_mo], h=h)\n",
+    "    \n",
+    "    mfh_model.plot_mo_eigenvector(i_mo, spin=0, ax=axs[2])\n",
+    "    mfh_model.plot_mo_eigenvector(i_mo, spin=1, ax=axs[3])\n",
+    "    mfh_model.plot_sts_map(axs[4], mfh_model.evals[0, i_mo], h=h)\n",
+    "    \n",
+    "    tb_model.plot_mo_eigenvector(i_mo, spin=0, ax=axs[5])\n",
+    "    tb_model.plot_sts_map(axs[6], tb_model.evals[0, i_mo], h=h)\n",
+    "    \n",
+    "    plt.subplots_adjust(wspace=0.0, hspace=0)\n",
+    "    plt.show()"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -179,6 +232,19 @@
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.7.9"
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
   }
  },
  "nbformat": 4,

geom/phenalenyl-dimer.xyz

@@ -0,0 +1,28 @@
+26
+Properties=species:S:1:pos:R:3:tags:I:1 pbc="F F F"
+C       22.69117931       8.32750000       3.99999999        2
+C       24.02485843       9.09750002       3.99999998        1
+C       20.02382107       8.32749996       3.99999999        2
+C       18.69014194       9.09749995       3.99999999        1
+C       18.69014191      10.63749995       3.99999999        2
+C       20.02382104      11.40749996       3.99999998        1
+C       21.35750017      10.63749999       3.99999998        2
+C       21.35750019       9.09750000       3.99999999        1
+C       24.02485841      10.63750003       3.99999998        2
+C       22.69117927      11.40750001       3.99999998        1
+C       22.69117926      12.94750001       3.99999997        2
+C       21.35750011      13.71749999       3.99999997        1
+C       20.02382101      12.94749998       3.99999998        2
+C       20.02382101      20.66249999       3.99999999        1
+C       18.69014189      19.89249997       3.99999998        2
+C       22.69117925      20.66250003       3.99999999        1
+C       24.02485838      19.89250004       3.99999999        2
+C       24.02485840      18.35250004       3.99999999        1
+C       22.69117928      17.58250003       3.99999998        2
+C       21.35750015      18.35250000       3.99999998        1
+C       21.35750013      19.89249999       3.99999999        2
+C       18.69014190      18.35249996       3.99999998        1
+C       20.02382104      17.58249998       3.99999998        2
+C       20.02382106      16.04249998       3.99999997        1
+C       21.35750020      15.27250000       3.99999997        2
+C       22.69117931      16.04250001       3.99999998        1

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="tb-mean-field-hubbard",
-    version="1.4.1",
+    version="1.5.0",
     author="Kristjan Eimre",
     author_email="[email protected]",
     description="Package to run tight-binding mean field hubbard calculations",

tb_mean_field_hubbard/__init__.py

@@ -1,3 +1,3 @@
 from .mfh import MeanFieldHubbardModel
 
-__version__ = "1.4.1"
+__version__ = "1.5.0"

tb_mean_field_hubbard/mfh.py

@@ -3,15 +3,11 @@
 import matplotlib.pyplot as plt
 import matplotlib.patches as mpatches
 
-import os
-
 import ase
 import ase.neighborlist
 
 import pythtb
 
-import sys
-
 from . import utils
 
 ### ------------------------------------------------------------------------------
@@ -36,7 +32,7 @@ def __init__(self, ase_geom, t_list=[2.7], charge=0, multiplicity='auto', bond_c
         self.num_atoms = len(ase_geom)
 
         self.figure_size = (np.ptp(self.ase_geom.positions, axis=0)[:2] + 1.0) / 4.0
-        self.figure_size[0] += 2.5
+        self.figure_size[0] = max([self.figure_size[0], 3.5])  # to have enough space for titles
 
         self.spin_guess = self._load_spin_guess(self.ase_geom)
 
@@ -58,6 +54,10 @@ def __init__(self, ase_geom, t_list=[2.7], charge=0, multiplicity='auto', bond_c
         self.evals = None
         self.evecs = None
 
+        # natural orbital evals and evecs
+        self.no_evals = None
+        self.no_evecs = None
+
     ### ------------------------------------------------------------------------------
     ### TB routines
     ### ------------------------------------------------------------------------------
@@ -131,7 +131,7 @@ def visualize_tb_model(self):
     def visualize_spin_guess(self):
         plt.figure(figsize=self.figure_size)
         spin_guess_plot = 0.5 * (self.spin_guess[:, 0] - self.spin_guess[:, 1])
-        utils.make_plot(plt.gca(), self.ase_geom, self.neighbor_list, spin_guess_plot)
+        utils.make_evec_plot(plt.gca(), self.ase_geom, self.neighbor_list, spin_guess_plot)
         plt.show()
 
     def print_parameters(self):
@@ -319,7 +319,7 @@ def _get_atoms_extent(self, atoms, edge_space):
         ymax = np.max(atoms.positions[:, 1]) + edge_space
         return xmin, xmax, ymin, ymax
 
-    def calc_orb_map(self, i_spin, i_orb, h=10.0, edge_space=5.0, dx=0.1, z_eff=3.25):
+    def calc_orb_map(self, evec, h=10.0, edge_space=5.0, dx=0.1, z_eff=3.25):
 
         extent = self._get_atoms_extent(self.ase_geom, edge_space)
 
@@ -329,15 +329,21 @@ def calc_orb_map(self, i_spin, i_orb, h=10.0, edge_space=5.0, dx=0.1, z_eff=3.25
 
         orb_map = np.zeros((len(x_arr), len(y_arr)), dtype=np.complex)
 
-        evc = self.evecs[i_spin][i_orb]
-        for at, coef in zip(self.ase_geom, evc):
+        for at, coef in zip(self.ase_geom, evec):
             p = at.position
             local_i, local_grid = utils.get_local_grid(x_arr, y_arr, p, cutoff=1.2 * h + 4.0)
             pz_orb = utils.carbon_2pz_slater(local_grid[0] - p[0], local_grid[1] - p[1], h, z_eff)
             orb_map[local_i[0]:local_i[1], local_i[2]:local_i[3]] += coef * pz_orb
 
         return orb_map
 
+    def plot_orb_squared_map(self, ax, evec, h=10.0, edge_space=5.0, dx=0.1, title=None, cmap='seismic', z_eff=3.25):
+        orb_map = np.abs(self.calc_orb_map(evec, h, edge_space, dx, z_eff))**2
+        extent = self._get_atoms_extent(self.ase_geom, edge_space)
+        ax.imshow(orb_map.T, origin='lower', cmap=cmap, extent=extent)
+        ax.axis('off')
+        ax.set_title(title)
+
     def calc_sts_map(self, energy, broadening=0.05, h=10.0, edge_space=5.0, dx=0.1, z_eff=3.25):
 
         extent = self._get_atoms_extent(self.ase_geom, edge_space)
@@ -351,7 +357,8 @@ def calc_sts_map(self, energy, broadening=0.05, h=10.0, edge_space=5.0, dx=0.1,
             for i_orb, evl in enumerate(self.evals[i_spin]):
                 if np.abs(energy - evl) <= 3.0 * broadening:
                     broad_coef = utils.gaussian(energy - evl, broadening)
-                    orb_ldos_map = np.abs(self.calc_orb_map(i_spin, i_orb, h, edge_space, dx, z_eff))**2
+                    evec = self.evecs[i_spin][i_orb]
+                    orb_ldos_map = np.abs(self.calc_orb_map(evec, h, edge_space, dx, z_eff))**2
                     final_map += broad_coef * orb_ldos_map
         return final_map
 
@@ -374,6 +381,28 @@ def plot_sts_map(self,
         ax.axis('off')
         ax.set_title(title)
 
+    def plot_eigenvector(self, ax, evec, title=None):
+        utils.make_evec_plot(ax, self.ase_geom, self.neighbor_list, evec, title=title)
+
+    def plot_mo_eigenvector(self, mo_index, spin=0, ax=None):
+        title = "mo%d s%d %s, en: %.2f" % (mo_index, spin, utils.orb_label(
+            mo_index, self.num_spin_el[spin]), self.evals[spin][mo_index])
+        if ax is None:
+            plt.figure(figsize=self.figure_size)
+            self.plot_eigenvector(plt.gca(), self.evecs[spin][mo_index], title=title)
+            plt.show()
+        else:
+            self.plot_eigenvector(ax, self.evecs[spin][mo_index], title=title)
+
+    def plot_no_eigenvector(self, no_index, ax=None):
+        title = "no%d, occ=%.4f" % (no_index, self.no_evals[no_index])
+        if ax is None:
+            plt.figure(figsize=self.figure_size)
+            self.plot_eigenvector(plt.gca(), self.no_evecs[no_index], title=title)
+            plt.show()
+        else:
+            self.plot_eigenvector(ax, self.no_evecs[no_index], title=title)
+
     def report(self, num_orb=2, sts_h=10.0, sts_broad=0.05):
 
         print(f"multiplicity:       {self.multiplicity:12d}")
@@ -382,7 +411,7 @@ def report(self, num_orb=2, sts_h=10.0, sts_broad=0.05):
         print("---")
         print("spin density:")
         plt.figure(figsize=self.figure_size)
-        utils.make_plot(plt.gca(), self.ase_geom, self.neighbor_list, self.spin_density)
+        utils.make_evec_plot(plt.gca(), self.ase_geom, self.neighbor_list, self.spin_density)
         plt.show()
 
         print("---")
@@ -406,23 +435,35 @@ def report(self, num_orb=2, sts_h=10.0, sts_broad=0.05):
 
             _fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(4 * self.figure_size[0], self.figure_size[1]))
 
-            title1 = "mo%d α %s, en: %.2f" % (i_mo, utils.orb_label(i_mo, self.num_spin_el[0]), self.evals[0][i_mo])
-            utils.make_plot(axs[0], self.ase_geom, self.neighbor_list, self.evecs[0][i_mo], title1)
-
-            title2 = "mo%d β %s, en: %.2f" % (i_mo, utils.orb_label(i_mo, self.num_spin_el[1]), self.evals[1][i_mo])
-            utils.make_plot(axs[1], self.ase_geom, self.neighbor_list, self.evecs[1][i_mo], title2)
+            self.plot_mo_eigenvector(i_mo, spin=0, ax=axs[0])
+            self.plot_mo_eigenvector(i_mo, spin=1, ax=axs[1])
 
-            title3 = "sts h=%.1f, en: %.2f" % (sts_h, self.evals[0][i_mo])
-            self.plot_sts_map(axs[2], self.evals[0][i_mo], broadening=sts_broad, h=sts_h, title=title3)
+            title1 = "sts h=%.1f, en: %.2f" % (sts_h, self.evals[0][i_mo])
+            self.plot_sts_map(axs[2], self.evals[0][i_mo], broadening=sts_broad, h=sts_h, title=title1)
 
-            title4 = "sts h=%.1f, en: %.2f" % (sts_h, self.evals[1][i_mo])
-            self.plot_sts_map(axs[3], self.evals[1][i_mo], broadening=sts_broad, h=sts_h, title=title4)
+            title2 = "sts h=%.1f, en: %.2f" % (sts_h, self.evals[1][i_mo])
+            self.plot_sts_map(axs[3], self.evals[1][i_mo], broadening=sts_broad, h=sts_h, title=title2)
 
             plt.show()
 
-    def plot_orbital(self, mo_index, spin=0):
-        title = "mo%d s%d %s, en: %.2f" % (mo_index, spin, utils.orb_label(
-            mo_index, self.num_spin_el[spin]), self.evals[spin][mo_index])
-        plt.figure(figsize=self.figure_size)
-        utils.make_plot(plt.gca(), self.ase_geom, self.neighbor_list, self.evecs[spin][mo_index], title=title)
-        plt.show()
+    def calculate_natural_orbitals(self):
+        # build the one particle reduced density matrix
+        dens_mat = None
+
+        for i_spin in range(2):
+            for i_el in range(self.num_spin_el[i_spin]):
+                evec = self.evecs[i_spin, i_el]
+                if dens_mat is None:
+                    dens_mat = np.outer(evec, np.conj(evec))
+                else:
+                    dens_mat += np.outer(evec, np.conj(evec))
+
+        # Diagonalize the density matrix
+        self.no_evals, self.no_evecs = np.linalg.eig(dens_mat)
+        self.no_evals = np.abs(self.no_evals)
+        self.no_evecs = self.no_evecs.T
+
+        # sort the natural orbitals based on occupations
+        sort_inds = (-1 * self.no_evals).argsort()
+        self.no_evals = self.no_evals[sort_inds]
+        self.no_evecs = self.no_evecs[sort_inds]

tb_mean_field_hubbard/utils.py

@@ -150,7 +150,7 @@ def visualize_evec(ax, atoms, evec):
     #print(" Vol: %.6f"%vol)
 
 
-def make_plot(ax, atoms, neighbor_list, data, title=None, filename=None):
+def make_evec_plot(ax, atoms, neighbor_list, data, title=None, filename=None):
 
     ax.set_aspect('equal')
     visualize_backbone(ax, atoms, neighbor_list)