diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
new file mode 100644
index 0000000..ccebecd
--- /dev/null
+++ b/.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
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..54a5444
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+__pycache__/
+strawberrypy/__pycache__/
+strawberrypy/_pythtb/__pycache__/
+strawberrypy/_tbmodels/__pycache__/
+.devcontainer/
+.pytest_cache
\ No newline at end of file
diff --git a/LICENSE.txt b/LICENSE.txt
new file mode 100644
index 0000000..5959d97
--- /dev/null
+++ b/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.
\ No newline at end of file
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..c9c5cdf
--- /dev/null
+++ b/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
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..4d98c09
--- /dev/null
+++ b/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='favata.roberta@gmail.com',
+    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',],
+)
\ No newline at end of file
diff --git a/strawberrypy/__init__.py b/strawberrypy/__init__.py
new file mode 100644
index 0000000..32efe72
--- /dev/null
+++ b/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
diff --git a/strawberrypy/_pythtb/__init__.py b/strawberrypy/_pythtb/__init__.py
new file mode 100644
index 0000000..7359c02
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/strawberrypy/_pythtb/finite_systems.py b/strawberrypy/_pythtb/finite_systems.py
new file mode 100644
index 0000000..cd435cc
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/strawberrypy/_pythtb/lattice.py b/strawberrypy/_pythtb/lattice.py
new file mode 100644
index 0000000..d93335b
--- /dev/null
+++ b/strawberrypy/_pythtb/lattice.py
@@ -0,0 +1,102 @@
+import numpy as np
+import pythtb
+
+def reciprocal_vec(model):
+    #returns reciprocal lattice vectors
+    lat = model.get_lat()     
+    a_matrix = np.array([[lat[1,1], -lat[0,1]],[-lat[1,0], lat[0,0]]])
+    b_matrix = (2.*np.pi / (lat[0,0]*lat[1,1]-lat[0,1]*lat[1,0])) * a_matrix
+
+    b1 = b_matrix[:,0].reshape(-1,1)
+    b2 = b_matrix[:,1].reshape(-1,1)
+    return b1, b2
+
+def orb_cart (model):
+    #returns position of orbitals in cartesian coordinates
+    n_orb = model.get_num_orbitals()
+    lat_super = model.get_lat()          
+    orb_red = model.get_orb()            
+
+    orb_c = []
+    for i in range (n_orb):
+        orb_c.append( (np.matmul(lat_super.transpose(),orb_red[i].reshape(-1,1))).squeeze() )   
+    orb_c = np.array(orb_c)
+    return orb_c
+
+def orb_cart_spin (model):
+    #returns position of orbitals in cartesian coordinates
+    n_occ = model.get_num_orbitals()
+    lat_super = model.get_lat()          
+    orb_red = model.get_orb()            
+
+    orb_c = []
+    for i in range (n_occ):
+        orb_c.append( (np.matmul(lat_super.transpose(),orb_red[i].reshape(-1,1))).squeeze() )  
+        orb_c.append( (np.matmul(lat_super.transpose(),orb_red[i].reshape(-1,1))).squeeze() )  
+    orb_c = np.array(orb_c)
+    return orb_c
+
+def periodic_gauge (u_n0, b, model):
+    #returns the matrix of occupied eigenvectors at the edge of the Brillouin zone imposing periodic gauge
+    n_orb = model.get_num_orbitals()
+    n_occ = n_orb//2
+
+    orb_c = orb_cart(model)
+    vec_scal_b = orb_c @ b
+    vec_exp_b = np.exp(-1.j*vec_scal_b)
+
+    u_nb = vec_exp_b.T * u_n0[:n_occ,:]
+
+    return u_nb
+
+def periodic_gauge_spin(u_n0, b, model):
+    #returns the matrix of occupied eigenvectors at the edge of the Brillouin zone imposing periodic gauge
+    n_orb = model.get_num_orbitals()
+    n_occ = n_orb
+
+    orb_c = orb_cart_spin(model)
+    vec_scal_b = orb_c @ b
+    vec_exp_b = np.exp(-1.j*vec_scal_b)
+
+    u_nb = vec_exp_b.T * u_n0[:n_occ,:]
+
+    return u_nb
+
+def onsite_disorder(source_model, w : float, spinstates : int = 2, seed : int = None):
+    """
+    Add onsite (Anderson) disorder to the specified model. The disorder amplitude per site is taken randomly in [-w/2, w/2].
+
+        Args:
+        - source_model : the model to add disorder to
+        - w : disorder amplitude
+        - spinstates : spin of the model
+        - seed : seed for random number generator
+
+        Returns:
+        - model : the disordered model
+    """
+    # Quick return for no disorder
+    if w == 0:
+        return source_model
+
+    if seed is not None:
+        np.random.seed(seed)
+
+    # Number of orbitals in the supercell model = norbs (original) x num
+    norbs = source_model.get_num_orbitals()
+    
+    # Onsite energies per unit cell (2 is by convention with TBModels)
+    disorder = 0.5 * w * (2 * np.random.rand(norbs // spinstates) - 1.0)
+    disorder = np.repeat(disorder, spinstates)
+    onsite = source_model._site_energies + disorder
+
+    newmodel = pythtb.tb_model(dim_k = 0, dim_r = source_model._dim_r, lat = source_model._lat, orb = source_model._orb, nspin = source_model._nspin)
+    newmodel.set_onsite(onsite, mode = 'set')
+
+    # Cycle over the rows of the hopping matrix
+    hoppings = source_model._hoppings
+    for k in range(len(hoppings)):
+        if np.absolute(hoppings[k][0]) < 1e-10: continue
+        newmodel.set_hop(hoppings[k][0], hoppings[k][1], hoppings[k][2], mode = "add")
+    
+    return newmodel
\ No newline at end of file
diff --git a/strawberrypy/_pythtb/localmarkers.py b/strawberrypy/_pythtb/localmarkers.py
new file mode 100644
index 0000000..553fa47
--- /dev/null
+++ b/strawberrypy/_pythtb/localmarkers.py
@@ -0,0 +1,93 @@
+import pythtb as ptb
+import numpy as np
+import scipy.linalg as la
+from opt_einsum import contract
+
+from .finite_systems import *
+from strawberrypy.common_func import *
+from .lattice import orb_cart
+
+def local_chern_marker(model, nx_sites : int, ny_sites : int, direction : int = None, start : int = 0, return_projector : bool = False, projector : np.ndarray = None, macroscopic_average : bool = False, cutoff : float = 0.8):
+    """
+    Evaluate the Z Chern marker on the whole lattice if direction is None. If direction is not None evaluates the Z Chern marker along direction starting from start.
+    
+        Args:
+        - model : instance of TBModels of the model
+        - nx_sites : number of unit cells in the x direction of the finite system
+        - ny_sites : number of unit cells in the y direction of the finite system
+        - direction : direction along which compute the local Chern marker, default is None (returns the whole lattice Chern marker), allowed values are 0 for 'x' direction and 1 for 'y' direction
+        - start : if direction is not None, is the coordinate of the unit cell at the start of the evaluation of the Chern marker
+        - return_projector : if True, returns the ground state projector at the end of the calculation, default is False
+        - projector : input the ground state projector to be used in the calculation. Default is None, which means it is computed from the model
+        - macroscopic_average : if True, returns the local Chern marker averaged over a radius equal to the cutoff
+        - cutoff : cutoff set for the calculation of averages
+
+        Returns:
+        - lattice_chern : local Chern marker of the whole lattice if direction is None
+        - lcm_direction : local Chern marker along direction starting from start
+        - projector : ground state projector, returned if return_projector is set True (default is False)
+    """
+
+    # Check input variables
+    if direction not in [None, 0, 1]:
+        raise RuntimeError("Direction allowed are None, 0 (which stands for x), and 1 (which stands for y)")
+    
+    if direction is not None:
+        if direction == 0:
+            if start not in range(ny_sites): raise RuntimeError("Invalid start parameter (must be within [0, ny_sites - 1])")
+        else:
+            if start not in range(nx_sites): raise RuntimeError("Invalid start parameter (must be within [0, nx_sites - 1])")
+
+    # Atoms in the unit cell
+    atoms_uc = int(model.get_num_orbitals() / (nx_sites * ny_sites))
+
+    # 2D unit cell volume
+    uc_vol = np.linalg.norm(np.cross(model._lat[0], model._lat[1]))
+
+    if projector is None:
+        # Eigenvectors at \Gamma
+        _, eigenvecs = model.solve_one(eig_vectors = True)
+        eigenvecs = eigenvecs.T
+
+        # Build the ground state projector
+        gs_projector = contract("ji,ki->jk", eigenvecs[:, :int(0.5 * len(eigenvecs))], eigenvecs[:, :int(0.5 * len(eigenvecs))].conjugate())
+    else:
+        gs_projector = projector
+
+    # Position of the sites in the lattice
+    positions = orb_cart(model)
+
+    # Position operator in tight-binding approximation (site orbitals basis)
+    rx = np.diag(positions[:, 0]); ry = np.diag(positions[:, 1])
+
+    # Chern marker operator
+    chern_operator = np.imag(gs_projector @ commutator(rx, gs_projector) @ commutator(ry, gs_projector))
+    chern_operator *= -4 * np.pi / uc_vol
+
+    # If macroscopic average I have to compute the lattice values with the averages first
+    if macroscopic_average:
+        chern_on_lattice = average_over_radius(np.diag(chern_operator), positions.T[0], positions.T[1], cutoff, nx_sites, ny_sites, atoms_uc)
+    
+    if direction is not None:
+        # Evaluate index of the selected direction
+        indices = xy_to_index('x' if direction == 1 else 'y', start, nx_sites, ny_sites, atoms_uc)
+
+        # If macroscopic average consider the averaged lattice, else the Chern operator
+        if macroscopic_average:
+            lcm_direction = [chern_on_lattice[int(indices[atoms_uc * i] / atoms_uc)] for i in range(int(len(indices) / atoms_uc))]
+        else:
+            lcm_direction = [np.sum([chern_operator[indices[atoms_uc * i + j], indices[atoms_uc * i + j]] for j in range(atoms_uc)]) for i in range(int(len(indices) / atoms_uc))]
+        
+        if not return_projector:
+            return np.array(lcm_direction)
+        else:
+            return np.array(lcm_direction), gs_projector
+
+    if not macroscopic_average:
+        chern_on_lattice = partialtrace(chern_operator, atoms_uc)
+
+    lattice_chern = chern_on_lattice.reshape(nx_sites, ny_sites).T
+    if not return_projector:
+        return np.array(lattice_chern)
+    else:
+        return np.array(lattice_chern), gs_projector
\ No newline at end of file
diff --git a/strawberrypy/_pythtb/singlepoint_invariant.py b/strawberrypy/_pythtb/singlepoint_invariant.py
new file mode 100644
index 0000000..5ca9c48
--- /dev/null
+++ b/strawberrypy/_pythtb/singlepoint_invariant.py
@@ -0,0 +1,113 @@
+import numpy as np
+import scipy.linalg as la
+from .lattice import *
+from strawberrypy.common_func import *
+
+def single_point_chern (model, formula='symmetric'):
+    n_orb = model.get_num_orbitals()   
+    n_occ = n_orb//2
+
+    # Gamma point is where to performe the single diagonalization
+    point=[0.,0.]
+
+    chern = []
+    _, u_n0 = model.solve_one(point, eig_vectors=True)
+
+    b1, b2 = reciprocal_vec(model)
+    u_nb1 = periodic_gauge(u_n0, b1, model)
+    u_nb2 = periodic_gauge(u_n0, b2, model)
+
+    udual_nb1 = dual_state(u_n0, u_nb1, n_occ, n_orb)
+    udual_nb2 = dual_state(u_n0, u_nb2, n_occ, n_orb)
+
+    if (formula=='asymmetric' or formula =='both'):
+        sum_occ = 0.
+        for i in range(n_occ):
+            sum_occ += np.vdot(udual_nb1[i],udual_nb2[i])
+
+        chern.append(-np.imag(sum_occ)/np.pi)
+
+    if (formula=='symmetric' or formula =='both'):
+        u_nmb1 = periodic_gauge(u_n0, -b1, model)
+        u_nmb2 = periodic_gauge(u_n0, -b2, model)
+
+        udual_nmb1 = dual_state(u_n0, u_nmb1, n_occ, n_orb)
+        udual_nmb2 = dual_state(u_n0, u_nmb2, n_occ, n_orb)
+
+        sum_occ = 0.
+        for i in range(n_occ):
+            sum_occ += np.vdot((udual_nmb1[i]-udual_nb1[i]),(udual_nmb2[i]-udual_nb2[i]))
+
+        chern.append(-np.imag(sum_occ)/(4*np.pi) )
+
+    return chern
+
+def single_point_spin_chern(model, spin='down', formula='symmetric', return_gap_pszp=False):
+    n_orb = model.get_num_orbitals()
+    n_occ = n_orb
+    n_bande = 2*n_occ
+    n_sub = n_occ//2
+
+    # Gamma point is where to performe the single diagonalization
+    point=[0.,0.]
+
+    closing_gap = False
+    spin_chern = []
+
+    b1, b2 = reciprocal_vec(model)
+
+    _, u_0 = model.solve_one(point, eig_vectors=True)
+
+    u_n0 = np.zeros([n_bande, n_bande], dtype = np.complex128)
+    for n in range (n_bande):
+        u_n0[n,:] = u_0[n].reshape((1, -1))
+
+    pszp = pszp_matrix(u_n0, n_occ)
+    eval_pszp, eig_pszp = la.eigh(pszp)
+    gap_pszp = eval_pszp[n_sub] - eval_pszp[n_sub-1]
+    eig_pszp = eig_pszp.T
+
+    if gap_pszp<10**(-14):
+        closing_gap = True
+
+    if closing_gap==True:
+        raise Exception('Closing PszP gap!!')
+    else:
+        #check symmetry of P Sz P spectrum 
+        if (eval_pszp[n_sub]*eval_pszp[n_sub-1]>0):
+            raise Exception('P Sz P spectrum NOT symmetric!!!')
+        
+        q_0 = np.zeros([n_occ, n_bande], dtype=complex)
+        q_0 = eig_pszp @ u_n0[:n_occ,:]
+
+        q_b1 = periodic_gauge_spin(q_0, b1, model)
+        q_b2 = periodic_gauge_spin(q_0, b2, model)
+
+        qdual_b1 = dual_state_spin(q_0, q_b1, n_sub, n_bande, spin)
+        qdual_b2 = dual_state_spin(q_0, q_b2, n_sub, n_bande, spin)
+
+        if (formula=='asymmetric' or formula=='both'):
+            sum_occ = 0.
+            for i in range(n_sub):
+                sum_occ += np.vdot(qdual_b1[i],qdual_b2[i])
+
+            spin_chern.append(-np.imag(sum_occ)/np.pi)
+
+
+        if (formula=='symmetric' or formula=='both'):
+            q_mb1 = periodic_gauge_spin(q_0,-b1,model)
+            q_mb2 = periodic_gauge_spin(q_0,-b2,model)
+
+            qdual_mb1 = dual_state_spin(q_0, q_mb1, n_sub, n_bande, spin)
+            qdual_mb2 = dual_state_spin(q_0, q_mb2, n_sub, n_bande, spin)
+
+            sum_occ = 0.
+            for i in range(n_sub):
+                sum_occ += np.vdot((qdual_mb1[i]-qdual_b1[i]),(qdual_mb2[i]-qdual_b2[i]))
+
+            spin_chern.append(-np.imag(sum_occ)/(4*np.pi)) 
+
+    if return_gap_pszp==True:
+        spin_chern.append(gap_pszp)
+    
+    return spin_chern
diff --git a/strawberrypy/_tbmodels/__init__.py b/strawberrypy/_tbmodels/__init__.py
new file mode 100644
index 0000000..7359c02
--- /dev/null
+++ b/strawberrypy/_tbmodels/__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
\ No newline at end of file
diff --git a/strawberrypy/_tbmodels/finite_systems.py b/strawberrypy/_tbmodels/finite_systems.py
new file mode 100644
index 0000000..1937e02
--- /dev/null
+++ b/strawberrypy/_tbmodels/finite_systems.py
@@ -0,0 +1,218 @@
+import numpy as np
+import tbmodels as tbm
+
+def get_positions(model, nx_sites : int, ny_sites : int):
+    """
+    Returns the cartesian coordinates of the states in a finite sample of tbmodels.Model
+    """
+    
+    positions = np.copy(model.pos)
+    for i in range(model.pos.shape[0]):
+        positions[i][0] *= nx_sites
+        positions[i][1] *= ny_sites
+        cartesian_pos = np.dot(positions[i], model.uc)
+        positions[i] = cartesian_pos
+    return positions
+
+def cut_piece_tbm(source_model, num : int, fin_dir : int, dimk : int, glue : bool = False):
+    """
+    Reduce by 1 the dimension of a model in TBModels in a given direction, building a supercell made by a given number of atoms. If a zero dimensional system, is needed it is required to run twice the function with dimk = 1 and then dimk = 0.
+    
+        Args:
+        - source_model: the model to reduce
+        - num: number of cell to keep in the finite direction
+        - fin_dir: finite direction (ie, 0 = x, 1 = y)
+        - dimk: number of periodic directions after the cut
+        - glue: glue edges to make periodicity
+
+        Returns:
+        - tbmodels.Model finite along fin_dir
+    """
+
+    # Ensure the condition the function has been built for
+    if num <= 0:
+        raise RuntimeError("Negative number of cells in the finite direction required.")
+    if fin_dir not in [0, 1]:
+        raise RuntimeError("Finite direction not allowed (only 2D systems).")
+    if dimk not in [0, 1]:
+        raise RuntimeError("Leftover k-space dimension not allowed.")
+    if num == 1 and glue == True:
+        raise RuntimeError("Cannot glue edges with one cell in the finite direction.")
+    
+    # Number of orbitals in the supercell model = norbs (original) x num
+    norbs = source_model.size
+
+    # Define the supercell
+    newpos = []
+    for i in range(num):
+        for j in range(norbs):
+            # Convert coordinates into cartesian coordinates
+            orb_tmp = np.copy(source_model.pos[j, :])
+
+            # One direction is fine but I need to map the other into the unit cell
+            orb_tmp[fin_dir] += float(i)
+            orb_tmp[fin_dir] /= num
+
+            # Convert the coordinates back into lattice coordinates
+            newpos.append(orb_tmp)
+
+    # Onsite energies per unit cell (2 is by convention with TBModels)
+    onsite = num * [2 * np.real(source_model.hop[source_model._zero_vec][j][j]) for j in range(norbs)]
+
+    # Hopping amplitudes and positions
+    hoppings = [[key, val] for key, val in iter(source_model.hop.items())]
+
+    # Hoppings to be added
+    hopping_list = []
+
+    # Cycle over the number of defined hoppings
+    for j in range(len(hoppings)):
+
+        # Set lattice vector of the current hopping matrix
+        objective = np.copy(hoppings[j][0])
+
+        # Maximum bond length
+        jump_fin = hoppings[j][0][fin_dir]
+
+        # If I have a finite direction I make the hopping vector finite, and
+        #   if I have no periodic direction, I but every hopping in the [zero] cell
+        if dimk != 0:
+            objective[fin_dir] = 0
+        else:
+            objective = np.array([0 for i in range(source_model.dim)])
+
+        # Cycle over the rows of the hopping matrix
+        for k in range(hoppings[j][1].shape[0]):
+
+            # Cycle over the columns of the hopping matrix
+            for l in range(hoppings[j][1].shape[1]):
+
+                # Hopping amplitudes
+                amplitude = hoppings[j][1][k][l]
+                if np.absolute(amplitude) < 1e-10:
+                    continue
+
+                # Cycle over the cells in the supercell
+                for i in range(num):
+                    starting = k + i * norbs
+                    ending = l + (i + jump_fin) * norbs
+
+                    # Decide wether to add the hopping or not
+                    to_add = True
+
+                    if not glue:
+                        if ending < 0 or ending >= norbs * num:
+                            to_add = False
+                    else:
+                        ending = int(ending) % int(norbs * num)
+
+                    # Avoid setting twice onsite energies
+                    if starting == ending and (objective == [0 for i in range(source_model.dim)]).all():
+                        continue
+
+                    if to_add == True:
+                        hopping_list.append([amplitude, int(starting), int(ending), objective])
+    
+    model = tbm.Model.from_hop_list(hop_list = hopping_list, on_site = onsite, size = norbs * num, dim = source_model.dim,
+        occ = source_model.occ * num, uc = source_model.uc, pos = newpos, contains_cc = False)
+    
+    return model
+
+def make_finite(model: tbm.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 = cut_piece_tbm(model, num = ny_sites, fin_dir = 1, dimk = 1, glue = False)
+    finite = cut_piece_tbm(ribbon, num = nx_sites, fin_dir = 0, dimk = 0, glue = False)
+    
+    return finite
+
+def make_heterostructure(model1 : tbm.Model , model2 : tbm.Model,  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.size / (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.size == model2.size and np.allclose(model1.pos, model2.pos) and np.allclose(model1.uc, model2.uc) and model1.occ == model2.occ):
+        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.diag(list(model1.hop.values())[0]).copy()
+    onsite2 = np.diag(list(model2.hop.values())[0]).copy()
+    
+    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]
+    onsite1 = [2 * term for term in onsite1]
+
+    # 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 = [[key, val] for key, val in iter(model1.hop.items())][0][1]
+    hoppings2 = [[key, val] for key, val in iter(model2.hop.items())][0][1]
+    hopping_list = []
+
+    # Cycle over the rows of the hopping matrix
+    for k in range(hoppings1.shape[0]):
+
+        # Cycle over the columns of the hopping matrix
+        for l in range(hoppings1.shape[1]):
+            if k == l: continue
+
+            # Hopping amplitudes
+            amplitude1 = hoppings1[k][l]
+            amplitude2 = hoppings2[k][l]
+                    
+            if k in indices:
+                if np.absolute(amplitude2) < 1e-10: continue
+                hopping_list.append([amplitude2, k, l, [0 for i in range(model1.dim)]])
+            else:
+                if np.absolute(amplitude1) < 1e-10: continue
+                hopping_list.append([amplitude1, k, l, [0 for i in range(model1.dim)]])
+
+    newmodel = tbm.Model.from_hop_list(hop_list = hopping_list, on_site = onsite1, size = model1.size, dim = model1.dim,
+        occ = model1.occ, uc = model1.uc, pos = model1.pos, contains_cc = False)
+
+    return newmodel
\ No newline at end of file
diff --git a/strawberrypy/_tbmodels/lattice.py b/strawberrypy/_tbmodels/lattice.py
new file mode 100644
index 0000000..7b370f6
--- /dev/null
+++ b/strawberrypy/_tbmodels/lattice.py
@@ -0,0 +1,96 @@
+import numpy as np
+import tbmodels
+
+def reciprocal_vec(model):
+    #returns reciprocal lattice vectors
+    b_matrix = model.reciprocal_lattice
+    b1 = b_matrix[0,:]
+    b2 = b_matrix[1,:]
+    return b1, b2
+
+def orb_cart(model):
+    #returns position of orbitals in cartesian coordinates
+    orb_red = model.pos
+    n_orb = len(model.pos)
+    lat_super = model.uc
+
+    orb_c = []
+    for i in range (n_orb):
+        orb_c.append((np.matmul(lat_super.transpose(),orb_red[i].reshape(-1,1))).squeeze())
+    orb_c = np.array(orb_c)
+    return orb_c
+
+def periodic_gauge(u_n0, b, model):
+    #returns the matrix of occupied eigenvectors at the edge of the Brillouin zone imposing periodic gauge
+    n_occ = model.occ
+
+    orb_c = orb_cart(model)
+    vec_scal_b = orb_c @ b
+    vec_exp_b = np.exp(-1.j*vec_scal_b)
+
+    u_nb = vec_exp_b.T * u_n0[:n_occ,:]
+
+    return u_nb
+
+def onsite_disorder(source_model, w : float = 0, spinstates : int = 2, seed : int = None):
+    """
+    Add onsite (Anderson) disorder to the specified model. The disorder amplitude per site is taken randomly in [-w/2, w/2].
+
+        Args:
+        - source_model : the model to add disorder to
+        - w : disorder amplitude
+        - spinstates : spin of the model
+        - seed : seed for random number generator
+
+        Returns:
+        - model : the disordered model
+    """
+    # Quick return for no disorder
+    if w == 0:
+        return source_model
+
+    # Set the seed for the random number generator
+    if seed is not None:
+        np.random.seed(seed)
+
+    # Number of orbitals in the supercell model = norbs (original) x num
+    norbs = source_model.size
+    
+    # Onsite energies per unit cell (2 is by convention with TBModels)
+    onsite = [2 * np.real(source_model.hop[source_model._zero_vec][j][j]) for j in range(norbs)]
+    disorder = 0.5 * w * (2 * np.random.rand(norbs // spinstates) - 1.0)
+    disorder = np.repeat(disorder, spinstates)
+    onsite += disorder
+
+    # Hopping amplitudes and positions
+    hoppings = [[key, val] for key, val in iter(source_model.hop.items())]
+
+    # Hoppings to be added
+    hopping_list = []
+
+    # Cycle over the number of defined hoppings
+    for j in range(len(hoppings)):
+
+        # Set lattice vector of the current hopping matrix
+        objective = np.copy(hoppings[j][0])
+
+        # Cycle over the rows of the hopping matrix
+        for k in range(hoppings[j][1].shape[0]):
+
+            # Cycle over the columns of the hopping matrix
+            for l in range(hoppings[j][1].shape[1]):
+
+                # Hopping amplitudes
+                amplitude = hoppings[j][1][k][l]
+                if np.absolute(amplitude) < 1e-10:
+                    continue
+
+                if k == l:
+                    continue
+                else:
+                    hopping_list.append([amplitude, int(k), int(l), objective])
+    
+    model = tbmodels.Model.from_hop_list(hop_list = hopping_list, on_site = onsite, size = norbs, dim = source_model.dim,
+        occ = source_model.occ, uc = source_model.uc, pos = source_model.pos, contains_cc = False)
+    
+    return model
\ No newline at end of file
diff --git a/strawberrypy/_tbmodels/localmarkers.py b/strawberrypy/_tbmodels/localmarkers.py
new file mode 100644
index 0000000..fbb6fa7
--- /dev/null
+++ b/strawberrypy/_tbmodels/localmarkers.py
@@ -0,0 +1,93 @@
+import tbmodels as tbm
+import numpy as np
+import scipy.linalg as la
+from opt_einsum import contract
+
+from .finite_systems import *
+from strawberrypy.common_func import *
+
+def local_chern_marker(model, nx_sites : int, ny_sites : int, direction : int = None, start : int = 0, return_projector : bool = False, projector : np.ndarray = None, macroscopic_average : bool = False, cutoff : float = 0.8):
+    """
+    Evaluate the Z Chern marker on the whole lattice if direction is None. If direction is not None evaluates the Z Chern marker along direction starting from start.
+    
+        Args:
+        - model : instance of TBModels of the model
+        - nx_sites : number of unit cells in the x direction of the finite system
+        - ny_sites : number of unit cells in the y direction of the finite system
+        - direction : direction along which compute the local Chern marker, default is None (returns the whole lattice Chern marker), allowed values are 0 for 'x' direction and 1 for 'y' direction
+        - start : if direction is not None, is the coordinate of the unit cell at the start of the evaluation of the Chern marker
+        - return_projector : if True, returns the ground state projector at the end of the calculation, default is False
+        - projector : input the ground state projector to be used in the calculation. Default is None, which means it is computed from the model
+        - macroscopic_average : if True, returns the local Chern marker averaged over a radius equal to the cutoff
+        - cutoff : cutoff set for the calculation of averages
+
+        Returns:
+        - lattice_chern : local Chern marker of the whole lattice if direction is None
+        - lcm_direction : local Chern marker along direction starting from start
+        - projector : ground state projector, returned if return_projector is set True (default is False)
+    """
+
+    # Check input variables
+    if direction not in [None, 0, 1]:
+        raise RuntimeError("Direction allowed are None, 0 (which stands for x), and 1 (which stands for y)")
+    
+    if direction is not None:
+        if direction == 0:
+            if start not in range(ny_sites): raise RuntimeError("Invalid start parameter (must be within [0, ny_sites - 1])")
+        else:
+            if start not in range(nx_sites): raise RuntimeError("Invalid start parameter (must be within [0, nx_sites - 1])")
+
+    # Atoms in the unit cell
+    atoms_uc = int(model.size / (nx_sites * ny_sites))
+
+    # 2D unit cell volume
+    uc_vol = np.linalg.norm(np.cross(model.uc[0], model.uc[1]))
+
+    if projector is None:
+        # Eigenvectors at \Gamma
+        hamiltonian = model.hamilton([0, 0], convention = 1)
+        _, eigenvecs = la.eigh(hamiltonian)
+        eigenvecs = eigenvecs
+
+        # Build the ground state projector
+        gs_projector = contract("ji,ki->jk", eigenvecs[:, :int(0.5 * len(eigenvecs))], eigenvecs[:, :int(0.5 * len(eigenvecs))].conjugate())
+    else:
+        gs_projector = projector
+
+    # Convert TBM-positions in actual positions
+    positions = get_positions(model, nx_sites, ny_sites)
+
+    # Position operator in tight-binding approximation (site orbitals basis)
+    rx = np.diag(positions[:, 0]); ry = np.diag(positions[:, 1])
+
+    # Chern marker operator
+    chern_operator = np.imag(gs_projector @ commutator(rx, gs_projector) @ commutator(ry, gs_projector))
+    chern_operator *= -4 * np.pi / uc_vol
+
+    # If macroscopic average I have to compute the lattice values with the averages first
+    if macroscopic_average:
+        chern_on_lattice = average_over_radius(np.diag(chern_operator), positions.T[0], positions.T[1], cutoff, nx_sites, ny_sites, atoms_uc)
+    
+    if direction is not None:
+        # Evaluate index of the selected direction
+        indices = xy_to_index('x' if direction == 1 else 'y', start, nx_sites, ny_sites, atoms_uc)
+
+        # If macroscopic average consider the averaged lattice, else the Chern operator
+        if macroscopic_average:
+            lcm_direction = [chern_on_lattice[int(indices[atoms_uc * i] / atoms_uc)] for i in range(int(len(indices) / atoms_uc))]
+        else:
+            lcm_direction = [np.sum([chern_operator[indices[atoms_uc * i + j], indices[atoms_uc * i + j]] for j in range(atoms_uc)]) for i in range(int(len(indices) / atoms_uc))]
+        
+        if not return_projector:
+            return np.array(lcm_direction)
+        else:
+            return np.array(lcm_direction), gs_projector
+
+    if not macroscopic_average:
+        chern_on_lattice = partialtrace(chern_operator, atoms_uc)
+
+    lattice_chern = chern_on_lattice.reshape(nx_sites, ny_sites).T
+    if not return_projector:
+        return np.array(lattice_chern)
+    else:
+        return np.array(lattice_chern), gs_projector
\ No newline at end of file
diff --git a/strawberrypy/_tbmodels/singlepoint_invariant.py b/strawberrypy/_tbmodels/singlepoint_invariant.py
new file mode 100644
index 0000000..ea60e7f
--- /dev/null
+++ b/strawberrypy/_tbmodels/singlepoint_invariant.py
@@ -0,0 +1,115 @@
+import numpy as np
+import scipy.linalg as la
+from .lattice import *
+from strawberrypy.common_func import *
+
+def solve_one(model, k):
+    hamiltonian = model.hamilton(k,convention=1)
+
+    eigval, eigvec = la.eigh(hamiltonian)
+
+    return eigval, eigvec.T
+
+def single_point_chern (model, formula='symmetric'):
+    n_orb = len(model.pos)
+    n_occ = model.occ
+
+    # Gamma point is where to performe the single diagonalization
+    point=[0.,0.]
+
+    chern = []
+    _, u_n0 = solve_one(model,point)
+
+    b1, b2 = reciprocal_vec(model)
+    u_nb1 = periodic_gauge(u_n0, b1, model)
+    u_nb2 = periodic_gauge(u_n0, b2, model)
+
+    udual_nb1 = dual_state(u_n0, u_nb1, n_occ, n_orb)
+    udual_nb2 = dual_state(u_n0, u_nb2, n_occ, n_orb)
+
+    if (formula=='asymmetric' or formula =='both'):
+        sum_occ = 0.
+        for i in range(n_occ):
+            sum_occ += np.vdot(udual_nb1[i],udual_nb2[i])
+
+        chern.append(-np.imag(sum_occ)/np.pi)
+
+    if (formula=='symmetric' or formula =='both'):
+        u_nmb1 = periodic_gauge(u_n0, -b1, model)
+        u_nmb2 = periodic_gauge(u_n0, -b2, model)
+
+        udual_nmb1 = dual_state(u_n0, u_nmb1, n_occ, n_orb)
+        udual_nmb2 = dual_state(u_n0, u_nmb2, n_occ, n_orb)
+
+        sum_occ = 0.
+        for i in range(n_occ):
+            sum_occ += np.vdot((udual_nmb1[i]-udual_nb1[i]),(udual_nmb2[i]-udual_nb2[i]))
+
+        chern.append(-np.imag(sum_occ)/(4*np.pi) )
+
+    return chern
+
+def single_point_spin_chern(model, spin='down', formula='symmetric', return_gap_pszp=False):
+    n_orb = len(model.pos)
+    n_occ = model.occ
+    n_sub = n_occ//2
+
+    # Gamma point is where to performe the single diagonalization
+    point=[0.,0.]
+
+    closing_gap = False
+    spin_chern = []
+
+    b1, b2 = reciprocal_vec(model)
+
+    _, u_n0 = solve_one(model,point)
+    
+    pszp = pszp_matrix(u_n0, n_occ)
+    eval_pszp, eig_pszp = la.eigh(pszp)
+    gap_pszp = eval_pszp[n_sub] - eval_pszp[n_sub-1]
+    eig_pszp = eig_pszp.T
+
+    if gap_pszp<10**(-14):
+        closing_gap = True
+
+    if closing_gap==True:
+        raise Exception('Closing PszP gap!!')
+    else:
+        #check symmetry of P Sz P spectrum 
+        if (eval_pszp[n_sub]*eval_pszp[n_sub-1]>0):
+            raise Exception('P Sz P spectrum NOT symmetric!!!')
+
+        q_0 = np.zeros([n_occ, n_orb], dtype=complex)
+        q_0 = eig_pszp @ u_n0[:n_occ,:]
+
+        q_b1 = periodic_gauge(q_0, b1, model)
+        q_b2 = periodic_gauge(q_0, b2, model)
+
+        qdual_b1 = dual_state_spin(q_0, q_b1, n_sub, n_orb, spin)
+        qdual_b2 = dual_state_spin(q_0, q_b2, n_sub, n_orb, spin)
+
+        if (formula=='asymmetric' or formula=='both'):
+            sum_occ = 0.
+            for i in range(n_sub):
+                sum_occ += np.vdot(qdual_b1[i],qdual_b2[i])
+
+            spin_chern.append(-np.imag(sum_occ)/np.pi)
+
+
+        if (formula=='symmetric' or formula=='both'):
+            q_mb1 = periodic_gauge(q_0,-b1,model)
+            q_mb2 = periodic_gauge(q_0,-b2,model)
+
+            qdual_mb1 = dual_state_spin(q_0, q_mb1, n_sub, n_orb, spin)
+            qdual_mb2 = dual_state_spin(q_0, q_mb2, n_sub, n_orb, spin)
+
+            sum_occ = 0.
+            for i in range(n_sub):
+                sum_occ += np.vdot((qdual_mb1[i]-qdual_b1[i]),(qdual_mb2[i]-qdual_b2[i]))
+
+            spin_chern.append(-np.imag(sum_occ)/(4*np.pi)) 
+
+    if return_gap_pszp==True:
+        spin_chern.append(gap_pszp)
+   
+    return spin_chern
diff --git a/strawberrypy/common_func.py b/strawberrypy/common_func.py
new file mode 100644
index 0000000..203a5cf
--- /dev/null
+++ b/strawberrypy/common_func.py
@@ -0,0 +1,121 @@
+import numpy as np
+import random
+import scipy.linalg as la
+import time
+
+def dual_state(un0, unb, n_occ, n_orb):
+    s_matrix_b = np.zeros([n_occ,n_occ], dtype=np.complex128)
+    s_matrix_b = np.conjugate(un0[:n_occ,:]) @ (unb[:n_occ,:]).T
+    
+    s_inv_b = np.linalg.inv(s_matrix_b)
+
+    udual_nb = np.zeros((n_occ, n_orb), dtype=np.complex128)
+    udual_nb = (s_inv_b.T) @ unb[:n_occ,:]
+
+    return udual_nb
+
+def dual_state_spin(q0, qb, n_sub, n_orb, spin):
+
+    s_matrix_b = np.zeros([n_sub,n_sub], dtype=np.complex128)
+    udual_nb = np.zeros((n_sub, n_orb), dtype=np.complex128)
+
+    if spin ==  'down':
+        s_matrix_b = np.conjugate(q0[:n_sub,:]) @ (qb[:n_sub,:]).T
+        s_inv_b = np.linalg.inv(s_matrix_b)
+        udual_nb = (s_inv_b.T) @ qb[:n_sub,:]
+    elif spin == 'up':
+        s_matrix_b = np.conjugate(q0[n_sub:,:]) @ (qb[n_sub:,:]).T
+        s_inv_b = np.linalg.inv(s_matrix_b)
+        udual_nb = (s_inv_b.T) @ qb[n_sub:,:]
+    return udual_nb
+
+def pszp_matrix (u_n0, n_occ):
+    sz = [1,-1]*n_occ
+    sz_un0 = (sz*u_n0).T
+    pszp = np.ndarray([n_occ,n_occ],dtype=complex)
+    pszp = u_n0[:n_occ,:].conjugate() @ sz_un0[:,:n_occ]
+    return pszp
+
+def commutator(A, B):
+    """
+    Computes the commutator [A, B]
+    """
+    return np.array(A @ B - B @ A)
+
+def xy_to_index(fixed_coordinate : str, xy : int, nx : int, ny : int, atoms_uc : int):
+    """
+    Returns the indices of all the sites (including indices of atoms within the untit cell) along fixed_coordinate starting from position xy
+
+        Args:
+        - fixed_coordinate : direction in which look for the indices, allowed values are 'x' and 'y'
+        - xy : starting position
+        - nx, ny : lattice dimension
+        - atoms_uc : number of atoms in the unit cell
+
+        Returns:
+        - indices : the list of indiecs along the given direction
+    """
+    if fixed_coordinate == 'x':
+        return np.array([atoms_uc * xy * ny + i for i in range(atoms_uc * ny)]).flatten().tolist()
+    elif fixed_coordinate == 'y':
+        return np.array([[atoms_uc * xy + atoms_uc * ny * i + j for j in range(atoms_uc)] for i in range(0, nx)]).flatten().tolist()
+    else:
+        raise RuntimeError("Direction not allowed, only 'x' or 'y'")
+    
+def partialtrace(matrix, num : int):
+    """
+    Computes the partial trace of a given matrix, num specifies the number of elements to include in the partial trace
+    """
+    if not int(len(matrix) / num) == len(matrix) / num:
+        raise RuntimeError("Number of element to take the partial trace does not divide the matrix dimension.")
+    return np.array([np.sum([matrix[i * num + k, i * num + k] for k in range(num)]) for i in range(int(len(matrix) / num))])
+
+def lattice_contraction(xcoord : np.ndarray, ycoord : np.ndarray, cutoff : float, nx : int, ny : int, atoms_uc : int):
+    """
+    Define the sites in the lattice that have to be contracted for each point in order to evaluate averages
+
+        Args:
+        - xcoord, ycoord : coordinates of the sites
+        - cutoff : maximum radius of the contraction window
+        - nx, ny : dimensions of the lattice
+        - atoms_uc : atoms in the unit cell
+
+        Returns:
+        - contraction : the list of sites that have to be contracted in a point to evaluate averages for each point
+    """
+    contraction = []
+    def within_range(current, trial):
+        return True if (xcoord[current] - xcoord[trial]) ** 2 + (ycoord[current] - ycoord[trial]) ** 2 < cutoff ** 2 else False
+
+    for current in range(nx * ny):
+        contraction.append(np.array([np.array([atoms_uc * trial + ind for ind in range(atoms_uc)], dtype = int) for trial in range(nx * ny) if within_range(atoms_uc * trial, atoms_uc * current)]).flatten())
+
+    return contraction
+
+def average_over_radius(vals, xcoord, ycoord, cutoff, nx, ny, atoms_uc, contraction = None):
+    """
+    Computes the average over a certain radius of vals in a lattice
+
+        Args:
+        - vals : the values to average
+        - xoord, ycoord : coordinates of the sites
+        - cutoff : maximum radius of the contraction window
+        - nx, ny : dimensions of the lattice
+        - atoms_uc : atoms in the unit cell
+        - contraction : if the contraction window has already been calculated, it is possible to pass it, default is None and the contraction window is evaluated here
+
+        Returns:
+        - contracted_vals : values averaged using the contraction
+    """
+    return_vals = []
+
+    if contraction is None:
+        contraction = lattice_contraction(xcoord, ycoord, cutoff, nx, ny, atoms_uc)
+
+    # Macroscopia average within a certain radius
+    for current in range(nx * ny):
+        tmp = [vals[ind] for ind in contraction[current]]
+        return_vals.append(np.sum(tmp) / (len(tmp) / atoms_uc))
+
+    # Sum the averages over the unit cell
+    return np.array(return_vals)
\ No newline at end of file
diff --git a/strawberrypy/example_models/__init__.py b/strawberrypy/example_models/__init__.py
new file mode 100644
index 0000000..5bf94b9
--- /dev/null
+++ b/strawberrypy/example_models/__init__.py
@@ -0,0 +1,2 @@
+from .kane_mele import kane_mele_pythtb, kane_mele_tbmodels, km_anderson_disorder_pythtb, km_anderson_disorder_tbmodels
+from .haldane import haldane_pythtb, haldane_tbmodels, h_anderson_disorder_pythtb, h_anderson_disorder_tbmodels
\ No newline at end of file
diff --git a/strawberrypy/example_models/haldane.py b/strawberrypy/example_models/haldane.py
new file mode 100644
index 0000000..1ae2f08
--- /dev/null
+++ b/strawberrypy/example_models/haldane.py
@@ -0,0 +1,53 @@
+from pythtb import tb_model
+from tbmodels import Model
+
+import numpy as np
+
+def haldane_pythtb(delta, t, t2, phi, L):
+    # From http://www.physics.rutgers.edu/pythtb/examples.html#haldane-model
+    lat=[[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
+    orb=[[0.0,0.0],[1./3.,1./3.]]
+    model=tb_model(2,2,lat,orb)
+    model.set_onsite([-delta,delta])
+    for lvec in ([ 0, 0], [-1, 0], [ 0,-1]):
+        model.set_hop(t, 0, 1, lvec)
+    for lvec in ([ 1, 0], [-1, 1], [ 0,-1]):
+        model.set_hop(t2*np.exp(1.j*phi), 0, 0, lvec)
+    for lvec in ([-1, 0], [ 1,-1], [ 0, 1]):
+        model.set_hop(t2*np.exp(1.j*phi), 1, 1, lvec)
+
+    sc_model = model.make_supercell([[L,0],[0,L]])
+    return sc_model
+
+def haldane_tbmodels(delta, t, t2, phi, L):
+    primitive_cell = [[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
+    orb = [[0.0,0.0],[1./3.,1./3.]]
+    h_model = Model(on_site=[-delta,delta], dim=2, occ=1, pos=orb, uc=primitive_cell)
+    for lvec in ([0,0],[-1,0],[0,-1]):
+        h_model.add_hop(t, 0,1,lvec)
+    for lvec in ([1,0],[-1,1],[0,-1]):
+        h_model.add_hop(t2*np.exp(1.j*phi), 0,0,lvec)
+    for lvec in ([-1,0],[1,-1],[0,1]):
+        h_model.add_hop(t2*np.exp(1.j*phi), 1,1,lvec) 
+
+    sc_model = h_model.supercell([L,L])
+       
+    return sc_model
+
+def h_anderson_disorder_pythtb(model, w):
+    n_orb = model.get_num_orbitals() 
+    if w != 0.:
+        for j in range (n_orb):
+            dis = 0.5*w*(2*np.random.random()-1.0)
+            model.set_onsite(dis, j, mode='add')
+
+    return model
+
+def h_anderson_disorder_tbmodels(model, w):
+    n_orb = len(model.pos)
+    if w != 0.:
+        dis = 0.5*w*(2*np.random.rand(n_orb)-1.0)                   
+
+        model.add_on_site(dis)
+
+    return model
\ No newline at end of file
diff --git a/strawberrypy/example_models/kane_mele.py b/strawberrypy/example_models/kane_mele.py
new file mode 100644
index 0000000..c86bfcf
--- /dev/null
+++ b/strawberrypy/example_models/kane_mele.py
@@ -0,0 +1,110 @@
+from pythtb import tb_model
+from tbmodels import Model
+
+import numpy as np
+
+def kane_mele_pythtb(rashba, esite, spin_orb, L):
+    # From http://www.physics.rutgers.edu/pythtb/examples.html#kane-mele-model-using-spinor-features
+
+    # define lattice vectors
+    lat=[[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
+    # define coordinates of orbitals
+    orb=[[0.0,0.0],[1./3.,1./3.]]
+
+    # make two dimensional tight-binding Kane-Mele model
+    km_model=tb_model(2,2,lat,orb,nspin=2)
+
+    # set other parameters of the model
+    thop=1.0
+
+    rashba = rashba*spin_orb
+    esite = esite*spin_orb
+
+
+    # set on-site energies
+    km_model.set_onsite([esite, -esite])
+
+    # set hoppings (one for each connected pair of orbitals)
+    # (amplitude, i, j, [lattice vector to cell containing j])
+
+    # useful definitions
+    sigma_x=np.array([0.,1.,0.,0])
+    sigma_y=np.array([0.,0.,1.,0])
+    sigma_z=np.array([0.,0.,0.,1])
+
+    # spin-independent first-neighbor hops
+    for lvec in ([ 0, 0], [-1, 0], [ 0,-1]):
+        km_model.set_hop(thop, 0, 1, lvec)
+        
+    # spin-dependent second-neighbor hops
+    for lvec in ([ 1, 0], [-1, 1], [ 0,-1]):
+       km_model.set_hop(1.j*spin_orb*sigma_z, 0, 0, lvec)
+    for lvec in ([-1, 0], [ 1,-1], [ 0, 1]):
+       km_model.set_hop(1.j*spin_orb*sigma_z, 1, 1, lvec)
+
+    # Rashba first-neighbor hoppings: (s_x)(dy)-(s_y)(d_x)
+    r3h =np.sqrt(3.0)/2.0
+    # bond unit vectors are (r3h,half) then (0,-1) then (-r3h,half)
+    km_model.set_hop(1.j*rashba*( 0.5*sigma_x-r3h*sigma_y), 0, 1, [ 0, 0], mode="add")
+    km_model.set_hop(1.j*rashba*(-1.0*sigma_x            ), 0, 1, [ 0,-1], mode="add")
+    km_model.set_hop(1.j*rashba*( 0.5*sigma_x+r3h*sigma_y), 0, 1, [-1, 0], mode="add")
+
+    sc_model = km_model.make_supercell([[L,0],[0,L]])
+
+    return sc_model
+   
+def kane_mele_tbmodels(rashba,esite,spin_orb,L):
+
+    t_hop = 1.0
+    r = rashba*spin_orb
+    e = esite*spin_orb
+
+    primitive_cell = [[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
+    orb = [[0.0,0.0],[0.0,0.0],[1./3.,1./3.],[1./3.,1./3.]]
+    #set energy on site
+    km_model = Model(on_site=[e,e,-e,-e], dim=2, occ=2, pos=orb, uc=primitive_cell)
+
+    #spin-independent first-neighbor hops
+    for lvec in ([0,0],[-1,0],[0,-1]):
+        km_model.add_hop(t_hop,0,2,lvec)
+        km_model.add_hop(t_hop,1,3,lvec)
+
+    #spin-dependent second-neighbor hops
+    for lvec in ([1,0],[-1,1],[0,-1]):
+        km_model.add_hop(1.j*spin_orb, 0,0,lvec)
+        km_model.add_hop(-1.j*spin_orb, 1,1,lvec)
+    for lvec in ([-1,0],[1,-1],[0,1]):
+        km_model.add_hop(1.j*spin_orb,2,2,lvec)
+        km_model.add_hop(-1.j*spin_orb, 3,3,lvec)
+
+    #Rashba first neighbor hoppings : (s_x)(dy)-(s_y)(dx)
+    r3h = np.sqrt(3.0)/2.0
+    #bond unit vectors are (r3h,half) then (0,-1) then (-r3h,half)
+    km_model.add_hop(1.j*r*(0.5-1.j*r3h), 1, 2, [ 0, 0])
+    km_model.add_hop(1.j*r*(0.5+1.j*r3h), 0, 3, [ 0, 0])
+    km_model.add_hop(-1.j*r, 1, 2, [ 0, -1])
+    km_model.add_hop(-1.j*r, 0, 3, [ 0, -1])
+    km_model.add_hop(1.j*r*(0.5+1.j*r3h), 1, 2, [ -1, 0])
+    km_model.add_hop(1.j*r*(0.5-1.j*r3h), 0, 3, [ -1, 0])
+
+    sc_model = km_model.supercell([L,L])
+    return sc_model
+
+def km_anderson_disorder_pythtb(model, w):
+    n_orb = model.get_num_orbitals()   
+    if w != 0.:
+        for j in range (n_orb):
+            dis = 0.5*w*(2*np.random.random()-1.0)
+            model.set_onsite(dis, j, mode='add')
+    
+    return model
+
+def km_anderson_disorder_tbmodels(model, w):
+    n_orb = len(model.pos)
+    if w != 0.:
+        d = 0.5*w*(2*np.random.rand(n_orb//2)-1.0)                   
+        dis = np.repeat(d,2)
+
+        model.add_on_site(dis)
+
+    return model
diff --git a/strawberrypy/package.py b/strawberrypy/package.py
new file mode 100644
index 0000000..8bfb78b
--- /dev/null
+++ b/strawberrypy/package.py
@@ -0,0 +1,113 @@
+from . import _pythtb
+from . import _tbmodels
+
+import pythtb 
+import tbmodels
+
+def single_point_chern(model, formula='symmetric'):
+    if isinstance(model, pythtb.tb_model):
+        return _pythtb.single_point_chern(model, formula)
+    elif isinstance(model, tbmodels.Model):
+        return _tbmodels.single_point_chern(model, formula)
+    else:
+        raise NotImplementedError('Invalid model.')
+
+def single_point_spin_chern(model, spin='down', formula='symmetric', return_gap_pszp=False):
+    if isinstance(model, pythtb.tb_model):
+        return _pythtb.single_point_spin_chern(model, spin, formula, return_gap_pszp)
+    elif isinstance(model, tbmodels.Model):
+        return _tbmodels.single_point_spin_chern(model, spin, formula, return_gap_pszp)
+    else:
+        raise NotImplementedError('Invalid model.')
+    
+# Various lattice functions
+def onsite_disorder(model, w : float, spinstates : int = 2, seed : int = None):
+    """
+    Add an Anderson onsite term in the hamiltonian, randomly distributed in [-w/2, w/2].
+
+        Args:
+        - model : the model to modify
+        - w : maximum amplitude of disorder
+        - spinstates : spin of the model
+        - seed : seed for random number generator
+
+        Returns:
+        - model : the model with added onsite disorder
+    """
+    if isinstance(model, tbmodels.Model):
+        return _tbmodels.onsite_disorder(model, w, spinstates, seed)
+    elif isinstance(model, pythtb.tb_model):
+        return _pythtb.onsite_disorder(model, w, spinstates, seed)
+    else:
+        raise NotImplementedError("Invalid model.")
+    
+# Building finite systems
+    
+def make_finite(model, nx_sites : int = 1, ny_sites : int = 1):
+    """
+    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 isinstance(model, tbmodels.Model):
+        return _tbmodels.make_finite(model, nx_sites, ny_sites)
+    elif isinstance(model, pythtb.tb_model):
+        return _pythtb.make_finite(model, nx_sites, ny_sites)
+    else:
+        raise NotImplementedError("Invalid model.")
+    
+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
+    """
+    if isinstance(model1, tbmodels.Model) and isinstance(model2, tbmodels.Model):
+        return _tbmodels.make_heterostructure(model1 , model2,  nx_sites, ny_sites, direction, start, stop)
+    elif isinstance(model1, pythtb.tb_model) and isinstance(model2, pythtb.tb_model):
+        return _pythtb.make_heterostructure(model1 , model2,  nx_sites, ny_sites, direction, start, stop)
+    else:
+        raise NotImplementedError('Invalid model.')
+
+# Local marker functions
+
+def local_chern_marker(model, nx_sites : int, ny_sites : int, direction : int = None, start : int = 0, return_projector : bool = None, projector = None, macroscopic_average : bool = False, cutoff : float = 0.8):
+    """
+    Evaluate the Z Chern marker on the whole lattice if direction is None. If direction is not None evaluates the Z Chern marker along direction starting from start.
+    
+        Args:
+        - model : instance of the model
+        - nx_sites : number of unit cells in the x direction of the finite system
+        - ny_sites : number of unit cells in the y direction of the finite system
+        - direction : direction along which compute the local Chern marker, default is None (returns the whole lattice Chern marker), allowed values are 0 for 'x' direction and 1 for 'y' direction
+        - start : if direction is not None, is the coordinate of the unit cell at the start of the evaluation of the Chern marker
+        - return_projector : if True, returns the ground state projector at the end of the calculation, default is False
+        - projector : input the ground state projector to be used in the calculation. Default is None, which means it is computed from the model
+        - macroscopic_average : if True, returns the local Chern marker averaged over a radius equal to the cutoff
+        - cutoff : cutoff set for the calculation of averages
+
+        Returns:
+        - lattice_chern : local Chern marker of the whole lattice if direction is None
+        - lcm_direction : local Chern marker along direction starting from start
+        - projector : ground state projector, returned if return_projector is set True (default is False)
+    """
+
+    if isinstance(model, tbmodels.Model):
+        return _tbmodels.local_chern_marker(model, nx_sites, ny_sites, direction, start, return_projector, projector, macroscopic_average, cutoff)
+    elif isinstance(model, pythtb.tb_model):
+        return _pythtb.local_chern_marker(model, nx_sites, ny_sites, direction, start, return_projector, projector, macroscopic_average, cutoff)
+    else:
+        raise NotImplementedError('Invalid model.')
\ No newline at end of file
diff --git a/tests/test_h.py b/tests/test_h.py
new file mode 100644
index 0000000..66ece68
--- /dev/null
+++ b/tests/test_h.py
@@ -0,0 +1,49 @@
+import numpy as np
+import math
+
+import os, sys
+currentdir = os.path.dirname(os.path.realpath(__file__))
+parentdir = os.path.dirname(currentdir)
+sys.path.append(parentdir)
+
+from strawberrypy import single_point_chern
+
+from strawberrypy.example_models import haldane_pythtb, haldane_tbmodels, h_anderson_disorder_pythtb, h_anderson_disorder_tbmodels
+
+def test_spcn(L=6, t=-4., t2=1., delta=2., pi_phi=-2., w=1.5, which_formula = 'symmetric'):
+#inputs are:    linear size of supercell LxL
+#               t = first neighbours real hopping
+#               t2 = second neighbours
+#               delta = energy on site
+#               pi_phi --> phi = pi/(pi_phi)  where phi = second neighbours hopping phase
+#               w = disorder stregth W/t
+#               which_formula = choice of single point formula 'asymmetric', 'symmetric' or 'both'     
+
+    #Haldane model parameters
+    phi = np.pi/pi_phi
+
+    #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)
+
+    #add Anderson disorder in PythTB model
+    np.random.seed(15)
+    h_pythtb_model = h_anderson_disorder_pythtb(h_pythtb_model, w)
+
+    #add Anderson disorder in TBmodels model
+    np.random.seed(15)
+    h_tbmodels_model = h_anderson_disorder_tbmodels(h_tbmodels_model, w)
+
+    # Single Point Chern Number (SPCN) calculation for models created with both packages, for the same disorder configuration
+
+    chern_pythtb = single_point_chern(h_pythtb_model, formula=which_formula)
+
+    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 )
+
+    assert math.isclose(chern_pythtb[0],chern_tbmodels[0],abs_tol=1e-10)
\ No newline at end of file
diff --git a/tests/test_km.py b/tests/test_km.py
new file mode 100644
index 0000000..3690b2c
--- /dev/null
+++ b/tests/test_km.py
@@ -0,0 +1,46 @@
+import numpy as np
+import math
+
+import os, sys
+currentdir = os.path.dirname(os.path.realpath(__file__))
+parentdir = os.path.dirname(currentdir)
+sys.path.append(parentdir)
+
+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
+
+def test_spscn(L=6, r=1., e=3., spin_o=0.3, w=2., spin_chern='down', which_formula='symmetric'):
+#inputs are:    linear size of supercell LxL
+#               r = rashba/spin_orb
+#               e = e_onsite/spin_orb
+#               w = disorder stregth W/t
+#               spin_chern = choice of spin Chern number  'up' or 'down' 
+#               which_formula = choice of single point formula 'asymmetric', 'symmetric' or 'both'        
+
+    #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)
+
+    #add Anderson disorder in PythTB model
+    np.random.seed(15)
+    km_pythtb_model = km_anderson_disorder_pythtb(km_pythtb_model,w)
+
+    #add Anderson disorder in TBmodels model
+    np.random.seed(15)
+    km_tbmodels_model = km_anderson_disorder_tbmodels(km_tbmodels_model, w)
+
+    # Single Point Spin Chern Number (SPSCN) calculation for models created with both packages, for the same disorder configuration
+
+    spin_chern_pythtb = single_point_spin_chern(km_pythtb_model, spin=spin_chern, formula=which_formula)
+
+    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 )
+
+    assert math.isclose(spin_chern_pythtb[0],spin_chern_tbmodels[0],abs_tol=1e-10)
+  
diff --git a/tests/test_lcm_h.py b/tests/test_lcm_h.py
new file mode 100644
index 0000000..4dbd15f
--- /dev/null
+++ b/tests/test_lcm_h.py
@@ -0,0 +1,37 @@
+import numpy as np
+import os, sys
+currentdir = os.path.dirname(os.path.realpath(__file__))
+parentdir = os.path.dirname(currentdir)
+sys.path.append(parentdir)
+
+from strawberrypy import local_chern_marker, make_finite, onsite_disorder
+from strawberrypy.example_models import haldane_pythtb, haldane_tbmodels
+
+def test_lcm_haldane():
+    # Construction of the Haldane model with TBmodels
+    hmodel_tbm = haldane_tbmodels(delta = 0.5, t = -1, t2 = 0.15, phi = np.pi / 2, L = 1)
+    hmodel_tbm = make_finite(model = hmodel_tbm, nx_sites = 10, ny_sites = 10)
+
+    # Evaluation of the local Chern marker along the x (0) direction at fixed y = 5
+    lcm_tbm = local_chern_marker(model = hmodel_tbm, nx_sites = 10, ny_sites = 10, direction = 0, start = 5)
+
+    # Construction of the Haldane model with PythTB
+    hmodel_pythtb = haldane_pythtb(delta = 0.5, t = -1, t2 = 0.15, phi = np.pi / 2, L = 1)
+    hmodel_pythtb = make_finite(model = hmodel_pythtb, nx_sites = 10, ny_sites = 10)
+
+    # Evaluation of the local Chern marker along the x (0) direction at fixed y = 5
+    lcm_pythtb = local_chern_marker(model = hmodel_pythtb, nx_sites = 10, ny_sites = 10, direction = 0, start = 5)
+
+    # Assert test for the two results
+    assert np.allclose(lcm_pythtb, lcm_tbm)
+
+    # Add disorder to the models
+    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)
+
+    # Evaluation of the local Chern marker along the x (0) direction at fixed y = 5, averaging over a region of radius 2
+    lcm_tbm = 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)
+
+    # Assert test for the two results
+    assert np.allclose(lcm_pythtb, lcm_tbm)
\ No newline at end of file