first commit from spinv

.github/workflows/build.yml

@@ -0,0 +1,22 @@
+# This is a basic workflow to help you get started with Actions
+name: Python package
+
+on: [push]
+
+# A workflow run is made up of one or more jobs that can run sequentially or in parallel
+jobs:
+  # This workflow contains a single job called "build"
+  build:
+    # The type of runner that the job will run on
+    runs-on: ubuntu-latest
+
+    # Steps represent a sequence of tasks that will be executed as part of the job
+    steps:
+      # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
+      - uses: actions/checkout@v2
+
+      - name: Check pytest tests
+        run: |
+          pip install . 
+          pip3 install pytest 
+          pytest

.gitignore

@@ -0,0 +1,6 @@
+__pycache__/
+strawberrypy/__pycache__/
+strawberrypy/_pythtb/__pycache__/
+strawberrypy/_tbmodels/__pycache__/
+.devcontainer/
+.pytest_cache

LICENSE.txt

@@ -0,0 +1,23 @@
+MIT License
+
+Copyright (c) 2023 Universita' di Trieste, Italy.
+All rights reserved.
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to 
+deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following condition:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,166 @@
+# StraWBerryPy
+StraWBerryPy (Single-poinT and local invaRiAnts for Wannier Berriologies in Python) is a Python package calculating topological invariants for non-crystalline 2D topological insulators. The single-point formulas in the supercell framework for the Chern number [[Ceresoli-Resta](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.76.012405)] and for the spin Chern number [[Favata-Marrazzo](https://iopscience.iop.org/article/10.1088/2516-1075/acba6f/meta)] are implemented in StraWBerryPy. 
+In addition, StraWBerryPy can handle finite systems (such as bounded samples and heterostructures) and compute the local Chern marker [[Bianco-Resta](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.84.241106)].
+The code provides dedicated interfaces to tight-binding packages [PythTB](http://www.physics.rutgers.edu/pythtb/) and [Tbmodels](https://tbmodels.greschd.ch/en/latest/).
+
+
+## Quick start
+Here, two examples for calculating the single-point topological invariant in the supercell framework for tight-binding models in presence of Anderson disorder. 
+
+### Kane-Mele example
+Let's consider as prototype of quantum spin Hall insulator the Kane-Mele model and calculate the topological invariant through the single-point spin Chern number formulation.  
+
+First create the tight-binding model in primitive cell using PythTB and Tbmodels packages. Create a supercell, for example including L lattice points in each lattice vector direction.  
+Then add disorder in the supercell, for instance Anderson disorder which is a uniformly distributed random on site potential.  
+For a reference implementation, look at [Kane-Mele example](strawberrypy/example_models/kane_mele.py).
+
+```python
+import numpy as np
+from strawberrypy import single_point_spin_chern
+from strawberrypy.example_models import kane_mele_pythtb, kane_mele_tbmodels, km_anderson_disorder_pythtb, km_anderson_disorder_tbmodels
+
+L = 12              # supercell LxL linear size
+
+# Topological phase
+r = 1.              # r = rashba/spin_orb
+e = 3.              # e = e_onsite/spin_orb
+spin_o = 0.3        # spin_orb = spin_orb/t (t=1)
+
+# # Trivial phase
+# r = 3. 
+# e = 5.5
+# spin_o = 0.3
+
+# Create Kane-Mele model in supercell LxL through PythTB package
+km_pythtb_model = kane_mele_pythtb(r, e, spin_o, L)
+
+# Create Kane-Mele model in supercell LxL through TBmodels package
+km_tbmodels_model = kane_mele_tbmodels(r, e, spin_o, L)
+
+w = 1.              # w = disorder strength such that random potential is in [-w/2, w/2]      
+
+# Add Anderson disorder in PythTB model
+np.random.seed(10)
+km_pythtb_model = km_anderson_disorder_pythtb(km_pythtb_model, w)
+
+# Add Anderson disorder in TBmodels model
+np.random.seed(10)
+km_tbmodels_model = km_anderson_disorder_tbmodels(km_tbmodels_model, w)
+```
+Finally, use single_point_spin_chern function to calculate the single-point invariant for the disordered model in the supercell framework.   
+The function takes as inputs:
+- the model created with PythTB or TBmodels packages 
+- variable "spin" indicating which spin Chern number we want to calculate : 'up' or 'down'
+- variable "formula" selecting which single-point formula we want to calculate : 'asymmetric', 'symmetric' or 'both'
+
+```python
+spin_chern = 'down'
+which_formula = 'symmetric'
+
+# Single Point Spin Chern Number (SPSCN) calculation for Pythtb model
+spin_chern_pythtb = single_point_spin_chern(km_pythtb_model, spin=spin_chern, formula=which_formula)
+
+# Single Point Spin Chern Number (SPSCN) calculation for TBmodels model
+spin_chern_tbmodels = single_point_spin_chern(km_tbmodels_model, spin=spin_chern, formula=which_formula)
+
+# If which_formula = 'both', then Single Point Spin Chern numbers are printed as follows : 'asymmetric' 'symmetric'
+print('PythTB package, supercell size L =', L, ' disorder strength = ', w,  ' SPSCN :', *spin_chern_pythtb )
+print('TBmodels package, supercell size L =', L, ' disorder strength = ', w,  ' SPSCN :', *spin_chern_tbmodels )
+```
+
+### Haldane example
+The single-point invariant calculation can be also performed for a Chern insulator, like the Haldane model, using single_point_chern function.   
+For the model implementation, look at [Haldane example](strawberrypy/example_models/haldane.py).
+```python
+import numpy as np
+from strawberrypy import single_point_chern
+from strawberrypy.example_models import haldane_pythtb, haldane_tbmodels, h_anderson_disorder_pythtb, h_anderson_disorder_tbmodels
+
+L = 12              # supercell LxL linear size
+
+# Topological phase
+t = -4.                   # t = first neighbours real hopping
+t2 = 1.                   # t2 = second neighbours hopping
+delta = 2.                # delta = energy on site
+phi = -np.pi/2.0          # phi = second neighbours hopping phase
+
+# # Trivial phase
+# t = -4.                   # t = first neighbours real hopping
+# t2 = 1.                   # t2 = second neighbours hopping
+# delta = 4.5               # delta = energy on site
+# phi = -np.pi/6.0          # phi = second neighbours hopping phase
+
+# Create Haldane model in supercell LxL through PythTB package
+h_pythtb_model = haldane_pythtb(delta, t, t2, phi, L)
+
+# Create Haldane model in supercell LxL through TBmodels package
+h_tbmodels_model = haldane_tbmodels(delta, t, t2, phi, L)
+
+w = 1.                      # w = disorder strength such that random potential is in [-w/2, w/2]
+
+# Add Anderson disorder in PythTB model
+np.random.seed(10)
+h_pythtb_model = h_anderson_disorder_pythtb(h_pythtb_model, w)
+
+# Add Anderson disorder in TBmodels model
+np.random.seed(10)
+h_tbmodels_model = h_anderson_disorder_tbmodels(h_tbmodels_model, w)
+
+which_formula = 'both'
+
+# Single Point Chern Number (SPCN) calculation for Pythtb model
+chern_pythtb = single_point_chern(h_pythtb_model, formula=which_formula)
+
+# Single Point Chern Number (SPCN) calculation for TBmodels model
+chern_tbmodels = single_point_chern(h_tbmodels_model, formula=which_formula)
+
+# If which_formula = 'both', then Single Point Chern Numbers are printed as follows : 'asymmetric' 'symmetric'
+print('PythTB package, supercell size L =', L, ' disorder strength = ', w,  ' SPCN :', *chern_pythtb )
+print('TBmodels package, supercell size L =', L, ' disorder strength = ', w,  ' SPCN :', *chern_tbmodels )
+```
+## Local invariants
+Here an example for calculating the local Chern marker for a bounded Haldane model in presence of Anderson disorder. First create the tight-binding model in the primitive cell using PythTB and TBmodels packages. Then, the models has to be cut in both x and y directions specifying the size of the resulting samples.
+To conclude the construction of the models the disorder can be added, for instance Anderson disorder which is a uniformly distributed random on site potential.
+```python
+import numpy as np
+from strawberrypy import local_chern_marker, make_finite, onsite_disorder
+from strawberrypy.example_models import haldane_pythtb, haldane_tbmodels
+
+# Create Haldane models through PythTB and TBmodels packages
+hmodel_pbc_tbmodels = haldane_tbmodels(delta = 0.5, t = -1, t2 = 0.15, phi = np.pi / 2, L = 1)
+hmodel_pbc_pythtb = haldane_pythtb(delta = 0.5, t = -1, t2 = 0.15, phi = np.pi / 2, L = 1)
+
+# Cut the models to make a sample of size 10 x 10
+hmodel_obc_tbmodels = make_finite(model = hmodel_pbc_tbmodels, nx_sites = 10, ny_sites = 10)
+hmodel_obc_pythtb = make_finite(model = hmodel_pbc_pythtb, nx_sites = 10, ny_sites = 10)
+
+# Add Anderson disorder within [-w/2, w/2] to the samples. The argument spinstates specifies the spin of the model
+hmodel_tbm_disorder = onsite_disorder(model = hmodel_tbm, w = 1, spinstates = 1, seed = 184)
+hmodel_pythtb_disorder = onsite_disorder(model = hmodel_pythtb, w = 1, spinstates = 1, seed = 184)
+```
+Finally the local Chern marker can be computed using ```local_chern_marker```, which takes as arguments:
+- the model created using PythTB or TBmodels (using ```model```)
+- the size of the bounded sample (via ```nx_sites``` and ```ny_sites```)
+- the direction along which to compute the local marker (if not specified otherwise the function computes the local marker for the whole lattice) using the argument ```direction```, and the cell along the orthogonal direction from which to start (via ```start```)
+- a boolean value, ```return_projector```, which specifies if return also the ground state projector of the model, and the argument ```projector``` to input it in order to avoid recalculation of heavy objects
+- the possibility to perform macroscopic averages for disordered samples, via the boolean ```macroscopic_average```, and the cutoff radius of the average using the argument ```cutoff```
+
+In the example below, the local Chern marker is evaluated along the x direction, starting from the cell whose y reduced coordinate is half the dimension of the lattice, and is averaged over a region with cutoff 2:
+```python
+# Compute the local Chern markers for TBmodels and PythTB
+lcm_tbmodels = local_chern_marker(model = hmodel_tbm_disorder, nx_sites = 10, ny_sites = 10, direction = 0, start = 5, macroscopic_average = True, cutoff = 2)
+lcm_pythtb = local_chern_marker(model = hmodel_pythtb_disorder, nx_sites = 10, ny_sites = 10, direction = 0, start = 5, macroscopic_average = True, cutoff = 2)
+
+print("Local Chern marker along the x direction starting from y=5, computed using TBmodels: ", lcm_tbmodels)
+print("Local Chern marker along the x direction starting from y=5, computed using PythTB: ", lcm_pythtb)
+```
+
+## Installation
+```
+git clone https://github.com/strawberrypy-developers/strawberrypy.git
+cd strawberrypy
+pip install .
+```
+
+## Authors
+Roberta Favata and Nicolas Baù and Antimo Marrazzo

setup.py

@@ -0,0 +1,14 @@
+import setuptools
+
+setuptools.setup(
+    name='strawberrypy',
+    version='0.2.0', 
+    author='Roberta Favata and Nicolas Baù and Antimo Marrazzo',   
+    author_email='[email protected]',
+    description='Python package for calculation of topological invariants through single-point formulas and local markers',
+    license='MIT License',
+    packages=setuptools.find_packages(),
+    install_requires=['pythtb==1.8.0',
+                      'tbmodels==1.4.3',
+                      'opt-einsum==3.3.0',],
+)

strawberrypy/__init__.py

@@ -0,0 +1,5 @@
+from .package import single_point_spin_chern, single_point_chern
+from .package import make_finite, make_heterostructure
+from .package import onsite_disorder
+from .package import local_chern_marker
+from . import example_models

strawberrypy/_pythtb/__init__.py

@@ -0,0 +1,4 @@
+from .singlepoint_invariant import single_point_chern, single_point_spin_chern
+from .finite_systems import make_finite, make_heterostructure
+from .lattice import onsite_disorder
+from .localmarkers import local_chern_marker

strawberrypy/_pythtb/finite_systems.py

@@ -0,0 +1,91 @@
+import numpy as np
+import pythtb as ptb
+
+def make_finite(model, nx_sites: int, ny_sites: int):
+    """
+    Make a finite mdoel along x and y direction by first cutting on the y direction and then on the x. This convention has been used to track the positions in the functions
+
+        Args:
+        - model : instance of the model, which should be periodic in both x and y direction
+        - nx_sites, ny_sites : number of sites of the finite sample in both directions
+
+        Returns:
+        - model : the finite model
+    """
+
+    if not (nx_sites > 0 and ny_sites > 0):
+        raise RuntimeError("Number of sites along finite direction must be greater than 0")
+
+    ribbon = model.cut_piece(num = ny_sites, fin_dir = 1, glue_edgs = False)
+    finite = ribbon.cut_piece(num = nx_sites, fin_dir = 0, glue_edgs = False)
+
+    return finite
+
+def make_heterostructure(model1 , model2,  nx_sites : int, ny_sites : int, direction : int, start : int, stop : int):
+    """
+    Modify a finite model by merging another system in it. The system will be split in the direction starting from start.
+
+        Args:
+        - model1, model2: the models which composes the heterostructure
+        - nx_sites, ny_sites : x and y length of the finite system
+        - direction : direction in which the splitting happen, allowed 0 for 'x' or 1 for 'y'
+        - start : starting point for the splitting in the 'direction' direction
+        - end : end point of the splitting in the 'direction' direction
+
+        Returns:
+        - model : the model composed my the two subsystems
+    """
+
+    # Number of states per unit cell
+    atoms_uc = int(model1.get_num_orbitals() / (nx_sites * ny_sites))
+
+    # Check validity of input data
+    if not start < stop:
+        raise RuntimeError("Starting point is greater or equal to the end point")
+    if not (start >= 0 and start < (nx_sites if direction == 0 else ny_sites)):
+        raise RuntimeError("Start point value not allowed")
+    if not (stop > 0 and stop < (nx_sites if direction == 0 else ny_sites)):
+        raise RuntimeError("End point value not allowed")
+    if direction not in [0, 1]:
+        raise RuntimeError("Direction not allowed: insert 0 for 'x' and 1 for 'y'")
+
+    # Assert the model is the same
+    if not (model1.get_num_orbitals() == model2.get_num_orbitals() and np.allclose(model1._orb, model2._orb) and np.allclose(model1._lat, model2._lat)):
+        raise RuntimeError("The models to merge must be the same model with different parameters")
+
+    # Make a copy of the model and remove all onsite terms
+    onsite1 = np.copy(model1._site_energies)
+    onsite2 = np.copy(model2._site_energies)
+
+    if direction == 0:
+        # Splitting along the x direction
+        ind = np.array([[(start + i) * ny_sites * atoms_uc + j * atoms_uc for j in range(ny_sites)] for i in range(stop - start + 1)]).flatten()
+    else:
+        # Splitting along the y direction
+        ind = np.array([[i * ny_sites * atoms_uc + start * atoms_uc + j * atoms_uc for j in range(stop - start + 1)] for i in range(nx_sites)]).flatten()
+
+    for i in ind:
+        for j in range(atoms_uc):
+            onsite1[i + j] = onsite2[i + j]
+
+    # Indices of every atom in the selected cells, not only of the initial atom of the cell
+    indices = np.array([[i + j for j in range(atoms_uc)] for i in ind]).flatten()
+
+    # Hopping amplitudes and positions
+    hoppings1 = model1._hoppings; hoppings2 = model2._hoppings
+
+    # The model to return
+    newmodel = ptb.tb_model(dim_k = 0, dim_r = model1._dim_r, lat = model1._lat, orb = model1._orb, nspin = model1._nspin)
+    newmodel.set_onsite(onsite1, mode = 'set')
+
+    # Cycle over the rows of the hopping matrix
+    for k in range(len(hoppings1)):
+
+        if k in indices:
+            if np.absolute(hoppings2[k][0]) < 1e-10: continue
+            newmodel.set_hop(hoppings2[k][0], hoppings2[k][1], hoppings2[k][2], mode = "add")
+        else:
+            if np.absolute(hoppings1[k][0]) < 1e-10: continue
+            newmodel.set_hop(hoppings1[k][0], hoppings1[k][1], hoppings1[k][2], mode = "add")
+
+    return newmodel