From c6fadb383100e70c8cc21f22ec52573b2ade082d Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Wed, 21 Feb 2024 17:10:18 +0100
Subject: [PATCH 01/67] md test

---
 .../extracted_declarations.h                  | 194 ++++++++++++++++++
 python/tutorial.md                            |  29 +++
 2 files changed, 223 insertions(+)
 create mode 100644 python/lattice_symmetries/extracted_declarations.h
 create mode 100644 python/tutorial.md

diff --git a/python/lattice_symmetries/extracted_declarations.h b/python/lattice_symmetries/extracted_declarations.h
new file mode 100644
index 0000000..170a8ab
--- /dev/null
+++ b/python/lattice_symmetries/extracted_declarations.h
@@ -0,0 +1,194 @@
+typedef struct {
+  void *elts;
+  uint64_t num_elts;
+
+  void *freer;
+} chpl_external_array;
+typedef struct ls_hs_scalar {
+  double _real;
+  double _imag;
+} ls_hs_scalar;
+void ls_hs_init(void);
+void ls_hs_exit(void);
+
+void ls_hs_set_exception_handler(void (*handler)(char const *message));
+void ls_hs_error(char const *message);
+typedef struct ls_hs_symmetry ls_hs_symmetry;
+
+typedef struct ls_hs_symmetries ls_hs_symmetries;
+
+typedef enum ls_hs_particle_type {
+  LS_HS_SPIN,
+  LS_HS_SPINFUL_FERMION,
+  LS_HS_SPINLESS_FERMION
+} ls_hs_particle_type;
+
+typedef void (*ls_hs_internal_state_index_kernel_type)(
+    ptrdiff_t batch_size, uint64_t const *alphas, ptrdiff_t alphas_stride,
+    ptrdiff_t *indices, ptrdiff_t indices_stride, void const *private_data);
+
+typedef void (*ls_hs_internal_is_representative_kernel_type)(
+    ptrdiff_t batch_size, uint64_t const *alphas, ptrdiff_t alphas_stride,
+    uint8_t *are_representatives, double *norms, void const *private_data);
+
+typedef void (*ls_hs_internal_state_info_kernel_type)(
+    ptrdiff_t batch_size, uint64_t const *alphas, ptrdiff_t alphas_stride,
+    uint64_t *betas, ptrdiff_t betas_stride, ls_hs_scalar *characters,
+    double *norms, void const *private_data);
+
+typedef struct ls_hs_basis_kernels {
+  ls_hs_internal_state_info_kernel_type state_info_kernel;
+  void *state_info_data;
+  ls_hs_internal_is_representative_kernel_type is_representative_kernel;
+  void *is_representative_data;
+  ls_hs_internal_state_index_kernel_type state_index_kernel;
+  void *state_index_data;
+} ls_hs_basis_kernels;
+
+typedef struct ls_hs_permutation_group {
+  _Atomic int refcount;
+  int number_bits;
+  int number_shifts;
+  int number_masks;
+  uint64_t *masks;
+  uint64_t *shifts;
+  double *eigvals_re;
+  double *eigvals_im;
+  void *haskell_payload;
+} ls_hs_permutation_group;
+
+typedef struct ls_hs_basis {
+  _Atomic int refcount;
+  int number_sites;
+  int number_particles;
+  int number_up;
+  ls_hs_particle_type particle_type;
+  int spin_inversion;
+  bool state_index_is_identity;
+  bool requires_projection;
+  ls_hs_basis_kernels *kernels;
+  chpl_external_array representatives;
+  void *haskell_payload;
+} ls_hs_basis;
+
+typedef struct ls_hs_expr ls_hs_expr;
+
+typedef void (*ls_hs_index_replacement_type)(int spin, int site, int *new_spin,
+                                             int *new_site);
+
+typedef struct ls_hs_nonbranching_terms {
+  int number_terms;
+  int number_bits;
+  // number_words = ceil(number_bits / 64)
+  ls_hs_scalar const *v; // array of shape [number_terms]
+  uint64_t const *m;     // array of shape [number_terms, number_words]
+  uint64_t const *l;     // array of shape [number_terms, number_words]
+  uint64_t const *r;     // array of shape [number_terms, number_words]
+  uint64_t const *x;     // array of shape [number_terms, number_words]
+  uint64_t const *s;     // array of shape [number_terms, number_words]
+  // all arrays are contiguous in row-major order
+} ls_hs_nonbranching_terms;
+
+typedef struct ls_hs_operator {
+  _Atomic int refcount;
+  ls_hs_basis const *basis;
+  ls_hs_nonbranching_terms const *off_diag_terms;
+  ls_hs_nonbranching_terms const *diag_terms;
+  // ls_internal_operator_kernel_data const *apply_off_diag_cxt;
+  // ls_internal_operator_kernel_data const *apply_diag_cxt;
+  void *haskell_payload;
+} ls_hs_operator;
+
+typedef struct ls_hs_yaml_config {
+  ls_hs_basis const *basis;
+  ls_hs_operator const *hamiltonian;
+  int number_observables;
+  ls_hs_operator const *const *observables;
+} ls_hs_yaml_config;
+
+typedef struct ls_chpl_kernels {
+  void (*enumerate_states)(ls_hs_basis const *, uint64_t, uint64_t,
+                           chpl_external_array *);
+  void (*operator_apply_off_diag)(ls_hs_operator *, int64_t, uint64_t *,
+                                  chpl_external_array *, chpl_external_array *,
+                                  chpl_external_array *, int64_t);
+  void (*operator_apply_diag)(ls_hs_operator *, int64_t, uint64_t *,
+                              chpl_external_array *, int64_t);
+  void (*matrix_vector_product)(ls_hs_operator *, int, double const *,
+                                double *);
+} ls_chpl_kernels;
+void ls_chpl_init(void);
+void ls_chpl_finalize(void);
+ls_chpl_kernels const *ls_hs_internal_get_chpl_kernels(void);
+void ls_hs_internal_set_chpl_kernels(ls_chpl_kernels const *kernels);
+void ls_hs_internal_destroy_external_array(chpl_external_array *arr);
+// The following should, in principle, not exist
+void ls_hs_state_index(ls_hs_basis const *const basis,
+                       ptrdiff_t const batch_size,
+                       uint64_t const *const restrict spins,
+                       ptrdiff_t const spins_stride,
+                       ptrdiff_t *const restrict indices,
+                       ptrdiff_t const indices_stride);
+void ls_hs_is_representative(ls_hs_basis const *basis, ptrdiff_t batch_size,
+                             uint64_t const *alphas, ptrdiff_t alphas_stride,
+                             uint8_t *are_representatives, double *norms);
+void ls_hs_state_info(ls_hs_basis const *basis, ptrdiff_t batch_size,
+                      uint64_t const *alphas, ptrdiff_t alphas_stride,
+                      uint64_t *betas, ptrdiff_t betas_stride,
+                      ls_hs_scalar *characters, double *norms);
+void ls_hs_build_representatives(ls_hs_basis *basis, uint64_t const lower,
+                                 uint64_t const upper);
+void ls_hs_unchecked_set_representatives(ls_hs_basis *basis,
+                                         chpl_external_array const *states,
+					 int cache_bits);
+ls_hs_symmetry * ls_hs_symmetry_from_json(char const *);
+void ls_hs_destroy_symmetry(ls_hs_symmetry *);
+int ls_hs_symmetry_sector(ls_hs_symmetry const*);
+double ls_hs_symmetry_phase(ls_hs_symmetry const*);
+int ls_hs_symmetry_length(ls_hs_symmetry const*);
+int * ls_hs_symmetry_permutation(ls_hs_symmetry const*);
+void ls_hs_destroy_permutation(int *);
+ls_hs_symmetries * ls_hs_symmetries_from_json(char const *);
+void ls_hs_destroy_symmetries(ls_hs_symmetries *);
+ls_hs_basis * ls_hs_clone_basis(ls_hs_basis const*);
+void ls_hs_destroy_basis(ls_hs_basis *);
+ls_hs_basis * ls_hs_basis_from_json(char const *);
+char const * ls_hs_basis_to_json(ls_hs_basis const*);
+void ls_hs_destroy_string(char const *);
+uint64_t ls_hs_min_state_estimate(ls_hs_basis const*);
+uint64_t ls_hs_max_state_estimate(ls_hs_basis const*);
+bool ls_hs_basis_has_fixed_hamming_weight(ls_hs_basis const*);
+bool ls_hs_basis_has_spin_inversion_symmetry(ls_hs_basis const*);
+bool ls_hs_basis_has_permutation_symmetries(ls_hs_basis const*);
+bool ls_hs_basis_requires_projection(ls_hs_basis const*);
+void ls_hs_basis_build(ls_hs_basis const*);
+bool ls_hs_basis_is_built(ls_hs_basis const*);
+int ls_hs_basis_number_bits(ls_hs_basis const*);
+int ls_hs_basis_number_words(ls_hs_basis const*);
+char const * ls_hs_basis_state_to_string(ls_hs_basis const*, uint64_t const*);
+ptrdiff_t ls_hs_fixed_hamming_state_to_index(uint64_t);
+uint64_t ls_hs_fixed_hamming_index_to_state(ptrdiff_t, int);
+char const * ls_hs_expr_to_json(ls_hs_expr const*);
+ls_hs_expr * ls_hs_expr_from_json(char const *);
+void ls_hs_destroy_expr(ls_hs_expr *);
+char const * ls_hs_expr_to_string(ls_hs_expr const*);
+ls_hs_expr * ls_hs_expr_plus(ls_hs_expr const*, ls_hs_expr const*);
+ls_hs_expr * ls_hs_expr_minus(ls_hs_expr const*, ls_hs_expr const*);
+ls_hs_expr * ls_hs_expr_times(ls_hs_expr const*, ls_hs_expr const*);
+ls_hs_expr * ls_hs_expr_scale(ls_hs_scalar const*, ls_hs_expr const*);
+ls_hs_expr * ls_hs_replace_indices(ls_hs_expr const*, ls_hs_index_replacement_type);
+bool ls_hs_expr_equal(ls_hs_expr const*, ls_hs_expr const*);
+ls_hs_expr * ls_hs_expr_adjoint(ls_hs_expr const*);
+bool ls_hs_expr_is_hermitian(ls_hs_expr const*);
+bool ls_hs_expr_is_real(ls_hs_expr const*);
+bool ls_hs_expr_is_identity(ls_hs_expr const*);
+ls_hs_operator * ls_hs_create_operator(ls_hs_basis const*, ls_hs_expr const*);
+ls_hs_operator * ls_hs_clone_operator(ls_hs_operator const*);
+void ls_hs_destroy_operator(ls_hs_operator *);
+int ls_hs_operator_max_number_off_diag(ls_hs_operator *);
+ls_hs_expr * ls_hs_operator_get_expr(ls_hs_operator const*);
+ls_hs_basis * ls_hs_operator_get_basis(ls_hs_operator const*);
+void ls_hs_prepare_hphi(ls_hs_operator const*, char const *);
+void ls_hs_prepare_mvmc(ls_hs_operator const*, char const *);
+ls_hs_yaml_config * ls_hs_load_yaml_config(char const *);
+void ls_hs_destroy_yaml_config(ls_hs_yaml_config *);
diff --git a/python/tutorial.md b/python/tutorial.md
new file mode 100644
index 0000000..a21650a
--- /dev/null
+++ b/python/tutorial.md
@@ -0,0 +1,29 @@
+# lattice_symmetries tutorial
+
+##Contents
+
+*[Installing]
+*[Introduction]
+*[Examples]
+  *[ED]
+
+## Installing
+
+The installation of lattice_symmetries requires Nix. 
+Therefore, at first you need to install Nix. Then you need to fork the last version from guthub.
+The last step is to configure Nix, so that every depencence works nicely. Voila, and you have a working package.
+
+## Introduction
+
+The lattice_symmetries is a powerful package that nicely works with many-body quantum spin systems
+and takes into account system symmetries.
+The lattice_symmetries implements fast matrix-vector multiplication that can be applied to various problems, 
+such as time evolution or exact diagonalization of many-body Hamiltonians. 
+
+## Examples
+
+Here we will take a look at different examples of lattice_symmetries applications.
+
+### ED
+
+Now we will consider the simplest example of exact diagonalization.

From 3ea2125a1a76f18e2adbbf507ad2fae00f979d93 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Wed, 21 Feb 2024 17:34:00 +0100
Subject: [PATCH 02/67] md test

---
 python/tutorial.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index a21650a..4d3cc58 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -2,10 +2,10 @@
 
 ##Contents
 
-*[Installing]
-*[Introduction]
-*[Examples]
-  *[ED]
+*[Installing] (#installing)
+*[Introduction] (#introduction)
+*[Examples] (#examples)
+  *[ED] (#ed)
 
 ## Installing
 

From 81b503c77003eee78ee4c307c6b88913f9cdaee0 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Wed, 21 Feb 2024 17:39:12 +0100
Subject: [PATCH 03/67] test md

---
 python/tutorial.md | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 4d3cc58..0e79e27 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -1,11 +1,11 @@
 # lattice_symmetries tutorial
 
-##Contents
+## Contents
 
-*[Installing] (#installing)
-*[Introduction] (#introduction)
-*[Examples] (#examples)
-  *[ED] (#ed)
+-[Installing] (#installing)
+-[Introduction] (#introduction)
+-[Examples] (#examples)
+-[ED] (#ed)
 
 ## Installing
 
@@ -15,9 +15,9 @@ The last step is to configure Nix, so that every depencence works nicely. Voila,
 
 ## Introduction
 
-The lattice_symmetries is a powerful package that nicely works with many-body quantum spin systems
+The `lattice_symmetries` is a powerful package that nicely works with many-body quantum spin systems
 and takes into account system symmetries.
-The lattice_symmetries implements fast matrix-vector multiplication that can be applied to various problems, 
+The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
 such as time evolution or exact diagonalization of many-body Hamiltonians. 
 
 ## Examples

From 20e7c7d98afbf13f9f6941a2033b2e3e535de977 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Wed, 21 Feb 2024 17:41:32 +0100
Subject: [PATCH 04/67] test md

---
 python/tutorial.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/python/tutorial.md b/python/tutorial.md
index 0e79e27..d5607df 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -7,12 +7,14 @@
 -[Examples] (#examples)
 -[ED] (#ed)
 
+<div id='installing'/>
 ## Installing
 
 The installation of lattice_symmetries requires Nix. 
 Therefore, at first you need to install Nix. Then you need to fork the last version from guthub.
 The last step is to configure Nix, so that every depencence works nicely. Voila, and you have a working package.
 
+<div id='introduction'/>
 ## Introduction
 
 The `lattice_symmetries` is a powerful package that nicely works with many-body quantum spin systems
@@ -20,10 +22,12 @@ and takes into account system symmetries.
 The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
 such as time evolution or exact diagonalization of many-body Hamiltonians. 
 
+<div id='examples'/>
 ## Examples
 
 Here we will take a look at different examples of lattice_symmetries applications.
 
+<div id='ed'/>
 ### ED
 
 Now we will consider the simplest example of exact diagonalization.

From 49c0a43af157fbebf5efa871bb43c0233d006817 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Wed, 21 Feb 2024 17:46:20 +0100
Subject: [PATCH 05/67] test md

---
 python/tutorial.md | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index d5607df..95b497d 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -2,19 +2,17 @@
 
 ## Contents
 
--[Installing] (#installing)
--[Introduction] (#introduction)
--[Examples] (#examples)
--[ED] (#ed)
+- [Installing](#Installing)
+- [Introduction](#Introduction)
+- [Examples](#Examples)
+    - [ED](#ED)
 
-<div id='installing'/>
 ## Installing
 
 The installation of lattice_symmetries requires Nix. 
 Therefore, at first you need to install Nix. Then you need to fork the last version from guthub.
 The last step is to configure Nix, so that every depencence works nicely. Voila, and you have a working package.
 
-<div id='introduction'/>
 ## Introduction
 
 The `lattice_symmetries` is a powerful package that nicely works with many-body quantum spin systems
@@ -22,12 +20,10 @@ and takes into account system symmetries.
 The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
 such as time evolution or exact diagonalization of many-body Hamiltonians. 
 
-<div id='examples'/>
 ## Examples
 
 Here we will take a look at different examples of lattice_symmetries applications.
 
-<div id='ed'/>
 ### ED
 
 Now we will consider the simplest example of exact diagonalization.

From 8dfaa3927f2714f7b2bb81fecc97f2949f3d22e4 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Thu, 22 Feb 2024 16:11:34 +0100
Subject: [PATCH 06/67] installing description

---
 python/tutorial.md | 49 ++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 41 insertions(+), 8 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 95b497d..65e1ce8 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -9,16 +9,49 @@
 
 ## Installing
 
-The installation of lattice_symmetries requires Nix. 
-Therefore, at first you need to install Nix. Then you need to fork the last version from guthub.
-The last step is to configure Nix, so that every depencence works nicely. Voila, and you have a working package.
+The first step before we can apply `lattice_symmetries` is to install Nix. The installation process slightly depends on your system and can be found in the [official documentaion](https://nix.dev/install-nix#install-nix).
+In the case of WSL (Windows Subsystem for Linux), it is preferred to install a single-user version:
 
-## Introduction
+```sh
+curl -L https://nixos.org/nix/install | sh -s -- --no-daemon
+```
+
+You can check that Nix is installed by opening a new terminal and typing:
+
+```sh
+nix --version
+```
+
+After that, it is necessary to add the support of flakes and experimental features into the configuration Nix file `nix.conf`. The file can be located in one of the two paths (starting from the root directory): `~/.config/nix/nix.conf` or `\etc\nix\nix.conf`. If the file is absent, you should create it in one of the mentioned directories. The file should contain the following line:
+
+```sh
+experimental-features = nix-command flakes
+```
+Now, the Nix is ready.
+
+The next step is to actually install `lattice_symmetries`. For that, you need to clone the `lattice_symmetries` from github:
+
+```sh
+git clone https://github.com/twesterhout/lattice-symmetries.git
+cd lattice-symmetries
+git submodule update --init --recursive
+```
+
+Then you should move to the `lattice-symmetries/python` folder and prepare Nix files:
+```sh
+cd python
+patchPhase || nix
+```
+
+Now you can build everything using Nix:
+
+```sh
+nix develop .#python
+```
+
+Voila, and you have the working package.
+If you open a new terminal, the last step should be repeated ,i.e., moving to  `lattice-symmetries/python` and typing `nix develop .#python`.
 
-The `lattice_symmetries` is a powerful package that nicely works with many-body quantum spin systems
-and takes into account system symmetries.
-The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
-such as time evolution or exact diagonalization of many-body Hamiltonians. 
 
 ## Examples
 

From 4713b1c02ac19237e0e0cd8f2e1e9fd210fb14a1 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Thu, 22 Feb 2024 16:16:55 +0100
Subject: [PATCH 07/67] intro returned

---
 python/tutorial.md | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 65e1ce8..10e00f3 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -1,4 +1,4 @@
-# lattice_symmetries tutorial
+# `lattice_symmetries` tutorial
 
 ## Contents
 
@@ -52,6 +52,12 @@ nix develop .#python
 Voila, and you have the working package.
 If you open a new terminal, the last step should be repeated ,i.e., moving to  `lattice-symmetries/python` and typing `nix develop .#python`.
 
+## Introduction
+
+The `lattice_symmetries` is a powerful package that nicely works with many-body quantum spin systems
+and takes into account system symmetries.
+The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
+such as time evolution or exact diagonalization of many-body Hamiltonians. 
 
 ## Examples
 

From 9babfa793fe4e40953c26de476d3fb46875b3552 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Tue, 27 Feb 2024 14:23:18 +0100
Subject: [PATCH 08/67] tutorial update

---
 python/tutorial.md | 40 +++++++++++++++++++++++++++++++++++-----
 1 file changed, 35 insertions(+), 5 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 10e00f3..02b0116 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -3,9 +3,15 @@
 ## Contents
 
 - [Installing](#Installing)
-- [Introduction](#Introduction)
+- [Introduction, basic concepts and functions](#Introduction, basic concepts and functions)
+    -[Basis](#Basis)
+    -[Expressions](#Expression)
+    -[Operators](#Operators)
+    -[Symmetries](#Symmetries)
 - [Examples](#Examples)
-    - [ED](#ED)
+    - [Simple ED](#Simple ED)
+    - [More complicated ED](#More complicated ED)
+    - [Time evolution](#Time evolution)
 
 ## Installing
 
@@ -52,17 +58,41 @@ nix develop .#python
 Voila, and you have the working package.
 If you open a new terminal, the last step should be repeated ,i.e., moving to  `lattice-symmetries/python` and typing `nix develop .#python`.
 
-## Introduction
+## Introduction, basic concepts and functions
 
 The `lattice_symmetries` is a powerful package that nicely works with many-body quantum spin systems
 and takes into account system symmetries.
 The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
-such as time evolution or exact diagonalization of many-body Hamiltonians. 
+such as time evolution or exact diagonalization of many-body Hamiltonians.
+
+The basic objects upon which the machinery is built are Basis, Expression, Operator, and Symmetry.
+- Basis objects allow to construct basis of many-body hamiltonians (spin and fermion) with appropriate symmetries.
+The Basis objects als
+- Expression objects is a nice way to work with symbolic representations of operators. It is possible to add
+
+### Basis
+
+### Expressions
+
+### Operators
+
+### Symmetries
+
+The full power of `lattice_symmetries` manifests if one use symmetries when constructing 
+symmetry adapted basis and linear operators acting on teh corresponding Hilbert space.
 
 ## Examples
 
 Here we will take a look at different examples of lattice_symmetries applications.
 
-### ED
+### Simple ED
 
 Now we will consider the simplest example of exact diagonalization.
+
+### More complicated ED
+
+Now let's consider a more complicated example of ED.
+
+### Time evolution
+
+Another example of capabilities of `lattice_symmetries` is time evolution.

From ece4ae53f4ecfe236c907ad86a1b0c05a343fbe7 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Tue, 27 Feb 2024 14:39:23 +0100
Subject: [PATCH 09/67] tutorial update

---
 python/tutorial.md | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 02b0116..7dbfd1c 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -3,15 +3,15 @@
 ## Contents
 
 - [Installing](#Installing)
-- [Introduction, basic concepts and functions](#Introduction, basic concepts and functions)
-    -[Basis](#Basis)
-    -[Expressions](#Expression)
-    -[Operators](#Operators)
-    -[Symmetries](#Symmetries)
+- [Introduction, basic concepts and functions](#Intro)
+    - [Basis](#Basis)
+    - [Expressions](#Expressions)
+    - [Operators](#Operators)
+    - [Symmetries](#Symmetries)
 - [Examples](#Examples)
-    - [Simple ED](#Simple ED)
-    - [More complicated ED](#More complicated ED)
-    - [Time evolution](#Time evolution)
+    - [Simple ED](#Simple_ED)
+    - [More complicated ED](#Complex_ED)
+    - [Time evolution](#Time_evolution)
 
 ## Installing
 

From 1563ee57f620735e1c8aa54b29c543996c6f1354 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Fri, 1 Mar 2024 11:57:55 +0100
Subject: [PATCH 10/67] some test tests

---
 python/test/casual_tests.py | 40 +++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)
 create mode 100644 python/test/casual_tests.py

diff --git a/python/test/casual_tests.py b/python/test/casual_tests.py
new file mode 100644
index 0000000..2591518
--- /dev/null
+++ b/python/test/casual_tests.py
@@ -0,0 +1,40 @@
+import glob
+import os
+import lattice_symmetries as ls
+import numpy as np
+import scipy.sparse.linalg
+#from pytest import approx
+
+
+def test_index():
+    basis = ls.SpinBasis(4)
+    basis.build()
+    print(basis.index(basis.states))
+    print(basis.number_states)
+    for i in range(basis.number_states):
+        print(i, basis.states[i], basis.state_to_string(i))
+    #assert np.array_equal(basis.index(basis.states), basis.states)
+    #assert np.array_equal(basis.index(basis.states[2]), 2)
+
+def test_operator_apply():
+    basis = ls.SpinBasis(2)
+    basis.build()
+    expr = ls.Expr("1.0 σᶻ₀ σᶻ₁ + 2.0 σ⁺₀ σ⁻₁ + 2.0 σ⁻₀ σ⁺₁")
+    hamiltonian = ls.Operator(basis, expr)
+    print(hamiltonian.expression)
+    print(hamiltonian.basis.state_to_string(0))
+    #len(hamiltonian.apply_off_diag_to_basis_state(basis.states[1]))
+
+
+def test_prepare_hphi():
+    basis = ls.SpinBasis(2)
+    expr = ls.Expr("σᶻ₀ σᶻ₁ + 2 σ⁺₀ σ⁻₁ + 2 σ⁻₀ σ⁺₁")
+    hamiltonian = ls.Operator(basis, expr)
+    hamiltonian.prepare_inputs_for_hphi("/tmp/lattice-symmetries-python/hphi")
+
+
+def main():
+    #test_index()
+    test_operator_apply()
+
+main()

From 6cb634f008d9263680a28abaf524f05cf90ba0df Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 15:11:26 +0100
Subject: [PATCH 11/67] tutorial update

---
 python/tutorial.md | 57 ++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 52 insertions(+), 5 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 7dbfd1c..b7cbf2e 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -5,6 +5,8 @@
 - [Installing](#Installing)
 - [Introduction, basic concepts and functions](#Intro)
     - [Basis](#Basis)
+        - [Spin basis](#Spin_basis)
+        - [Fermionic basis](Fermionic basis)
     - [Expressions](#Expressions)
     - [Operators](#Operators)
     - [Symmetries](#Symmetries)
@@ -65,17 +67,62 @@ and takes into account system symmetries.
 The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
 such as time evolution or exact diagonalization of many-body Hamiltonians.
 
-The basic objects upon which the machinery is built are Basis, Expression, Operator, and Symmetry.
-- Basis objects allow to construct basis of many-body hamiltonians (spin and fermion) with appropriate symmetries.
-The Basis objects als
-- Expression objects is a nice way to work with symbolic representations of operators. It is possible to add
+The basic objects upon which the machinery is built are `Basis`, `Expression`, `Operator`, and `Symmetry`.
+- `Basis` objects allow to construct basis of many-body Hilbert space consisting of spins or fermions with appropriate symmetries.
+Each basis state is represented as a sequence of $0$s and $1$s (i.e. bits), which can be interpreted as a binary number.
+- `Expression` objects is a nice way to work with symbolic representations of operators. It is possible to sum expressions, multiply them to numbers and each other.
+`Expressions` allow not to think about an explicit matrix representation of an operator, and user can work directly with analytical formulae as if they are written in paper. 
+- `Operator` object is an actual operator that can act on individual basis states and their linear combinations. 
+- `Symmetry` is a method to work with symmetry adapted basises.
+If an operator has symmmetries, it is useful to work in symmetry-adapted basis, since it can drastically decrease the dimension of the Hilbert space.
+
+Now we will take a look at basic functions and methods for these objects. 
 
 ### Basis
+#### Spin basis
+Let's take a look at simple examples, and at first we will not consider additional symmetries.
+The simplest example would be a spin basis:
+
+`
+import lattice_symmetries as ls
+ 
+basis = ls.SpinBasis(3) #We create basis consisting of three $\frac{1}{2}$-spins (each spin can be |0⟩ or |1⟩) 
+basis.build() #build basis
+print(basis.index(basis.states)) #Print array of indices of basis states
+print(basis.number_states) #Print total number of states in the basis
+for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
+    print(basis.states[i], basis.state_to_string(i)) #The last function represents basis states as a sequence of $0$s and $1$s  
+`
+The result looks like:
+`
+[0 1 2 3 4 5 6 7]
+8
+0 |000⟩
+1 |001⟩
+2 |010⟩
+3 |011⟩
+4 |100⟩
+5 |101⟩
+6 |110⟩
+7 |111⟩
+`
+The basis states are equal to their indices and binary represenations as they should.
+
+#### Fermionic basis
+We can also consider fermionic basis:
+
 
 ### Expressions
 
+Let's take a look at a few examples.
+
+
 ### Operators
 
+Based on an expression and a basis, we can build the corresponding operator acting on the chosen basis. Let's consider an example.
+
+
+
 ### Symmetries
 
 The full power of `lattice_symmetries` manifests if one use symmetries when constructing 
@@ -83,7 +130,7 @@ symmetry adapted basis and linear operators acting on teh corresponding Hilbert
 
 ## Examples
 
-Here we will take a look at different examples of lattice_symmetries applications.
+Here we will take a look at different examples of `lattice_symmetries` applications.
 
 ### Simple ED
 

From 54100783d14708f10c1fdc44c12f9c09752fc73c Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 15:14:52 +0100
Subject: [PATCH 12/67] tutorial update

---
 python/tutorial.md | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index b7cbf2e..c21e064 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -6,7 +6,7 @@
 - [Introduction, basic concepts and functions](#Intro)
     - [Basis](#Basis)
         - [Spin basis](#Spin_basis)
-        - [Fermionic basis](Fermionic basis)
+        - [Fermionic basis](#Fermionic basis)
     - [Expressions](#Expressions)
     - [Operators](#Operators)
     - [Symmetries](#Symmetries)
@@ -83,7 +83,7 @@ Now we will take a look at basic functions and methods for these objects.
 Let's take a look at simple examples, and at first we will not consider additional symmetries.
 The simplest example would be a spin basis:
 
-`
+```sh
 import lattice_symmetries as ls
  
 basis = ls.SpinBasis(3) #We create basis consisting of three $\frac{1}{2}$-spins (each spin can be |0⟩ or |1⟩) 
@@ -92,9 +92,10 @@ print(basis.index(basis.states)) #Print array of indices of basis states
 print(basis.number_states) #Print total number of states in the basis
 for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
     print(basis.states[i], basis.state_to_string(i)) #The last function represents basis states as a sequence of $0$s and $1$s  
-`
+```
+
 The result looks like:
-`
+```sh
 [0 1 2 3 4 5 6 7]
 8
 0 |000⟩
@@ -105,7 +106,7 @@ The result looks like:
 5 |101⟩
 6 |110⟩
 7 |111⟩
-`
+```
 The basis states are equal to their indices and binary represenations as they should.
 
 #### Fermionic basis

From 1ed938ce654a482589b4704d53d63debf188ce15 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 15:30:56 +0100
Subject: [PATCH 13/67] tutorial update

---
 python/tutorial.md | 25 +++++++++++++------------
 1 file changed, 13 insertions(+), 12 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index c21e064..4f3db98 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -5,15 +5,15 @@
 - [Installing](#Installing)
 - [Introduction, basic concepts and functions](#Intro)
     - [Basis](#Basis)
-        - [Spin basis](#Spin_basis)
-        - [Fermionic basis](#Fermionic basis)
+        - [Spin basis](#Spin-basis)
+        - [Fermionic basis](#Fermionic-basis)
     - [Expressions](#Expressions)
     - [Operators](#Operators)
     - [Symmetries](#Symmetries)
 - [Examples](#Examples)
-    - [Simple ED](#Simple_ED)
-    - [More complicated ED](#Complex_ED)
-    - [Time evolution](#Time_evolution)
+    - [Simple ED](#Simple-ED)
+    - [More complicated ED](#More-complicated-ED)
+    - [Time evolution](#Time-evolution)
 
 ## Installing
 
@@ -68,10 +68,10 @@ The `lattice_symmetries` implements fast matrix-vector multiplication that can b
 such as time evolution or exact diagonalization of many-body Hamiltonians.
 
 The basic objects upon which the machinery is built are `Basis`, `Expression`, `Operator`, and `Symmetry`.
-- `Basis` objects allow to construct basis of many-body Hilbert space consisting of spins or fermions with appropriate symmetries.
-Each basis state is represented as a sequence of $0$s and $1$s (i.e. bits), which can be interpreted as a binary number.
-- `Expression` objects is a nice way to work with symbolic representations of operators. It is possible to sum expressions, multiply them to numbers and each other.
-`Expressions` allow not to think about an explicit matrix representation of an operator, and user can work directly with analytical formulae as if they are written in paper. 
+- `Basis` object stores a basis of a many-body Hilbert space consisting of spins or fermions with appropriate symmetries.
+Each basis state is represented as a sequence of 0s and 1s (i.e. a sequence of bits), which can be interpreted as a binary number.
+- `Expression` object is a nice way to work with symbolic representations of operators. It is possible to sum expressions, and multiply them by numbers and each other.
+`Expression` allows not to think about an explicit matrix representation of an operator, and the user can work directly with analytical formulae. 
 - `Operator` object is an actual operator that can act on individual basis states and their linear combinations. 
 - `Symmetry` is a method to work with symmetry adapted basises.
 If an operator has symmmetries, it is useful to work in symmetry-adapted basis, since it can drastically decrease the dimension of the Hilbert space.
@@ -80,7 +80,8 @@ Now we will take a look at basic functions and methods for these objects.
 
 ### Basis
 #### Spin basis
-Let's take a look at simple examples, and at first we will not consider additional symmetries.
+Let's look at simple examples; at first we will not consider additional symmetries.
+
 The simplest example would be a spin basis:
 
 ```sh
@@ -91,7 +92,7 @@ basis.build() #build basis
 print(basis.index(basis.states)) #Print array of indices of basis states
 print(basis.number_states) #Print total number of states in the basis
 for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
-    print(basis.states[i], basis.state_to_string(i)) #The last function represents basis states as a sequence of $0$s and $1$s  
+    print(basis.states[i], basis.state_to_string(i)) #The last function represents basis states as a sequence of 0s and 1s  
 ```
 
 The result looks like:
@@ -107,7 +108,7 @@ The result looks like:
 6 |110⟩
 7 |111⟩
 ```
-The basis states are equal to their indices and binary represenations as they should.
+The basis states are equal to their indices and binary representations as they should be.
 
 #### Fermionic basis
 We can also consider fermionic basis:

From 45486f9d4877d1ebce5a32628aeff7b8b2e1a18e Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 15:32:33 +0100
Subject: [PATCH 14/67] tutorial update

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 4f3db98..23b5cb6 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -3,7 +3,7 @@
 ## Contents
 
 - [Installing](#Installing)
-- [Introduction, basic concepts and functions](#Intro)
+- [Introduction, basic concepts and functions](#Introduction,-basic-concepts-and-functions)
     - [Basis](#Basis)
         - [Spin basis](#Spin-basis)
         - [Fermionic basis](#Fermionic-basis)

From 5a8905de5da22d39bc51e6175c7ef39448537914 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 15:33:53 +0100
Subject: [PATCH 15/67] tutorial update

---
 python/tutorial.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 23b5cb6..2186bea 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -3,7 +3,7 @@
 ## Contents
 
 - [Installing](#Installing)
-- [Introduction, basic concepts and functions](#Introduction,-basic-concepts-and-functions)
+- [Basic concepts and functions](#Basic-concepts-and-functions)
     - [Basis](#Basis)
         - [Spin basis](#Spin-basis)
         - [Fermionic basis](#Fermionic-basis)
@@ -60,7 +60,7 @@ nix develop .#python
 Voila, and you have the working package.
 If you open a new terminal, the last step should be repeated ,i.e., moving to  `lattice-symmetries/python` and typing `nix develop .#python`.
 
-## Introduction, basic concepts and functions
+## Basic concepts and functions
 
 The `lattice_symmetries` is a powerful package that nicely works with many-body quantum spin systems
 and takes into account system symmetries.

From 0f4bd02438f8eae7680adc2b647ea340936849d3 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 17:39:25 +0100
Subject: [PATCH 16/67] tutorial update

---
 python/tutorial.md | 115 +++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 111 insertions(+), 4 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 2186bea..20be5b3 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -6,7 +6,8 @@
 - [Basic concepts and functions](#Basic-concepts-and-functions)
     - [Basis](#Basis)
         - [Spin basis](#Spin-basis)
-        - [Fermionic basis](#Fermionic-basis)
+        - [Spinful Fermionic basis](#Spinful-fermionic-basis)
+        - [Spinless Fermionic basis](#Spinless-fermionic-basis)
     - [Expressions](#Expressions)
     - [Operators](#Operators)
     - [Symmetries](#Symmetries)
@@ -92,7 +93,7 @@ basis.build() #build basis
 print(basis.index(basis.states)) #Print array of indices of basis states
 print(basis.number_states) #Print total number of states in the basis
 for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
-    print(basis.states[i], basis.state_to_string(i)) #The last function represents basis states as a sequence of 0s and 1s  
+for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
 ```
 
 The result looks like:
@@ -110,8 +111,75 @@ The result looks like:
 ```
 The basis states are equal to their indices and binary representations as they should be.
 
-#### Fermionic basis
-We can also consider fermionic basis:
+We can also consider only a part of the whole Hilbert space with given number of spin ups (i.e. hamming weight of binary representaions).
+We can specify it as follows:
+```sh
+basis = ls.SpinBasis(3, hamming_weight=2) #We want the subspace with only 2 spins up 
+basis.build() #build basis
+for i in range(basis.number_states): #Print the states in the basis
+    print(basis.states[i], basis.state_to_string(basis.states[i]))  
+```
+The output is:
+```sh
+3 |011⟩
+5 |101⟩
+6 |110⟩
+```
+We see that basis states include only spins with 2 spins up. It is also interesting to note that in this case a basis state is not equal to its index.
+
+Sometimes our system has spin inversion symmetry, and we can additionally shorten our basis. In this case, we can specify it as follows:
+```sh
+basis = ls.SpinBasis(4, spin_inversion=-1) #Spin inversion is present. It is not nessecary to specify hamming weight here, since it is fixed by spin inversion symmetry. 
+basis.build() #build basis
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i]))  
+```
+We have the basis states:
+```sh
+3 |011⟩
+5 |101⟩
+6 |110⟩
+```
+
+#### Spinless Fermionic basis
+We can also consider basis of fermions without spins. The basis states are stored as integers as for spin basis, however the binary representation has a second quantization interpretation.
+Each basis state is given by the sequence of 0s and 1s, where 1 means that there is a fermion on the corresponding site, and 0 means that the site is vacant.
+Let's consider the simplest example of fermions on two sites:
+```sh
+basis = ls.SpinlessFermionBasis(2) #We create fermionic basis on 2 sites
+basis.build() #build basis
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+```
+which gives:
+```sh
+0 |00⟩
+1 |01⟩
+2 |10⟩
+3 |11⟩
+```
+as one would expect.
+
+We can specify the number of particles as well:
+```sh
+basis = ls.SpinlessFermionBasis(4, number_particles=3) #We create fermionic basis on 4 sites with only 3 fermions
+basis.build() #build basis
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+```
+which gives:
+```sh
+7 |0111⟩
+11 |1011⟩
+13 |1101⟩
+14 |1110⟩
+```
+We can see that the basis consists of states with three fermions.
+
+#### Spinful Fermionic basis
+
+The last case includes fermions with spin. The binary representations of basis states can be read as a pair of (fermions with spin up on lattice, fermions with spin down on a lattice).
+We can create a basis of spinful fermions as follows:
 
 
 ### Expressions
@@ -122,8 +190,47 @@ Let's take a look at a few examples.
 ### Operators
 
 Based on an expression and a basis, we can build the corresponding operator acting on the chosen basis. Let's consider an example.
+```sh
+basis = ls.SpinfulFermionBasis(2) #We create basis of spinful fermions on 2 sites
+basis.build() #build basis
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+```
+which gives:
+```sh
+0 |00⟩|00⟩
+1 |00⟩|01⟩
+2 |00⟩|10⟩
+3 |00⟩|11⟩
+4 |01⟩|00⟩
+5 |01⟩|01⟩
+6 |01⟩|10⟩
+7 |01⟩|11⟩
+8 |10⟩|00⟩
+9 |10⟩|01⟩
+10 |10⟩|10⟩
+11 |10⟩|11⟩
+12 |11⟩|00⟩
+13 |11⟩|01⟩
+14 |11⟩|10⟩
+15 |11⟩|11⟩
+```
+We see that now binary representation has the meaning of second quantization of fermions with two spins.
 
+As before, we can specify a sector of the total Hilbert space with given number of fermions with spin down and spin up:
 
+```sh
+basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins down and up (N_down, N_up)=(2,1)
+basis.build() #build basis
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+```
+which gives:
+```sh
+7 |01⟩|11⟩
+11 |10⟩|11⟩
+```
+as expected.
 
 ### Symmetries
 

From 57d887706d3279f29f64c7417d679e12a787e750 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 17:47:47 +0100
Subject: [PATCH 17/67] tutorial update

---
 python/tutorial.md | 39 ++++++++++++++++++++-------------------
 1 file changed, 20 insertions(+), 19 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 20be5b3..e6599f8 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -109,9 +109,9 @@ The result looks like:
 6 |110⟩
 7 |111⟩
 ```
-The basis states are equal to their indices and binary representations as they should be.
+The basis states are equal to their indices and binary representations, as they should be.
 
-We can also consider only a part of the whole Hilbert space with given number of spin ups (i.e. hamming weight of binary representaions).
+We can also consider only a part of the whole Hilbert space with a given number of spin ups (i.e. hamming weight of binary representations).
 We can specify it as follows:
 ```sh
 basis = ls.SpinBasis(3, hamming_weight=2) #We want the subspace with only 2 spins up 
@@ -127,7 +127,7 @@ The output is:
 ```
 We see that basis states include only spins with 2 spins up. It is also interesting to note that in this case a basis state is not equal to its index.
 
-Sometimes our system has spin inversion symmetry, and we can additionally shorten our basis. In this case, we can specify it as follows:
+Sometimes, our system has spin inversion symmetry, and we can additionally shorten our basis. In this case, we can specify it as follows:
 ```sh
 basis = ls.SpinBasis(4, spin_inversion=-1) #Spin inversion is present. It is not nessecary to specify hamming weight here, since it is fixed by spin inversion symmetry. 
 basis.build() #build basis
@@ -142,8 +142,8 @@ We have the basis states:
 ```
 
 #### Spinless Fermionic basis
-We can also consider basis of fermions without spins. The basis states are stored as integers as for spin basis, however the binary representation has a second quantization interpretation.
-Each basis state is given by the sequence of 0s and 1s, where 1 means that there is a fermion on the corresponding site, and 0 means that the site is vacant.
+We can also consider the basis of fermions without spins. The basis states are stored as integers as for spin basis. However the binary representation has a second quantization interpretation.
+Each basis state is given by the sequence of 0s and 1s, where 1 means a fermion on the corresponding site, and 0 means that the site is vacant.
 Let's consider the simplest example of fermions on two sites:
 ```sh
 basis = ls.SpinlessFermionBasis(2) #We create fermionic basis on 2 sites
@@ -156,7 +156,7 @@ which gives:
 0 |00⟩
 1 |01⟩
 2 |10⟩
-3 |11⟩
+3 |11⟩111
 ```
 as one would expect.
 
@@ -178,18 +178,8 @@ We can see that the basis consists of states with three fermions.
 
 #### Spinful Fermionic basis
 
-The last case includes fermions with spin. The binary representations of basis states can be read as a pair of (fermions with spin up on lattice, fermions with spin down on a lattice).
+The last case includes fermions with spin. The binary representations of basis states can be read as a pair of (fermions with spin up on a lattice, fermions with spin down on a lattice).
 We can create a basis of spinful fermions as follows:
-
-
-### Expressions
-
-Let's take a look at a few examples.
-
-
-### Operators
-
-Based on an expression and a basis, we can build the corresponding operator acting on the chosen basis. Let's consider an example.
 ```sh
 basis = ls.SpinfulFermionBasis(2) #We create basis of spinful fermions on 2 sites
 basis.build() #build basis
@@ -215,9 +205,9 @@ which gives:
 14 |11⟩|10⟩
 15 |11⟩|11⟩
 ```
-We see that now binary representation has the meaning of second quantization of fermions with two spins.
+We see that binary representation now means the second quantization of fermions with two spins.
 
-As before, we can specify a sector of the total Hilbert space with given number of fermions with spin down and spin up:
+As before, we can specify a sector of the total Hilbert space with a given number of fermions with spin down and spin up:
 
 ```sh
 basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins down and up (N_down, N_up)=(2,1)
@@ -232,6 +222,17 @@ which gives:
 ```
 as expected.
 
+### Expressions
+
+Let's take a look at a few examples.
+
+
+### Operators
+
+Based on an expression and a basis, we can build the corresponding operator acting on the chosen basis. Let's consider an example.
+
+
+
 ### Symmetries
 
 The full power of `lattice_symmetries` manifests if one use symmetries when constructing 

From 39dc1236178e1c7e40fc2a3d5ce780e9f5f4eb96 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 17:50:41 +0100
Subject: [PATCH 18/67] tutorial update

---
 python/tutorial.md | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index e6599f8..9d3b2ce 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -93,7 +93,7 @@ basis.build() #build basis
 print(basis.index(basis.states)) #Print array of indices of basis states
 print(basis.number_states) #Print total number of states in the basis
 for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
-for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
+    print(basis.states[i], basis.state_to_string(basis.states[i]))  
 ```
 
 The result looks like:
@@ -115,7 +115,7 @@ We can also consider only a part of the whole Hilbert space with a given number
 We can specify it as follows:
 ```sh
 basis = ls.SpinBasis(3, hamming_weight=2) #We want the subspace with only 2 spins up 
-basis.build() #build basis
+basis.build()
 for i in range(basis.number_states): #Print the states in the basis
     print(basis.states[i], basis.state_to_string(basis.states[i]))  
 ```
@@ -130,7 +130,7 @@ We see that basis states include only spins with 2 spins up. It is also interest
 Sometimes, our system has spin inversion symmetry, and we can additionally shorten our basis. In this case, we can specify it as follows:
 ```sh
 basis = ls.SpinBasis(4, spin_inversion=-1) #Spin inversion is present. It is not nessecary to specify hamming weight here, since it is fixed by spin inversion symmetry. 
-basis.build() #build basis
+basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i]))  
 ```
@@ -147,7 +147,7 @@ Each basis state is given by the sequence of 0s and 1s, where 1 means a fermion
 Let's consider the simplest example of fermions on two sites:
 ```sh
 basis = ls.SpinlessFermionBasis(2) #We create fermionic basis on 2 sites
-basis.build() #build basis
+basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 
 ```
@@ -163,7 +163,7 @@ as one would expect.
 We can specify the number of particles as well:
 ```sh
 basis = ls.SpinlessFermionBasis(4, number_particles=3) #We create fermionic basis on 4 sites with only 3 fermions
-basis.build() #build basis
+basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 
 ```
@@ -182,7 +182,7 @@ The last case includes fermions with spin. The binary representations of basis s
 We can create a basis of spinful fermions as follows:
 ```sh
 basis = ls.SpinfulFermionBasis(2) #We create basis of spinful fermions on 2 sites
-basis.build() #build basis
+basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 
 ```
@@ -211,7 +211,7 @@ As before, we can specify a sector of the total Hilbert space with a given numbe
 
 ```sh
 basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins down and up (N_down, N_up)=(2,1)
-basis.build() #build basis
+basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 
 ```

From dd21e01d3fe591fd400bc23d799b5bf1dbfc0724 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 17:52:21 +0100
Subject: [PATCH 19/67] tutorial update

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 9d3b2ce..6dba906 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -6,8 +6,8 @@
 - [Basic concepts and functions](#Basic-concepts-and-functions)
     - [Basis](#Basis)
         - [Spin basis](#Spin-basis)
-        - [Spinful Fermionic basis](#Spinful-fermionic-basis)
         - [Spinless Fermionic basis](#Spinless-fermionic-basis)
+        - [Spinful Fermionic basis](#Spinful-fermionic-basis)
     - [Expressions](#Expressions)
     - [Operators](#Operators)
     - [Symmetries](#Symmetries)

From 9f2cc769d9c66aa8f0c3b1b6163e1d9e27c0a2fc Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Fri, 8 Mar 2024 17:54:24 +0100
Subject: [PATCH 20/67] tutorial update

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 6dba906..a4e8271 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -210,7 +210,7 @@ We see that binary representation now means the second quantization of fermions
 As before, we can specify a sector of the total Hilbert space with a given number of fermions with spin down and spin up:
 
 ```sh
-basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins down and up (N_down, N_up)=(2,1)
+basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins up and down (N_up, N_down)=(2,1)
 basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 

From b6014e5d3f2d1b7309e7e2db2de2a5f60e617b83 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Thu, 14 Mar 2024 14:49:27 +0100
Subject: [PATCH 21/67] a little bit added to Tom's readme

---
 python/tutorial.md | 92 +++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 91 insertions(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index a4e8271..615889d 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -224,7 +224,97 @@ as expected.
 
 ### Expressions
 
-Let's take a look at a few examples.
+Expressions is an easy way to work with operators (for example, Hamiltonians) on symbolic level using second quantization formalism.
+This means that you can use primitive operators such as $\sigma^x$, $\sigma^y$, and $\sigma^z$ to build expressions for your Hamiltonian and observables.
+It is also possible to sum different expressions and multiply them to each other to compose more complicated expressions.
+Let's consider at first the simplest examples. 
+
+#### Primitive operators
+
+At first we need to import `Expr` from lattice-symmetries:
+```sh
+from lattice_symmetries import Expr
+```
+Now we will consider primitives operators defined on site with number 0 as an example, but you can also construct primitive operators residing on other lattice sites:
+
+ - $\sigma^x$ and $S^x$:
+   ```pycon
+   >>> Expr("σˣ₀") == Expr("\\sigma^x_0") #It is possible to use different notations for Expr for primitive operators.
+   #Here we check that they agree. Index 0 means the index of the corresponding site.
+   True
+   >>> Expr("Sˣ₀") == 0.5 * Expr("σˣ₀") #We check that Sˣ₀ is the same as 0.5*σˣ₀.
+   True
+   >>> Expr("σˣ₀").to_dense() #It is also possible to visualize the expressions as matrices.
+   [[0, 1],
+    [1, 0]]
+
+   ```
+   Now, we will take a look at other Pauli and spin matrices.
+
+ - $\sigma^y$ and $S^y$:
+   ```pycon
+   >>> Expr("σʸ₀") == Expr("\\sigma^y_0")
+   True
+   >>> Expr("Sʸ₀") == 0.5 * ls.Expr("σʸ₀")
+   True
+   >>> Expr("σʸ₀").to_dense()
+   [[0, -1j],
+    [1j, 0]]
+   ```
+
+ - $\sigma^z$ and $S^z$:
+   ```pycon
+   >>> Expr("σᶻ₀") == Expr("\\sigma^z_0")
+   True
+   >>> Expr("Sᶻ₀") == 0.5 * Expr("σᶻ₀")
+   True
+   >>> Expr("σᶻ₀").to_dense()
+   [[1, 0],
+    [0, -1]]
+   ```
+    We see that everything works as one would expect.
+
+ - $\mathbb{1}$:
+   ```pycon
+   >>> Expr("I", particle="spin-1/2").to_dense()
+   [[1, 0],
+    [0, 1]]
+   ```
+   (*Note:* in this case, we need to explicitly specify the particle type because it cannot be deduced from the expression)
+
+#### Operator algebra
+
+Primitives can be combined using the `+`, `-`, and `*` operations to build more complex operators.
+Furthermore, expressions can also be multiplied by scalars from the left using the `*` operator.
+
+For example, here are a few ways to write down the Heisenberg interaction between sites 0 and 1
+
+$$
+\mathbf{S}_i \cdot \mathbf{S}_j = S^x_i S^x_j + S^y_i S^y_j + S^z_i S^z_j
+$$
+
+```pycon
+>>> Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁") == Expr("Sx0 Sx1 + Sy0 Sy1 + Sz0 Sz1")
+True
+>>> Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁") == \
+             Expr("Sˣ₀") * Expr("Sˣ₁") + Expr("Sʸ₀") * Expr("Sʸ₁") + Expr("Sᶻ₀") * Expr("Sᶻ₁")
+True
+>>> Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁") == Expr("0.25 (σˣ₀ σˣ₁ + σʸ₀ σʸ₁ + σᶻ₀ σᶻ₁)")
+True
+>>> Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁") == Expr("0.5 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + 0.25 σᶻ₀ σᶻ₁")
+True
+```
+
+Under the hood, lattice-symmetries rewrites all the expressions into the canonical form, simplifying stuff along the way:
+
+```pycon
+>>> str(Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁"))
+"0.25 σᶻ₀ σᶻ₁ + 0.5 σ⁺₀ σ⁻₁ + 0.5 σ⁻₀ σ⁺₁"
+>>> str(Expr("0.5 (σˣ₁ + 1im σʸ₁) - σ⁺₁"))
+"0.0 I"
+>>> str(Expr("σ⁺₁ σ⁺₁"))
+"0.0 I"
+```
 
 
 ### Operators

From 51111aa3959f89b3da8994d47bdd6ec3fdd5eb88 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Thu, 14 Mar 2024 15:06:13 +0100
Subject: [PATCH 22/67] few things

---
 python/tutorial.md | 30 +++++++++++++++++++++++++++---
 1 file changed, 27 insertions(+), 3 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 615889d..7f96ef2 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -9,6 +9,8 @@
         - [Spinless Fermionic basis](#Spinless-fermionic-basis)
         - [Spinful Fermionic basis](#Spinful-fermionic-basis)
     - [Expressions](#Expressions)
+        - [Primitive operators](#Primitive-operators)
+        - [Operator algebra](#Operator-algebra)
     - [Operators](#Operators)
     - [Symmetries](#Symmetries)
 - [Examples](#Examples)
@@ -224,7 +226,7 @@ as expected.
 
 ### Expressions
 
-Expressions is an easy way to work with operators (for example, Hamiltonians) on symbolic level using second quantization formalism.
+Expressions are an easy way to work with operators (for example, Hamiltonians) on a symbolic level using second quantization formalism.
 This means that you can use primitive operators such as $\sigma^x$, $\sigma^y$, and $\sigma^z$ to build expressions for your Hamiltonian and observables.
 It is also possible to sum different expressions and multiply them to each other to compose more complicated expressions.
 Let's consider at first the simplest examples. 
@@ -274,7 +276,7 @@ Now we will consider primitives operators defined on site with number 0 as an ex
    ```
     We see that everything works as one would expect.
 
- - $\mathbb{1}$:
+ - Identity operator $\mathbb{1}$:
    ```pycon
    >>> Expr("I", particle="spin-1/2").to_dense()
    [[1, 0],
@@ -282,6 +284,20 @@ Now we will consider primitives operators defined on site with number 0 as an ex
    ```
    (*Note:* in this case, we need to explicitly specify the particle type because it cannot be deduced from the expression)
 
+It is also possible to check various properties of expressions:
+   ```pycon
+   >>> a=Expr("\\sigma^z_0")
+   >>> a.is_real #If the corresponding matrix real
+   True
+   >>> a.is_hermitian #If the corresponding matrix hermitian
+   True
+   >>> a.is_identity #If the corresponding matrix identity
+   False
+   >>> a.number_sites #Number of sites in the expression
+   1
+   ```
+
+
 #### Operator algebra
 
 Primitives can be combined using the `+`, `-`, and `*` operations to build more complex operators.
@@ -290,7 +306,7 @@ Furthermore, expressions can also be multiplied by scalars from the left using t
 For example, here are a few ways to write down the Heisenberg interaction between sites 0 and 1
 
 $$
-\mathbf{S}_i \cdot \mathbf{S}_j = S^x_i S^x_j + S^y_i S^y_j + S^z_i S^z_j
+\mathbf{S}_i \cdot \mathbf{S}_j = S^x_0 S^x_1 + S^y_0 S^y_1 + S^z_0 S^z_1
 $$
 
 ```pycon
@@ -303,6 +319,14 @@ True
 True
 >>> Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁") == Expr("0.5 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + 0.25 σᶻ₀ σᶻ₁")
 True
+
+#Above we checked that different definitions of this interaction agree. Let's take a look at the corresponding matrix:
+>>> Expr("Sx0 Sx1 + Sy0 Sy1 + Sz0 Sz1").to_dense()
+
+[[ 0.25  0.    0.    0.  ]
+ [ 0.   -0.25  0.5   0.  ]
+ [ 0.    0.5  -0.25  0.  ]
+ [ 0.    0.    0.    0.25]]
 ```
 
 Under the hood, lattice-symmetries rewrites all the expressions into the canonical form, simplifying stuff along the way:

From 4f7c80e256af4eafcb2d9ef42f8a3d00dd3065c7 Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Thu, 14 Mar 2024 15:19:53 +0100
Subject: [PATCH 23/67] a little more added

---
 python/tutorial.md | 33 ++++++++++++++++++++-------------
 1 file changed, 20 insertions(+), 13 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 7f96ef2..dab2c6b 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -11,6 +11,8 @@
     - [Expressions](#Expressions)
         - [Primitive operators](#Primitive-operators)
         - [Operator algebra](#Operator-algebra)
+        - [Complex expressions](#Complex-expressions)
+        - [Properties](#Properties)
     - [Operators](#Operators)
     - [Symmetries](#Symmetries)
 - [Examples](#Examples)
@@ -158,7 +160,7 @@ which gives:
 0 |00⟩
 1 |01⟩
 2 |10⟩
-3 |11⟩111
+3 |11⟩
 ```
 as one would expect.
 
@@ -284,18 +286,6 @@ Now we will consider primitives operators defined on site with number 0 as an ex
    ```
    (*Note:* in this case, we need to explicitly specify the particle type because it cannot be deduced from the expression)
 
-It is also possible to check various properties of expressions:
-   ```pycon
-   >>> a=Expr("\\sigma^z_0")
-   >>> a.is_real #If the corresponding matrix real
-   True
-   >>> a.is_hermitian #If the corresponding matrix hermitian
-   True
-   >>> a.is_identity #If the corresponding matrix identity
-   False
-   >>> a.number_sites #Number of sites in the expression
-   1
-   ```
 
 
 #### Operator algebra
@@ -339,7 +329,24 @@ Under the hood, lattice-symmetries rewrites all the expressions into the canonic
 >>> str(Expr("σ⁺₁ σ⁺₁"))
 "0.0 I"
 ```
+#### Complex expressions
+
+So far, we defined expressions only on a few number of sites. Here we will consider more complicated expressions.
 
+#### Properties
+
+It is also possible to check various properties of expressions:
+   ```pycon
+   >>> a=Expr("\\sigma^z_0")
+   >>> a.is_real #If the corresponding matrix real
+   True
+   >>> a.is_hermitian #If the corresponding matrix hermitian
+   True
+   >>> a.is_identity #If the corresponding matrix identity
+   False
+   >>> a.number_sites #Number of sites in the expression
+   1
+   ```
 
 ### Operators
 

From bfae425216503390fdb69c2113eb1c04c772784b Mon Sep 17 00:00:00 2001
From: asjosik1991 <asjosik@LAPTOP-F94RAJDK.localdomain>
Date: Mon, 25 Mar 2024 11:14:31 +0100
Subject: [PATCH 24/67] update tutorial

---
 python/tutorial.md | 134 +++++++++++++++++++++++++++++++++++++++------
 1 file changed, 118 insertions(+), 16 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index dab2c6b..33b54b0 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -296,7 +296,7 @@ Furthermore, expressions can also be multiplied by scalars from the left using t
 For example, here are a few ways to write down the Heisenberg interaction between sites 0 and 1
 
 $$
-\mathbf{S}_i \cdot \mathbf{S}_j = S^x_0 S^x_1 + S^y_0 S^y_1 + S^z_0 S^z_1
+\mathbf{S}_0 \cdot \mathbf{S}_1 = S^x_0 S^x_1 + S^y_0 S^y_1 + S^z_0 S^z_1
 $$
 
 ```pycon
@@ -331,33 +331,101 @@ Under the hood, lattice-symmetries rewrites all the expressions into the canonic
 ```
 #### Complex expressions
 
-So far, we defined expressions only on a few number of sites. Here we will consider more complicated expressions.
+So far, we defined expressions only on a few number of sites, which can be indeed easily written explicitly. Here we will consider more complicated expressions, which require other techniques.
+One of the ways is to use the function `replace_indices`. For example, we can construct the sum of $\sigma^z$ operators.
+
+```pycon
+import operator #import relevant methods
+from functools import reduce
+
+expr1 = ls.Expr("σᶻ₁") #define an elementary expression
+many_exprs = [expr1.replace_indices({1: i}) for i in range(4)] #we apply the permutation of vertices and make the corresponding array
+expr2 = reduce(operator.add, many_exprs) #sum all elementary operations
+print(expr2)
+>>>
+σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
+```
+
+Another way is to define an expression on a graph defined by edges:
+
+```pycon
+edges = [(i,) for i in range(4)] #define graph of 4 vertices 
+expr = ls.Expr("σᶻ₁", sites=edges) #the elementary expression is applied for all vertices
+>>>
+σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
+```
+One can see that the results are the same.
 
 #### Properties
 
+One can make standart operations on expressions, such make an adjoint, as well
+
+```pycon
+a = ls.Expr("c†₀ c₁")
+b=a.adjoint() #Makes an adjoint expression
+```
+
 It is also possible to check various properties of expressions:
-   ```pycon
-   >>> a=Expr("\\sigma^z_0")
-   >>> a.is_real #If the corresponding matrix real
-   True
-   >>> a.is_hermitian #If the corresponding matrix hermitian
-   True
-   >>> a.is_identity #If the corresponding matrix identity
-   False
-   >>> a.number_sites #Number of sites in the expression
-   1
-   ```
+```pycon
+>>> a=Expr("\\sigma^z_0")
+>>> a.is_real #If the corresponding matrix real
+True
+>>> a.is_hermitian #If the corresponding matrix hermitian
+True
+>>> a.is_identity #If the corresponding matrix identity
+False
+>>> a.number_sites #Number of sites in the expression
+1
+```
 
 ### Operators
 
-Based on an expression and a basis, we can build the corresponding operator acting on the chosen basis. Let's consider an example.
+Based on an expression and a basis, we can build the corresponding operator acting on the chosen basis. Let's at first consider the Hubbard model on two sites without spins.
+One of the ways to contruct an operator is to create an expression and then create the operator:
+
+```pycon
+op_expr = ls.Expr("-c†₀ c₁-c†₁ c₀+2.0 n₀ n₀+2.0 n₁ n₁")
+basis=ls.SpinlessFermionBasis(2)
+basis.build() #It is nessecary to build the basis before creating an operator
+opr=ls.Operator(basis, op_expr) #Here we create the operator
+```
+Another way is to create the expression for each small expression and then add them to each other:
+```pycon
+op_expr = ls.Expr("-c†₀ c₁")
+op_expr -= ls.Expr("c†₁ c₀")
+op_expr += 2*ls.Expr("n₀ n₀")
+op_expr += 2*ls.Expr("n₁ n₁")
+opr=ls.Operator(basis, op_expr)
+opr.shape() #We can check the shape of the operator
+>>>
+(4,4)
+```
 
+We can make a multiplication of the operator on a vector using standart python notation.
 
+```pycon
+vec1=np.array([1.0,0,0,0]) #One should work with floats
+opr@vec1
+>>>
+[0. 0. 0. 0.]
+
+vec2=np.array([0,1.0,0,0])
+opr@vec2
+>>>
+[ 0.  2. -1.  0.]
+
+vec3=np.array([1.0,2.0,0,0])
+opr@vec3
+>>>
+[ 0.  4. -2.  0.]
+```
+From this operations one can build the matrix form of an operator.
 
 ### Symmetries
 
 The full power of `lattice_symmetries` manifests if one use symmetries when constructing 
-symmetry adapted basis and linear operators acting on teh corresponding Hilbert space.
+symmetry adapted basis and linear operators acting on teh corresponding Hilbert space. 
+The symmetries can be constructed with the help of expressions as well. 
 
 ## Examples
 
@@ -365,7 +433,41 @@ Here we will take a look at different examples of `lattice_symmetries` applicati
 
 ### Simple ED
 
-Now we will consider the simplest example of exact diagonalization.
+We will strat with the simplest example of exact diagonalization. We will consider Heisenberg chain on 10 sites and diagonalization with the help of symmetries.
+For that, we combine all the considerd methods.
+At first we create the expression for Heisenberg chain on 10 sites using the methods from the section [Complex expressions](#Complex-expressions).
+
+```pycon
+# Constructing the basis
+    basis = ls.SpinBasis(
+        number_spins=number_spins,
+        # NOTE: we don't actually need to specify hamming_weight when spin_inversion
+        # is set. The library will guess that hamming_weight = number_spins / 2.
+        spin_inversion=-1,
+        symmetries=symmetries,
+    )
+    basis.build()  # Build the list of representatives, we need it since we're doing ED
+
+# Constructinf the expression
+edges = [(i, (i + 1) % number_spins) for i in range(number_spins)]
+expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁", sites=edges)
+print("Expression:", expr)
+
+# Alternatively, we can create expr in an algebraic way:
+two_site_expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁")
+many_exprs = [two_site_expr.replace_indices({0: i, 1: j}) for (i, j) in edges]
+expr2 = reduce(operator.add, many_exprs)
+assert expr == expr2 #we check that the expressions are equal
+
+# Construct the Hamiltonian
+hamiltonian = ls.Operator(basis, expr)
+
+
+# Diagonalize the Hamiltonian using ARPACK
+eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA")
+print("Ground state energy is {}".format(eigenvalues[0]))
+assert np.isclose(eigenvalues[0], -18.06178542)
+```
 
 ### More complicated ED
 

From c88e1d521bc813311f340df534f5a16c02b7f898 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Wed, 10 Apr 2024 15:41:52 +0200
Subject: [PATCH 25/67] Update tutorial.md

---
 python/tutorial.md | 403 +++++++++++++++++++++++++++------------------
 1 file changed, 243 insertions(+), 160 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 33b54b0..3ead8b2 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -4,17 +4,21 @@
 
 - [Installing](#Installing)
 - [Basic concepts and functions](#Basic-concepts-and-functions)
-    - [Basis](#Basis)
-        - [Spin basis](#Spin-basis)
-        - [Spinless Fermionic basis](#Spinless-fermionic-basis)
-        - [Spinful Fermionic basis](#Spinful-fermionic-basis)
     - [Expressions](#Expressions)
         - [Primitive operators](#Primitive-operators)
         - [Operator algebra](#Operator-algebra)
         - [Complex expressions](#Complex-expressions)
         - [Properties](#Properties)
+	- [Basis](#Basis)
+        - [Spin basis](#Spin-basis)
+        - [Spinless Fermionic basis](#Spinless-fermionic-basis)
+        - [Spinful Fermionic basis](#Spinful-fermionic-basis)
+		- [Basis from Expressions](#Basis-from-Expressions)
     - [Operators](#Operators)
     - [Symmetries](#Symmetries)
+		- [Symmetries as Permutations] (Symmetries-as-Permutations)
+		- [Symmetries from Expressions](Symmetries-from-Expressions)
+		- [Symmetry adapted basis](Symmetry-adapted-basis)
 - [Examples](#Examples)
     - [Simple ED](#Simple-ED)
     - [More complicated ED](#More-complicated-ED)
@@ -72,160 +76,18 @@ and takes into account system symmetries.
 The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
 such as time evolution or exact diagonalization of many-body Hamiltonians.
 
-The basic objects upon which the machinery is built are `Basis`, `Expression`, `Operator`, and `Symmetry`.
-- `Basis` object stores a basis of a many-body Hilbert space consisting of spins or fermions with appropriate symmetries.
-Each basis state is represented as a sequence of 0s and 1s (i.e. a sequence of bits), which can be interpreted as a binary number.
+The basic objects upon which the machinery is built are `Basis`, `Expression`, `Operator`, and `Symmetries`.
+The `Symmetries` is not separate class, however, they lie in the core of `lattice_symmetries`, so we will consider them separately.
 - `Expression` object is a nice way to work with symbolic representations of operators. It is possible to sum expressions, and multiply them by numbers and each other.
 `Expression` allows not to think about an explicit matrix representation of an operator, and the user can work directly with analytical formulae. 
+- `Basis` object stores a basis of a many-body Hilbert space consisting of spins or fermions with appropriate symmetries.
+Each basis state is represented as a sequence of 0s and 1s (i.e. a sequence of bits), which can be interpreted as a binary number.
 - `Operator` object is an actual operator that can act on individual basis states and their linear combinations. 
-- `Symmetry` is a method to work with symmetry adapted basises.
+- `Symmetries` is a method to work with symmetry adapted basises.
 If an operator has symmmetries, it is useful to work in symmetry-adapted basis, since it can drastically decrease the dimension of the Hilbert space.
 
 Now we will take a look at basic functions and methods for these objects. 
 
-### Basis
-#### Spin basis
-Let's look at simple examples; at first we will not consider additional symmetries.
-
-The simplest example would be a spin basis:
-
-```sh
-import lattice_symmetries as ls
- 
-basis = ls.SpinBasis(3) #We create basis consisting of three $\frac{1}{2}$-spins (each spin can be |0⟩ or |1⟩) 
-basis.build() #build basis
-print(basis.index(basis.states)) #Print array of indices of basis states
-print(basis.number_states) #Print total number of states in the basis
-for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
-    print(basis.states[i], basis.state_to_string(basis.states[i]))  
-```
-
-The result looks like:
-```sh
-[0 1 2 3 4 5 6 7]
-8
-0 |000⟩
-1 |001⟩
-2 |010⟩
-3 |011⟩
-4 |100⟩
-5 |101⟩
-6 |110⟩
-7 |111⟩
-```
-The basis states are equal to their indices and binary representations, as they should be.
-
-We can also consider only a part of the whole Hilbert space with a given number of spin ups (i.e. hamming weight of binary representations).
-We can specify it as follows:
-```sh
-basis = ls.SpinBasis(3, hamming_weight=2) #We want the subspace with only 2 spins up 
-basis.build()
-for i in range(basis.number_states): #Print the states in the basis
-    print(basis.states[i], basis.state_to_string(basis.states[i]))  
-```
-The output is:
-```sh
-3 |011⟩
-5 |101⟩
-6 |110⟩
-```
-We see that basis states include only spins with 2 spins up. It is also interesting to note that in this case a basis state is not equal to its index.
-
-Sometimes, our system has spin inversion symmetry, and we can additionally shorten our basis. In this case, we can specify it as follows:
-```sh
-basis = ls.SpinBasis(4, spin_inversion=-1) #Spin inversion is present. It is not nessecary to specify hamming weight here, since it is fixed by spin inversion symmetry. 
-basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i]))  
-```
-We have the basis states:
-```sh
-3 |011⟩
-5 |101⟩
-6 |110⟩
-```
-
-#### Spinless Fermionic basis
-We can also consider the basis of fermions without spins. The basis states are stored as integers as for spin basis. However the binary representation has a second quantization interpretation.
-Each basis state is given by the sequence of 0s and 1s, where 1 means a fermion on the corresponding site, and 0 means that the site is vacant.
-Let's consider the simplest example of fermions on two sites:
-```sh
-basis = ls.SpinlessFermionBasis(2) #We create fermionic basis on 2 sites
-basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i])) 
-```
-which gives:
-```sh
-0 |00⟩
-1 |01⟩
-2 |10⟩
-3 |11⟩
-```
-as one would expect.
-
-We can specify the number of particles as well:
-```sh
-basis = ls.SpinlessFermionBasis(4, number_particles=3) #We create fermionic basis on 4 sites with only 3 fermions
-basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i])) 
-```
-which gives:
-```sh
-7 |0111⟩
-11 |1011⟩
-13 |1101⟩
-14 |1110⟩
-```
-We can see that the basis consists of states with three fermions.
-
-#### Spinful Fermionic basis
-
-The last case includes fermions with spin. The binary representations of basis states can be read as a pair of (fermions with spin up on a lattice, fermions with spin down on a lattice).
-We can create a basis of spinful fermions as follows:
-```sh
-basis = ls.SpinfulFermionBasis(2) #We create basis of spinful fermions on 2 sites
-basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i])) 
-```
-which gives:
-```sh
-0 |00⟩|00⟩
-1 |00⟩|01⟩
-2 |00⟩|10⟩
-3 |00⟩|11⟩
-4 |01⟩|00⟩
-5 |01⟩|01⟩
-6 |01⟩|10⟩
-7 |01⟩|11⟩
-8 |10⟩|00⟩
-9 |10⟩|01⟩
-10 |10⟩|10⟩
-11 |10⟩|11⟩
-12 |11⟩|00⟩
-13 |11⟩|01⟩
-14 |11⟩|10⟩
-15 |11⟩|11⟩
-```
-We see that binary representation now means the second quantization of fermions with two spins.
-
-As before, we can specify a sector of the total Hilbert space with a given number of fermions with spin down and spin up:
-
-```sh
-basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins up and down (N_up, N_down)=(2,1)
-basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i])) 
-```
-which gives:
-```sh
-7 |01⟩|11⟩
-11 |10⟩|11⟩
-```
-as expected.
-
 ### Expressions
 
 Expressions are an easy way to work with operators (for example, Hamiltonians) on a symbolic level using second quantization formalism.
@@ -338,7 +200,7 @@ One of the ways is to use the function `replace_indices`. For example, we can co
 import operator #import relevant methods
 from functools import reduce
 
-expr1 = ls.Expr("σᶻ₁") #define an elementary expression
+expr1 = Expr("σᶻ₁") #define an elementary expression
 many_exprs = [expr1.replace_indices({1: i}) for i in range(4)] #we apply the permutation of vertices and make the corresponding array
 expr2 = reduce(operator.add, many_exprs) #sum all elementary operations
 print(expr2)
@@ -350,18 +212,39 @@ Another way is to define an expression on a graph defined by edges:
 
 ```pycon
 edges = [(i,) for i in range(4)] #define graph of 4 vertices 
-expr = ls.Expr("σᶻ₁", sites=edges) #the elementary expression is applied for all vertices
+expr = Expr("σᶻ₁", sites=edges) #the elementary expression is applied for all vertices
 >>>
 σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
 ```
+
+One can use the function `on` for the same effect:
+```
+e=Expr("σᶻ₁")
+expt=e.on(edges)
+>>>
+σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
+```
+
+It is also possible to use `iGraph`:
+```
+import igraph as ig
+
+e=Expr("σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀")
+expt=e.on(ig.Graph.Lattice(dim=[2, 3], circular=True))
+>>>
+σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
+```
+
+
+
 One can see that the results are the same.
 
 #### Properties
 
-One can make standart operations on expressions, such make an adjoint, as well
+One can make standard operations on expressions, such make an adjoint, as well
 
 ```pycon
-a = ls.Expr("c†₀ c₁")
+a = Expr("c†₀ c₁")
 b=a.adjoint() #Makes an adjoint expression
 ```
 
@@ -378,6 +261,156 @@ False
 1
 ```
 
+### Basis
+
+#### Spin basis
+Let's look at simple examples; at first we will not consider additional symmetries.
+
+The simplest example would be a spin basis:
+
+```sh
+import lattice_symmetries as ls
+ 
+basis = ls.SpinBasis(3) #We create basis consisting of three $\frac{1}{2}$-spins (each spin can be |0⟩ or |1⟩) 
+basis.build() #build basis
+print(basis.index(basis.states)) #Print array of indices of basis states
+print(basis.number_states) #Print total number of states in the basis
+for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
+    print(basis.states[i], basis.state_to_string(basis.states[i]))  
+```
+
+The result looks like:
+```sh
+[0 1 2 3 4 5 6 7]
+8
+0 |000⟩
+1 |001⟩
+2 |010⟩
+3 |011⟩
+4 |100⟩
+5 |101⟩
+6 |110⟩
+7 |111⟩
+```
+The basis states are equal to their indices and binary representations, as they should be.
+
+We can also consider only a part of the whole Hilbert space with a given number of spin ups (i.e. hamming weight of binary representations).
+We can specify it as follows:
+```sh
+basis = ls.SpinBasis(3, hamming_weight=2) #We want the subspace with only 2 spins up 
+basis.build()
+for i in range(basis.number_states): #Print the states in the basis
+    print(basis.states[i], basis.state_to_string(basis.states[i]))  
+```
+The output is:
+```sh
+3 |011⟩
+5 |101⟩
+6 |110⟩
+```
+We see that basis states include only spins with 2 spins up. It is also interesting to note that in this case a basis state is not equal to its index.
+
+Sometimes, our system has spin inversion symmetry, and we can additionally shorten our basis. In this case, we can specify it as follows:
+```sh
+basis = ls.SpinBasis(4, spin_inversion=-1) #Spin inversion is present. It is not nessecary to specify hamming weight here, since it is fixed by spin inversion symmetry. 
+basis.build()
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i]))  
+```
+We have the basis states:
+```sh
+3 |011⟩
+5 |101⟩
+6 |110⟩
+```
+
+#### Spinless Fermionic basis
+We can also consider the basis of fermions without spins. The basis states are stored as integers as for spin basis. However the binary representation has a second quantization interpretation.
+Each basis state is given by the sequence of 0s and 1s, where 1 means a fermion on the corresponding site, and 0 means that the site is vacant.
+Let's consider the simplest example of fermions on two sites:
+```sh
+basis = ls.SpinlessFermionBasis(2) #We create fermionic basis on 2 sites
+basis.build()
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+```
+which gives:
+```sh
+0 |00⟩
+1 |01⟩
+2 |10⟩
+3 |11⟩
+```
+as one would expect.
+
+We can specify the number of particles as well:
+```sh
+basis = ls.SpinlessFermionBasis(4, number_particles=3) #We create fermionic basis on 4 sites with only 3 fermions
+basis.build()
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+```
+which gives:
+```sh
+7 |0111⟩
+11 |1011⟩
+13 |1101⟩
+14 |1110⟩
+```
+We can see that the basis consists of states with three fermions.
+
+#### Spinful Fermionic basis
+
+The last case includes fermions with spin. The binary representations of basis states can be read as a pair of (fermions with spin up on a lattice, fermions with spin down on a lattice).
+We can create a basis of spinful fermions as follows:
+```sh
+basis = ls.SpinfulFermionBasis(2) #We create basis of spinful fermions on 2 sites
+basis.build()
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+```
+which gives:
+```sh
+0 |00⟩|00⟩
+1 |00⟩|01⟩
+2 |00⟩|10⟩
+3 |00⟩|11⟩
+4 |01⟩|00⟩
+5 |01⟩|01⟩
+6 |01⟩|10⟩
+7 |01⟩|11⟩
+8 |10⟩|00⟩
+9 |10⟩|01⟩
+10 |10⟩|10⟩
+11 |10⟩|11⟩
+12 |11⟩|00⟩
+13 |11⟩|01⟩
+14 |11⟩|10⟩
+15 |11⟩|11⟩
+```
+We see that binary representation now means the second quantization of fermions with two spins.
+
+As before, we can specify a sector of the total Hilbert space with a given number of fermions with spin down and spin up:
+
+```sh
+basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins up and down (N_up, N_down)=(2,1)
+basis.build()
+for i in range(basis.number_states):
+    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+```
+which gives:
+```sh
+7 |01⟩|11⟩
+11 |10⟩|11⟩
+```
+as expected.
+
+####Basis from expressions
+
+Before we created basises explicitly, however, there is a way to construct basis directly from expressions.
+EXAMPLE
+
+
 ### Operators
 
 Based on an expression and a basis, we can build the corresponding operator acting on the chosen basis. Let's at first consider the Hubbard model on two sites without spins.
@@ -421,11 +454,62 @@ opr@vec3
 ```
 From this operations one can build the matrix form of an operator.
 
+It is also possible to apply an operator directly to basis vectors.
+APPLICATION
+
 ### Symmetries
 
+#### Symmetries as Permutations
+
 The full power of `lattice_symmetries` manifests if one use symmetries when constructing 
-symmetry adapted basis and linear operators acting on teh corresponding Hilbert space. 
-The symmetries can be constructed with the help of expressions as well. 
+symmetry adapted basis and linear operators acting on the corresponding Hilbert space. 
+The symmetries are constructed with the help of expressions as well, and are represented as a permutation group of indices (realized as simpy PermutationGroup.)
+
+Let's take a look:
+```
+n=4
+translation = ls.Permutation([(1 + i) % n for i in range(n)])
+
+```
+
+####Symmetries from Expressions
+
+
+Since later we will need the characters to construct symmetry adapted basises, there are two possibilities.
+One can find all permutation symmetries of an expression (this option can be used to study specific sectors):
+```
+from sympy.combinatorics import Permutation, PermutationGroup
+from sympy.core import Rational
+
+ #let's use a simple chain of 4 spins, which has the symmetry group D_4
+symmetries=expr.permutation_group()
+>>>
+ #The result is a sympy PermutationGroup
+```
+
+or the maximum abelian subgroup of the symmetry group:
+```
+symmetries=expr.abelian_permutation_group()
+>>>
+
+```
+
+####Symmetry-adapted basis
+
+In order to make calculations with the help of symmetries, we need to construct a symmetry adapted basis.
+The basis is constructed with the help of one dimensional representations of the symmetry group of the Hamiltonian,
+or in other words, the symmetry (permutation) group of the corresponding expression.
+
+The simplest example would be:
+```
+n=4
+translation = ls.Permutation([(1 + i) % n for i in range(n)]) #we consider one-dimensional translations as before
+b = ls.SpinBasis(number_spins=n, symmetries=[(translation, ls.Rational(k, n))])
+b.build()
+```
+
+An `Operator` for symmetry adapted basis can be build in the same way as without symmetries.
+
 
 ## Examples
 
@@ -448,7 +532,7 @@ At first we create the expression for Heisenberg chain on 10 sites using the met
     )
     basis.build()  # Build the list of representatives, we need it since we're doing ED
 
-# Constructinf the expression
+# Constructing the expression
 edges = [(i, (i + 1) % number_spins) for i in range(number_spins)]
 expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁", sites=edges)
 print("Expression:", expr)
@@ -462,8 +546,7 @@ assert expr == expr2 #we check that the expressions are equal
 # Construct the Hamiltonian
 hamiltonian = ls.Operator(basis, expr)
 
-
-# Diagonalize the Hamiltonian using ARPACK
+# Diagonalize the Hamiltonian using ARPACK. The fast matrix-vector multiplication of `lattice-symmetries` will be used
 eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA")
 print("Ground state energy is {}".format(eigenvalues[0]))
 assert np.isclose(eigenvalues[0], -18.06178542)

From 20222defc3f26355355c764676ccdc7a21be4e8d Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Wed, 10 Apr 2024 15:47:23 +0200
Subject: [PATCH 26/67] Update tutorial.md

---
 python/tutorial.md | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 3ead8b2..f86fe51 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -15,10 +15,10 @@
         - [Spinful Fermionic basis](#Spinful-fermionic-basis)
 		- [Basis from Expressions](#Basis-from-Expressions)
     - [Operators](#Operators)
-    - [Symmetries](#Symmetries)
-		- [Symmetries as Permutations] (Symmetries-as-Permutations)
-		- [Symmetries from Expressions](Symmetries-from-Expressions)
-		- [Symmetry adapted basis](Symmetry-adapted-basis)
+    - [Symmetry](#Symmetry)
+		- [Symmetry as Permutations](#Symmetries-as-Permutations)
+		- [Symmetry from Expressions](#Symmetries-from-Expressions)
+		- [Symmetry adapted basis](#Symmetry-adapted-basis)
 - [Examples](#Examples)
     - [Simple ED](#Simple-ED)
     - [More complicated ED](#More-complicated-ED)
@@ -76,14 +76,14 @@ and takes into account system symmetries.
 The `lattice_symmetries` implements fast matrix-vector multiplication that can be applied to various problems, 
 such as time evolution or exact diagonalization of many-body Hamiltonians.
 
-The basic objects upon which the machinery is built are `Basis`, `Expression`, `Operator`, and `Symmetries`.
-The `Symmetries` is not separate class, however, they lie in the core of `lattice_symmetries`, so we will consider them separately.
+The basic objects upon which the machinery is built are `Basis`, `Expression`, `Operator`, and `Symmetry`.
+The `Symmetries` is not a separate class, however, it lies in the core of `lattice_symmetries`, so we will consider it separately.
 - `Expression` object is a nice way to work with symbolic representations of operators. It is possible to sum expressions, and multiply them by numbers and each other.
 `Expression` allows not to think about an explicit matrix representation of an operator, and the user can work directly with analytical formulae. 
 - `Basis` object stores a basis of a many-body Hilbert space consisting of spins or fermions with appropriate symmetries.
 Each basis state is represented as a sequence of 0s and 1s (i.e. a sequence of bits), which can be interpreted as a binary number.
 - `Operator` object is an actual operator that can act on individual basis states and their linear combinations. 
-- `Symmetries` is a method to work with symmetry adapted basises.
+- `Symmetry` is a method to work with symmetry adapted basises.
 If an operator has symmmetries, it is useful to work in symmetry-adapted basis, since it can drastically decrease the dimension of the Hilbert space.
 
 Now we will take a look at basic functions and methods for these objects. 
@@ -457,9 +457,9 @@ From this operations one can build the matrix form of an operator.
 It is also possible to apply an operator directly to basis vectors.
 APPLICATION
 
-### Symmetries
+### Symmetry
 
-#### Symmetries as Permutations
+#### Symmetry as Permutations
 
 The full power of `lattice_symmetries` manifests if one use symmetries when constructing 
 symmetry adapted basis and linear operators acting on the corresponding Hilbert space. 
@@ -472,7 +472,7 @@ translation = ls.Permutation([(1 + i) % n for i in range(n)])
 
 ```
 
-####Symmetries from Expressions
+#### Symmetry from Expressions
 
 
 Since later we will need the characters to construct symmetry adapted basises, there are two possibilities.
@@ -494,7 +494,7 @@ symmetries=expr.abelian_permutation_group()
 
 ```
 
-####Symmetry-adapted basis
+#### Symmetry-adapted basis
 
 In order to make calculations with the help of symmetries, we need to construct a symmetry adapted basis.
 The basis is constructed with the help of one dimensional representations of the symmetry group of the Hamiltonian,

From b349cd69b258a1022627a0a4b91f3f1fc010e4e8 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Wed, 10 Apr 2024 15:48:29 +0200
Subject: [PATCH 27/67] Update tutorial.md

---
 python/tutorial.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index f86fe51..bea50f8 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -13,11 +13,11 @@
         - [Spin basis](#Spin-basis)
         - [Spinless Fermionic basis](#Spinless-fermionic-basis)
         - [Spinful Fermionic basis](#Spinful-fermionic-basis)
-		- [Basis from Expressions](#Basis-from-Expressions)
+		- [Basis from Expressions](#Basis-from-expressions)
     - [Operators](#Operators)
     - [Symmetry](#Symmetry)
-		- [Symmetry as Permutations](#Symmetries-as-Permutations)
-		- [Symmetry from Expressions](#Symmetries-from-Expressions)
+		- [Symmetry as Permutations](#Symmetries-as-permutations)
+		- [Symmetry from Expressions](#Symmetries-from-expressions)
 		- [Symmetry adapted basis](#Symmetry-adapted-basis)
 - [Examples](#Examples)
     - [Simple ED](#Simple-ED)

From c844478d64d696891feeb2530e4fbce37612b9a8 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 14:44:32 +0200
Subject: [PATCH 28/67] Update tutorial.md

---
 python/tutorial.md | 119 ++++++++++++++++++++++++++++-----------------
 1 file changed, 75 insertions(+), 44 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index bea50f8..8c989cb 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -5,8 +5,8 @@
 - [Installing](#Installing)
 - [Basic concepts and functions](#Basic-concepts-and-functions)
     - [Expressions](#Expressions)
-        - [Primitive operators](#Primitive-operators)
-        - [Operator algebra](#Operator-algebra)
+        - [Primitive expressions](#Primitive-expressions)
+        - [Algebra of expressions](#Algebra-of-expressions)
         - [Complex expressions](#Complex-expressions)
         - [Properties](#Properties)
 	- [Basis](#Basis)
@@ -16,8 +16,8 @@
 		- [Basis from Expressions](#Basis-from-expressions)
     - [Operators](#Operators)
     - [Symmetry](#Symmetry)
-		- [Symmetry as Permutations](#Symmetries-as-permutations)
-		- [Symmetry from Expressions](#Symmetries-from-expressions)
+		- [Symmetry as Permutations](#Symmetry-as-permutations)
+		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry adapted basis](#Symmetry-adapted-basis)
 - [Examples](#Examples)
     - [Simple ED](#Simple-ED)
@@ -91,21 +91,21 @@ Now we will take a look at basic functions and methods for these objects.
 ### Expressions
 
 Expressions are an easy way to work with operators (for example, Hamiltonians) on a symbolic level using second quantization formalism.
-This means that you can use primitive operators such as $\sigma^x$, $\sigma^y$, and $\sigma^z$ to build expressions for your Hamiltonian and observables.
+This means that you can use primitive symbols of well-known operators such as $\sigma^x$, $\sigma^y$, and $\sigma^z$ to build expressions for your Hamiltonian and observables.
 It is also possible to sum different expressions and multiply them to each other to compose more complicated expressions.
 Let's consider at first the simplest examples. 
 
-#### Primitive operators
+#### Primitive expressions
 
 At first we need to import `Expr` from lattice-symmetries:
 ```sh
 from lattice_symmetries import Expr
 ```
-Now we will consider primitives operators defined on site with number 0 as an example, but you can also construct primitive operators residing on other lattice sites:
+Now we will consider primitive symbols defined on site with number 0 as an example, but you can also construct them residing on other lattice sites:
 
  - $\sigma^x$ and $S^x$:
    ```pycon
-   >>> Expr("σˣ₀") == Expr("\\sigma^x_0") #It is possible to use different notations for Expr for primitive operators.
+   >>> Expr("σˣ₀") == Expr("\\sigma^x_0") #It is possible to use different notations for Expr for primitives.
    #Here we check that they agree. Index 0 means the index of the corresponding site.
    True
    >>> Expr("Sˣ₀") == 0.5 * Expr("σˣ₀") #We check that Sˣ₀ is the same as 0.5*σˣ₀.
@@ -150,9 +150,9 @@ Now we will consider primitives operators defined on site with number 0 as an ex
 
 
 
-#### Operator algebra
+#### Algebra of expressions
 
-Primitives can be combined using the `+`, `-`, and `*` operations to build more complex operators.
+Primitives can be combined using the `+`, `-`, and `*` operations to build more complex expressions.
 Furthermore, expressions can also be multiplied by scalars from the left using the `*` operator.
 
 For example, here are a few ways to write down the Heisenberg interaction between sites 0 and 1
@@ -194,25 +194,25 @@ Under the hood, lattice-symmetries rewrites all the expressions into the canonic
 #### Complex expressions
 
 So far, we defined expressions only on a few number of sites, which can be indeed easily written explicitly. Here we will consider more complicated expressions, which require other techniques.
-One of the ways is to use the function `replace_indices`. For example, we can construct the sum of $\sigma^z$ operators.
+One of the ways is to use the function `replace_indices`. For example, we can construct the sum of $\sigma^z$s defined on different sites: 
 
 ```pycon
-import operator #import relevant methods
+import operator #Import relevant methods
 from functools import reduce
 
-expr1 = Expr("σᶻ₁") #define an elementary expression
-many_exprs = [expr1.replace_indices({1: i}) for i in range(4)] #we apply the permutation of vertices and make the corresponding array
-expr2 = reduce(operator.add, many_exprs) #sum all elementary operations
+expr1 = Expr("σᶻ₁") #Define an elementary expression
+many_exprs = [expr1.replace_indices({1: i}) for i in range(4)] #We apply the permutation of vertices and make the corresponding array
+expr2 = reduce(operator.add, many_exprs) #Sum all elementary operations
 print(expr2)
 >>>
 σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
 ```
 
-Another way is to define an expression on a graph defined by edges:
+Another way is to define an expression on a (hyper)graph defined by edges:
 
 ```pycon
-edges = [(i,) for i in range(4)] #define graph of 4 vertices 
-expr = Expr("σᶻ₁", sites=edges) #the elementary expression is applied for all vertices
+edges = [(i,) for i in range(4)] #Define graph of 4 vertices. Edges consist only of one element, because the elementary expression is linear
+expr = Expr("σᶻ₁", sites=edges) #The elementary expression is applied for all vertices
 >>>
 σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
 ```
@@ -225,19 +225,17 @@ expt=e.on(edges)
 σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
 ```
 
-It is also possible to use `iGraph`:
+It is also possible to use `iGraph` to define an underlying (hyper)graph:
 ```
 import igraph as ig
 
-e=Expr("σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀")
-expt=e.on(ig.Graph.Lattice(dim=[2, 3], circular=True))
+e=Expr("σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀") #Here we have two-sites interaction, so we need an actual graph with edges consist of two vertices
+expt=e.on(ig.Graph.Lattice(dim=[2, 2])) # The graph is square 2x2
 >>>
-σᶻ₀ + σᶻ₁ + σᶻ₂ + σᶻ₃
+σ⁺₀ σ⁻₁ + σ⁺₀ σ⁻₂ + σ⁻₀ σ⁺₁ + σ⁻₀ σ⁺₂ + σ⁺₁ σ⁻₃ + σ⁻₁ σ⁺₃ + σ⁺₂ σ⁻₃ + σ⁻₂ σ⁺₃
 ```
-
-
-
-One can see that the results are the same.
+However, the method of defining a Hamiltonian on a graph works correctly only if the elementary expression (defined on one site) is homogeneous.
+It is also important that edges of a (hyper)graph should have the same number of vertices as the degree of the elementary expression.
 
 #### Properties
 
@@ -263,6 +261,9 @@ False
 
 ### Basis
 
+In this section we will consider generation of basises of various types.
+We will not consider symmetry adapted basises yet, we will do it in the section [Symmetry adapted basis](#Symmetry-adapted-basis).
+
 #### Spin basis
 Let's look at simple examples; at first we will not consider additional symmetries.
 
@@ -417,10 +418,12 @@ Based on an expression and a basis, we can build the corresponding operator acti
 One of the ways to contruct an operator is to create an expression and then create the operator:
 
 ```pycon
+import lattice_symmetries as ls
+
 op_expr = ls.Expr("-c†₀ c₁-c†₁ c₀+2.0 n₀ n₀+2.0 n₁ n₁")
 basis=ls.SpinlessFermionBasis(2)
 basis.build() #It is nessecary to build the basis before creating an operator
-opr=ls.Operator(basis, op_expr) #Here we create the operator
+opr=ls.Operator(op_expr, basis) #Here we create the operator
 ```
 Another way is to create the expression for each small expression and then add them to each other:
 ```pycon
@@ -428,7 +431,7 @@ op_expr = ls.Expr("-c†₀ c₁")
 op_expr -= ls.Expr("c†₁ c₀")
 op_expr += 2*ls.Expr("n₀ n₀")
 op_expr += 2*ls.Expr("n₁ n₁")
-opr=ls.Operator(basis, op_expr)
+opr=ls.Operator(op_expr, basis)
 opr.shape() #We can check the shape of the operator
 >>>
 (4,4)
@@ -463,7 +466,7 @@ APPLICATION
 
 The full power of `lattice_symmetries` manifests if one use symmetries when constructing 
 symmetry adapted basis and linear operators acting on the corresponding Hilbert space. 
-The symmetries are constructed with the help of expressions as well, and are represented as a permutation group of indices (realized as simpy PermutationGroup.)
+The symmetries are constructed with the help of expressions, and are represented as a permutation group of indices.
 
 Let's take a look:
 ```
@@ -494,7 +497,7 @@ symmetries=expr.abelian_permutation_group()
 
 ```
 
-#### Symmetry-adapted basis
+#### Symmetry adapted basis
 
 In order to make calculations with the help of symmetries, we need to construct a symmetry adapted basis.
 The basis is constructed with the help of one dimensional representations of the symmetry group of the Hamiltonian,
@@ -518,21 +521,33 @@ Here we will take a look at different examples of `lattice_symmetries` applicati
 ### Simple ED
 
 We will strat with the simplest example of exact diagonalization. We will consider Heisenberg chain on 10 sites and diagonalization with the help of symmetries.
-For that, we combine all the considerd methods.
-At first we create the expression for Heisenberg chain on 10 sites using the methods from the section [Complex expressions](#Complex-expressions).
+For that, we will combine methods described in the previous sections.
 
 ```pycon
+import lattice_symmetries as ls
+import numpy as np
+import scipy
+
+number_spins = 10  # System size
+
+# Constructing symmetries
+sites = np.arange(number_spins)
+# Momentum in x direction with phase φ=1/2 (i.e., with the eigenvalue λ=exp(-2πiφ)=-1)
+T = (ls.Permutation((sites + 1) % number_spins), ls.Rational(1, 2))
+# Parity with eigenvalue λ=-1
+P = (ls.Permutation(sites[::-1]), ls.Rational(1, 2))
+
 # Constructing the basis
-    basis = ls.SpinBasis(
-        number_spins=number_spins,
-        # NOTE: we don't actually need to specify hamming_weight when spin_inversion
-        # is set. The library will guess that hamming_weight = number_spins / 2.
-        spin_inversion=-1,
-        symmetries=symmetries,
-    )
-    basis.build()  # Build the list of representatives, we need it since we're doing ED
-
-# Constructing the expression
+basis = ls.SpinBasis(
+	number_spins=number_spins,
+	# NOTE: we don't actually need to specify hamming_weight when spin_inversion
+	# is set. The library will guess that hamming_weight = number_spins // 2.
+	spin_inversion=-1,
+	symmetries=[T, P], #Here we use the characters of the whole symmetry group (it's worth noting that, in general, translations and parity don't commute)
+)
+basis.build()  # Build the list of representatives, we need it since we're doing ED
+
+# Constructing the expression of the hamiltonian
 edges = [(i, (i + 1) % number_spins) for i in range(number_spins)]
 expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁", sites=edges)
 print("Expression:", expr)
@@ -544,12 +559,28 @@ expr2 = reduce(operator.add, many_exprs)
 assert expr == expr2 #we check that the expressions are equal
 
 # Construct the Hamiltonian
-hamiltonian = ls.Operator(basis, expr)
+hamiltonian = ls.Operator(expr, basis)
+
+# Diagonalize the Hamiltonian using ARPACK. The fast matrix-vector multiplication of `lattice-symmetries` will be used,
+#because hamiltonian is a lattice-symmetry operator
+eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA")
+print("Ground state energy is {}".format(eigenvalues[0]))
+assert np.isclose(eigenvalues[0], -18.06178542)
 
-# Diagonalize the Hamiltonian using ARPACK. The fast matrix-vector multiplication of `lattice-symmetries` will be used
+#Here we used a symmetry adapted basis, and therefore we were restricted to a specific sector of the total Hilbert space.
+#Since the considered system is relatively small, we can make the diagonalization in the full basis and check that we chose the right sector.
+
+new_basis = ls.SpinBasis(
+	number_spins=number_spins,
+	spin_inversion=-1,
+)
+new_basis.build()
+new_hamiltonian = ls.Operator(expr, new_basis)
 eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA")
 print("Ground state energy is {}".format(eigenvalues[0]))
 assert np.isclose(eigenvalues[0], -18.06178542)
+
+#Indeed, the outcome is the same within machine precision. 
 ```
 
 ### More complicated ED

From 3857627a7fdc5fb4deeceaf13c7cd218cb47666a Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 14:46:19 +0200
Subject: [PATCH 29/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 8c989cb..7549e64 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -84,7 +84,7 @@ The `Symmetries` is not a separate class, however, it lies in the core of `latti
 Each basis state is represented as a sequence of 0s and 1s (i.e. a sequence of bits), which can be interpreted as a binary number.
 - `Operator` object is an actual operator that can act on individual basis states and their linear combinations. 
 - `Symmetry` is a method to work with symmetry adapted basises.
-If an operator has symmmetries, it is useful to work in symmetry-adapted basis, since it can drastically decrease the dimension of the Hilbert space.
+If an operator has symmetries, it is useful to work in symmetry-adapted basis, since it can drastically decrease the dimension of the Hilbert space.
 
 Now we will take a look at basic functions and methods for these objects. 
 

From 72d0abfc4e3453e01d861d51287a90012f1db49e Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:16:28 +0200
Subject: [PATCH 30/67] Update tutorial.md

---
 python/tutorial.md | 82 ++++++++++++++++++++++++++++++++--------------
 1 file changed, 57 insertions(+), 25 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 7549e64..0b34d2c 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -16,7 +16,7 @@
 		- [Basis from Expressions](#Basis-from-expressions)
     - [Operators](#Operators)
     - [Symmetry](#Symmetry)
-		- [Symmetry as Permutations](#Symmetry-as-permutations)
+		- [Symmetry as Permutations](#Symmetries-as-permutations)
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry adapted basis](#Symmetry-adapted-basis)
 - [Examples](#Examples)
@@ -230,7 +230,7 @@ It is also possible to use `iGraph` to define an underlying (hyper)graph:
 import igraph as ig
 
 e=Expr("σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀") #Here we have two-sites interaction, so we need an actual graph with edges consist of two vertices
-expt=e.on(ig.Graph.Lattice(dim=[2, 2])) # The graph is square 2x2
+expr=e.on(ig.Graph.Lattice(dim=[2, 2])) # The graph is square 2x2
 >>>
 σ⁺₀ σ⁻₁ + σ⁺₀ σ⁻₂ + σ⁻₀ σ⁺₁ + σ⁻₀ σ⁺₂ + σ⁺₁ σ⁻₃ + σ⁻₁ σ⁺₃ + σ⁺₂ σ⁻₃ + σ⁻₂ σ⁺₃
 ```
@@ -406,10 +406,16 @@ which gives:
 ```
 as expected.
 
-####Basis from expressions
+####Basis from Expressions
 
-Before we created basises explicitly, however, there is a way to construct basis directly from expressions.
-EXAMPLE
+Before we created basises explicitly, however, there is a way to construct a basis directly from expressions.
+However, in this case, one can obtain only full Fock basis, without restriction on the number of particles.
+
+```
+expr=ls.Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁")
+basis=expr.full_basis()
+basis.build()
+```
 
 
 ### Operators
@@ -455,10 +461,6 @@ opr@vec3
 >>>
 [ 0.  4. -2.  0.]
 ```
-From this operations one can build the matrix form of an operator.
-
-It is also possible to apply an operator directly to basis vectors.
-APPLICATION
 
 ### Symmetry
 
@@ -466,37 +468,51 @@ APPLICATION
 
 The full power of `lattice_symmetries` manifests if one use symmetries when constructing 
 symmetry adapted basis and linear operators acting on the corresponding Hilbert space. 
-The symmetries are constructed with the help of expressions, and are represented as a permutation group of indices.
+The symmetries are constructed with the help of expressions, and are represented as a permutation group of indices;
+`lattice_symmetries` uses sympy to represent permutations, therefore one can take a look at sympy documentation for more details.
 
-Let's take a look:
+Let's take a look at a couple of examples:
+```
+p = ls.Permutation([1,2,3,0]) #This permutation shifts the indices, so that 0->1, 1->2, 2->3, 3->0
+>>> (0 1 2 3) #The same permutation written in the cycle representation
 ```
-n=4
-translation = ls.Permutation([(1 + i) % n for i in range(n)])
-
+```
+p = ls.Permutation([0,1,2,3]) #This is the identity permutation
+>>> (3) #No cycles, therefore identity. If the largest index doesn't move, it is shown in brackets, the format is used by sympy
+```
+```
+p=ls.Permutation([1,0,3,2]) #Exchange indices 0<->1 and 2<->3
+>>> (0 1)(2 3) #Two cycles
 ```
 
 #### Symmetry from Expressions
 
+In the previous section, we constructed the symmetries by hand, however, this can be done only for relatively easy and small systems.
+In most of the cases, the more powerful and convenient way is to ask `lattice_symmetries` to calculate the simmetries of a given expression. 
 
-Since later we will need the characters to construct symmetry adapted basises, there are two possibilities.
-One can find all permutation symmetries of an expression (this option can be used to study specific sectors):
-```
-from sympy.combinatorics import Permutation, PermutationGroup
-from sympy.core import Rational
+Since we need the characters to construct symmetry adapted basises, there are two possibilities.
+
+- One can find all permutation symmetries of an expression:
 
- #let's use a simple chain of 4 spins, which has the symmetry group D_4
-symmetries=expr.permutation_group()
+```
+ #let's use a simple chain of 3 spins, which has the symmetry group D_3
+e=ls.Expr("σ^x_0 σ^x_1")
+expr=e.on(ig.Graph.Lattice(dim=[3], circular=True)) # The periodic chain with 3 sites
+sym=expr.permutation_group()
 >>>
- #The result is a sympy PermutationGroup
+
 ```
+This option can be used to study specific sectors and does not cover the whole Hilbert space.
 
-or the maximum abelian subgroup of the symmetry group:
+- Another option is to find the maximum abelian subgroup of the symmetry group:
 ```
-symmetries=expr.abelian_permutation_group()
+ab_sym=expr.abelian_permutation_group()
 >>>
 
 ```
 
+In this case the sectors cover the whole Hilbert space.
+
 #### Symmetry adapted basis
 
 In order to make calculations with the help of symmetries, we need to construct a symmetry adapted basis.
@@ -505,7 +521,7 @@ or in other words, the symmetry (permutation) group of the corresponding express
 
 The simplest example would be:
 ```
-n=4
+n=3
 translation = ls.Permutation([(1 + i) % n for i in range(n)]) #we consider one-dimensional translations as before
 b = ls.SpinBasis(number_spins=n, symmetries=[(translation, ls.Rational(k, n))])
 b.build()
@@ -587,6 +603,22 @@ assert np.isclose(eigenvalues[0], -18.06178542)
 
 Now let's consider a more complicated example of ED.
 
+```pycon
+import lattice_symmetries as ls
+import numpy as np
+import scipy
+```
+
 ### Time evolution
 
 Another example of capabilities of `lattice_symmetries` is time evolution.
+In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
+$$
+e^{iHt}=e^{-i(E^{*}_g+aW')t}\Big[J_0(at)+\sum^{\infty}_{n=1}2(-i)^n J_n(at)T_n(H') \Big]
+$$
+where we rescale the original Hamiltonian $H$ with bandwidth $[E_g, E_s]$ to the Hamiltonian $H'$ with bandwidth $[-1+\epsilon,1-\epsilon]$, so that the series converges.
+The rescaling takes the form:
+$$
+H'=\frac{H-b}{a}
+$$
+where $a=1/(2-\epsilon)$, $b=(E_g+E_s)/2, and $\epsilon$ is a safety factor to make the spectrum of $H'$ liying within $[-1,1]$.
\ No newline at end of file

From f30fe6f553f380e501e18c02ed7aacb19c33bb16 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:20:43 +0200
Subject: [PATCH 31/67] Update tutorial.md

---
 python/tutorial.md | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 0b34d2c..c81e7c3 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -406,7 +406,7 @@ which gives:
 ```
 as expected.
 
-####Basis from Expressions
+#### Basis from Expressions
 
 Before we created basises explicitly, however, there is a way to construct a basis directly from expressions.
 However, in this case, one can obtain only full Fock basis, without restriction on the number of particles.
@@ -478,7 +478,8 @@ p = ls.Permutation([1,2,3,0]) #This permutation shifts the indices, so that 0->1
 ```
 ```
 p = ls.Permutation([0,1,2,3]) #This is the identity permutation
->>> (3) #No cycles, therefore identity. If the largest index doesn't move, it is shown in brackets, the format is used by sympy
+>>> (3) #No cycles, therefore identity. 
+#If the largest index doesn't move, it is shown in brackets, the format is used by sympy
 ```
 ```
 p=ls.Permutation([1,0,3,2]) #Exchange indices 0<->1 and 2<->3
@@ -613,12 +614,16 @@ import scipy
 
 Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
+
 $$
 e^{iHt}=e^{-i(E^{*}_g+aW')t}\Big[J_0(at)+\sum^{\infty}_{n=1}2(-i)^n J_n(at)T_n(H') \Big]
 $$
+
 where we rescale the original Hamiltonian $H$ with bandwidth $[E_g, E_s]$ to the Hamiltonian $H'$ with bandwidth $[-1+\epsilon,1-\epsilon]$, so that the series converges.
 The rescaling takes the form:
+
 $$
 H'=\frac{H-b}{a}
 $$
-where $a=1/(2-\epsilon)$, $b=(E_g+E_s)/2, and $\epsilon$ is a safety factor to make the spectrum of $H'$ liying within $[-1,1]$.
\ No newline at end of file
+
+where $a=1/(2-\epsilon)$, $b=(E_g+E_s)/2$, and $\epsilon$ is a safety factor to make the spectrum of $H'$ liying within $[-1,1]$.
\ No newline at end of file

From 0d5d78dbba6915d14dfd1bced58fdcba5d4a6f34 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:22:13 +0200
Subject: [PATCH 32/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index c81e7c3..55d970f 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}\Big[J_0(at)+\sum^{\infty}_{n=1}2(-i)^n J_n(at)T_n(H') \Big]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum^{\infty}_{n=1}2(-i)^n J_n(at)T_n(H')]
 $$
 
 where we rescale the original Hamiltonian $H$ with bandwidth $[E_g, E_s]$ to the Hamiltonian $H'$ with bandwidth $[-1+\epsilon,1-\epsilon]$, so that the series converges.

From b9896650579e163b2c33364faa09ed2d70dc8a33 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:23:35 +0200
Subject: [PATCH 33/67] Update tutorial.md

---
 python/tutorial.md | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 55d970f..43850d8 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -614,9 +614,12 @@ import scipy
 
 Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
+$$
+e^{iHt}=e^{-i(E^{*}_g+aW')t
+$$
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum^{\infty}_{n=1}2(-i)^n J_n(at)T_n(H')]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum^{\infty}_{n=1} 2(-i)^n J_n(at)T_n(H')]
 $$
 
 where we rescale the original Hamiltonian $H$ with bandwidth $[E_g, E_s]$ to the Hamiltonian $H'$ with bandwidth $[-1+\epsilon,1-\epsilon]$, so that the series converges.

From fe03c1656ee4e05a14675d1f15f3ea246e86e32d Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:23:52 +0200
Subject: [PATCH 34/67] Update tutorial.md

---
 python/tutorial.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/tutorial.md b/python/tutorial.md
index 43850d8..e33bec6 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -614,6 +614,7 @@ import scipy
 
 Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
+
 $$
 e^{iHt}=e^{-i(E^{*}_g+aW')t
 $$

From d6c96125639b6e7c11bc56af20aab9fc1f238bdd Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:24:32 +0200
Subject: [PATCH 35/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index e33bec6..7cd2c23 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t
+e^{iHt}=e^{-i(E^{*}_g+aW')t}
 $$
 
 $$

From 71e4cef4c26919baf1f7c85be262fd4dba0ac9a6 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:25:01 +0200
Subject: [PATCH 36/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 7cd2c23..6e9d26b 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)]
 $$
 
 $$

From b3dcb7ffc530dd9268630fbe7c55be8c4a150384 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:25:28 +0200
Subject: [PATCH 37/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 6e9d26b..bb919ef 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+2(-i)^n J_n(at)T_n(H')]
 $$
 
 $$

From cb0760dddb1e0a8b11301ec905f18340a88001b7 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:26:02 +0200
Subject: [PATCH 38/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index bb919ef..d25b12e 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+2(-i)^n J_n(at)T_n(H')]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum_{n=1}2(-i)^n J_n(at)T_n(H')]
 $$
 
 $$

From 8fd7ae57472ebf40171c5d0f8b3fc2398f18e69e Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:26:26 +0200
Subject: [PATCH 39/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index d25b12e..6c33ae7 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum_{n=1}2(-i)^n J_n(at)T_n(H')]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum 2(-i)^n J_n(at)T_n(H')]
 $$
 
 $$

From 5b4be76551e9e3a360d3cbcf009220585088c841 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:27:12 +0200
Subject: [PATCH 40/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 6c33ae7..601eea9 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum 2(-i)^n J_n(at)T_n(H')]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum^{\infty} 2(-i)^n J_n(at)T_n(H')]
 $$
 
 $$

From d6464e72c839f4ba804dafbb5cfc1182ae10d2db Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:27:37 +0200
Subject: [PATCH 41/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 601eea9..183a839 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum^{\infty} 2(-i)^n J_n(at)T_n(H')]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum_{n=1}^{\infty} 2(-i)^n J_n(at)T_n(H')]
 $$
 
 $$

From b6ec0c8979dd85b931ca602e80ee277b2e280745 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:29:49 +0200
Subject: [PATCH 42/67] Update tutorial.md

---
 python/tutorial.md | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 183a839..1084bfb 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,11 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum_{n=1}^{\infty} 2(-i)^n J_n(at)T_n(H')]
-$$
-
-$$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum^{\infty}_{n=1} 2(-i)^n J_n(at)T_n(H')]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum\limits_{n=1}^{\infty} 2(-i)^n J_n(at)T_n(H')]
 $$
 
 where we rescale the original Hamiltonian $H$ with bandwidth $[E_g, E_s]$ to the Hamiltonian $H'$ with bandwidth $[-1+\epsilon,1-\epsilon]$, so that the series converges.

From 1be41303a80859e5b115d08950438fff26e49b35 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:30:59 +0200
Subject: [PATCH 43/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 1084bfb..5ebe125 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -616,7 +616,7 @@ Another example of capabilities of `lattice_symmetries` is time evolution.
 In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum\limits_{n=1}^{\infty} 2(-i)^n J_n(at)T_n(H')]
+e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum\limits^{\infty} 2(-i)^n J_n(at)T_n(H')]
 $$
 
 where we rescale the original Hamiltonian $H$ with bandwidth $[E_g, E_s]$ to the Hamiltonian $H'$ with bandwidth $[-1+\epsilon,1-\epsilon]$, so that the series converges.

From 4b90875fccb4772a4d93bb50aa6d8f03bd18c354 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 17:37:41 +0200
Subject: [PATCH 44/67] Update tutorial.md

---
 python/tutorial.md | 66 ++++++++++++++++++----------------------------
 1 file changed, 25 insertions(+), 41 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 5ebe125..e64f018 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -98,7 +98,7 @@ Let's consider at first the simplest examples.
 #### Primitive expressions
 
 At first we need to import `Expr` from lattice-symmetries:
-```sh
+```pycon
 from lattice_symmetries import Expr
 ```
 Now we will consider primitive symbols defined on site with number 0 as an example, but you can also construct them residing on other lattice sites:
@@ -218,7 +218,7 @@ expr = Expr("σᶻ₁", sites=edges) #The elementary expression is applied for a
 ```
 
 One can use the function `on` for the same effect:
-```
+```pycon
 e=Expr("σᶻ₁")
 expt=e.on(edges)
 >>>
@@ -226,7 +226,7 @@ expt=e.on(edges)
 ```
 
 It is also possible to use `iGraph` to define an underlying (hyper)graph:
-```
+```pycon
 import igraph as ig
 
 e=Expr("σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀") #Here we have two-sites interaction, so we need an actual graph with edges consist of two vertices
@@ -269,7 +269,7 @@ Let's look at simple examples; at first we will not consider additional symmetri
 
 The simplest example would be a spin basis:
 
-```sh
+```pycon
 import lattice_symmetries as ls
  
 basis = ls.SpinBasis(3) #We create basis consisting of three $\frac{1}{2}$-spins (each spin can be |0⟩ or |1⟩) 
@@ -278,10 +278,7 @@ print(basis.index(basis.states)) #Print array of indices of basis states
 print(basis.number_states) #Print total number of states in the basis
 for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
     print(basis.states[i], basis.state_to_string(basis.states[i]))  
-```
-
-The result looks like:
-```sh
+>>>
 [0 1 2 3 4 5 6 7]
 8
 0 |000⟩
@@ -297,14 +294,12 @@ The basis states are equal to their indices and binary representations, as they
 
 We can also consider only a part of the whole Hilbert space with a given number of spin ups (i.e. hamming weight of binary representations).
 We can specify it as follows:
-```sh
+```pycon
 basis = ls.SpinBasis(3, hamming_weight=2) #We want the subspace with only 2 spins up 
 basis.build()
 for i in range(basis.number_states): #Print the states in the basis
     print(basis.states[i], basis.state_to_string(basis.states[i]))  
-```
-The output is:
-```sh
+>>>
 3 |011⟩
 5 |101⟩
 6 |110⟩
@@ -312,14 +307,12 @@ The output is:
 We see that basis states include only spins with 2 spins up. It is also interesting to note that in this case a basis state is not equal to its index.
 
 Sometimes, our system has spin inversion symmetry, and we can additionally shorten our basis. In this case, we can specify it as follows:
-```sh
+```pycon
 basis = ls.SpinBasis(4, spin_inversion=-1) #Spin inversion is present. It is not nessecary to specify hamming weight here, since it is fixed by spin inversion symmetry. 
 basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i]))  
-```
-We have the basis states:
-```sh
+>>>
 3 |011⟩
 5 |101⟩
 6 |110⟩
@@ -329,14 +322,12 @@ We have the basis states:
 We can also consider the basis of fermions without spins. The basis states are stored as integers as for spin basis. However the binary representation has a second quantization interpretation.
 Each basis state is given by the sequence of 0s and 1s, where 1 means a fermion on the corresponding site, and 0 means that the site is vacant.
 Let's consider the simplest example of fermions on two sites:
-```sh
+```pycon
 basis = ls.SpinlessFermionBasis(2) #We create fermionic basis on 2 sites
 basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 
-```
-which gives:
-```sh
+>>>
 0 |00⟩
 1 |01⟩
 2 |10⟩
@@ -345,14 +336,12 @@ which gives:
 as one would expect.
 
 We can specify the number of particles as well:
-```sh
+```pycon
 basis = ls.SpinlessFermionBasis(4, number_particles=3) #We create fermionic basis on 4 sites with only 3 fermions
 basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 
-```
-which gives:
-```sh
+>>>
 7 |0111⟩
 11 |1011⟩
 13 |1101⟩
@@ -364,14 +353,12 @@ We can see that the basis consists of states with three fermions.
 
 The last case includes fermions with spin. The binary representations of basis states can be read as a pair of (fermions with spin up on a lattice, fermions with spin down on a lattice).
 We can create a basis of spinful fermions as follows:
-```sh
+```pycon
 basis = ls.SpinfulFermionBasis(2) #We create basis of spinful fermions on 2 sites
 basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 
-```
-which gives:
-```sh
+>>>
 0 |00⟩|00⟩
 1 |00⟩|01⟩
 2 |00⟩|10⟩
@@ -393,25 +380,22 @@ We see that binary representation now means the second quantization of fermions
 
 As before, we can specify a sector of the total Hilbert space with a given number of fermions with spin down and spin up:
 
-```sh
+```pycon
 basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins up and down (N_up, N_down)=(2,1)
 basis.build()
 for i in range(basis.number_states):
     print(basis.states[i], basis.state_to_string(basis.states[i])) 
-```
-which gives:
-```sh
+>>>
 7 |01⟩|11⟩
 11 |10⟩|11⟩
 ```
-as expected.
 
 #### Basis from Expressions
 
 Before we created basises explicitly, however, there is a way to construct a basis directly from expressions.
 However, in this case, one can obtain only full Fock basis, without restriction on the number of particles.
 
-```
+```pycon
 expr=ls.Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁")
 basis=expr.full_basis()
 basis.build()
@@ -472,16 +456,16 @@ The symmetries are constructed with the help of expressions, and are represented
 `lattice_symmetries` uses sympy to represent permutations, therefore one can take a look at sympy documentation for more details.
 
 Let's take a look at a couple of examples:
-```
+```pycon
 p = ls.Permutation([1,2,3,0]) #This permutation shifts the indices, so that 0->1, 1->2, 2->3, 3->0
 >>> (0 1 2 3) #The same permutation written in the cycle representation
 ```
-```
+```pycon
 p = ls.Permutation([0,1,2,3]) #This is the identity permutation
 >>> (3) #No cycles, therefore identity. 
 #If the largest index doesn't move, it is shown in brackets, the format is used by sympy
 ```
-```
+```pycon
 p=ls.Permutation([1,0,3,2]) #Exchange indices 0<->1 and 2<->3
 >>> (0 1)(2 3) #Two cycles
 ```
@@ -495,18 +479,18 @@ Since we need the characters to construct symmetry adapted basises, there are tw
 
 - One can find all permutation symmetries of an expression:
 
-```
+```pycon
  #let's use a simple chain of 3 spins, which has the symmetry group D_3
 e=ls.Expr("σ^x_0 σ^x_1")
 expr=e.on(ig.Graph.Lattice(dim=[3], circular=True)) # The periodic chain with 3 sites
 sym=expr.permutation_group()
 >>>
 
-```
+```pycon
 This option can be used to study specific sectors and does not cover the whole Hilbert space.
 
 - Another option is to find the maximum abelian subgroup of the symmetry group:
-```
+```pycon
 ab_sym=expr.abelian_permutation_group()
 >>>
 
@@ -521,7 +505,7 @@ The basis is constructed with the help of one dimensional representations of the
 or in other words, the symmetry (permutation) group of the corresponding expression.
 
 The simplest example would be:
-```
+```pycon
 n=3
 translation = ls.Permutation([(1 + i) % n for i in range(n)]) #we consider one-dimensional translations as before
 b = ls.SpinBasis(number_spins=n, symmetries=[(translation, ls.Rational(k, n))])

From 30e43858ea0c4c755730feec283ef765749f08e0 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 18:11:34 +0200
Subject: [PATCH 45/67] grammar check

---
 python/tutorial.md | 134 ++++++++++++++++++++++-----------------------
 1 file changed, 66 insertions(+), 68 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index e64f018..0e80ca6 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -18,7 +18,7 @@
     - [Symmetry](#Symmetry)
 		- [Symmetry as Permutations](#Symmetries-as-permutations)
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
-		- [Symmetry adapted basis](#Symmetry-adapted-basis)
+		- [Symmetry-adapted basis](#Symmetry-adapted-basis)
 - [Examples](#Examples)
     - [Simple ED](#Simple-ED)
     - [More complicated ED](#More-complicated-ED)
@@ -26,7 +26,7 @@
 
 ## Installing
 
-The first step before we can apply `lattice_symmetries` is to install Nix. The installation process slightly depends on your system and can be found in the [official documentaion](https://nix.dev/install-nix#install-nix).
+The first step before we can apply `lattice_symmetries` is to install Nix. The installation process slightly depends on your system and can be found in the [official documentation](https://nix.dev/install-nix#install-nix).
 In the case of WSL (Windows Subsystem for Linux), it is preferred to install a single-user version:
 
 ```sh
@@ -66,8 +66,8 @@ Now you can build everything using Nix:
 nix develop .#python
 ```
 
-Voila, and you have the working package.
-If you open a new terminal, the last step should be repeated ,i.e., moving to  `lattice-symmetries/python` and typing `nix develop .#python`.
+Voila and you have the working package.
+If you open a new terminal, the last step should be repeated i.e., moving to  `lattice-symmetries/python` and typing `nix develop .#python`.
 
 ## Basic concepts and functions
 
@@ -78,35 +78,35 @@ such as time evolution or exact diagonalization of many-body Hamiltonians.
 
 The basic objects upon which the machinery is built are `Basis`, `Expression`, `Operator`, and `Symmetry`.
 The `Symmetries` is not a separate class, however, it lies in the core of `lattice_symmetries`, so we will consider it separately.
-- `Expression` object is a nice way to work with symbolic representations of operators. It is possible to sum expressions, and multiply them by numbers and each other.
-`Expression` allows not to think about an explicit matrix representation of an operator, and the user can work directly with analytical formulae. 
+- `Expression` object is a nice way to work with symbolic representations of operators. It is possible to sum expressions and multiply them by numbers and each other.
+`Expression` does not require the user to consider an explicit matrix representation of an operator and allows the user to work directly with analytical formulae. 
 - `Basis` object stores a basis of a many-body Hilbert space consisting of spins or fermions with appropriate symmetries.
 Each basis state is represented as a sequence of 0s and 1s (i.e. a sequence of bits), which can be interpreted as a binary number.
 - `Operator` object is an actual operator that can act on individual basis states and their linear combinations. 
-- `Symmetry` is a method to work with symmetry adapted basises.
-If an operator has symmetries, it is useful to work in symmetry-adapted basis, since it can drastically decrease the dimension of the Hilbert space.
+- `Symmetry` is a method to work with symmetry-adapted basis.
+If an operator has symmetries, it is useful to work with a symmetry-adapted basis since it can drastically decrease the dimension of the Hilbert space.
 
-Now we will take a look at basic functions and methods for these objects. 
+Now, we will take a look at basic functions and methods for these objects. 
 
 ### Expressions
 
-Expressions are an easy way to work with operators (for example, Hamiltonians) on a symbolic level using second quantization formalism.
+Expressions are an easy way to work with operators (such as Hamiltonians) on a symbolic level using the second quantization formalism.
 This means that you can use primitive symbols of well-known operators such as $\sigma^x$, $\sigma^y$, and $\sigma^z$ to build expressions for your Hamiltonian and observables.
 It is also possible to sum different expressions and multiply them to each other to compose more complicated expressions.
-Let's consider at first the simplest examples. 
+Let's consider the simplest examples first. 
 
 #### Primitive expressions
 
-At first we need to import `Expr` from lattice-symmetries:
+At first, we need to import `Expr` from lattice-symmetries:
 ```pycon
 from lattice_symmetries import Expr
 ```
-Now we will consider primitive symbols defined on site with number 0 as an example, but you can also construct them residing on other lattice sites:
+Now we will consider primitive symbols defined on a site with index 0 as an example, but you can also construct them residing on other lattice sites:
 
  - $\sigma^x$ and $S^x$:
    ```pycon
    >>> Expr("σˣ₀") == Expr("\\sigma^x_0") #It is possible to use different notations for Expr for primitives.
-   #Here we check that they agree. Index 0 means the index of the corresponding site.
+   #Here, we check that they agree. Index 0 means the index of the corresponding site.
    True
    >>> Expr("Sˣ₀") == 0.5 * Expr("σˣ₀") #We check that Sˣ₀ is the same as 0.5*σˣ₀.
    True
@@ -172,7 +172,7 @@ True
 >>> Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁") == Expr("0.5 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + 0.25 σᶻ₀ σᶻ₁")
 True
 
-#Above we checked that different definitions of this interaction agree. Let's take a look at the corresponding matrix:
+#Above, we checked that different definitions of this interaction agree. Let's take a look at the corresponding matrix:
 >>> Expr("Sx0 Sx1 + Sy0 Sy1 + Sz0 Sz1").to_dense()
 
 [[ 0.25  0.    0.    0.  ]
@@ -193,8 +193,8 @@ Under the hood, lattice-symmetries rewrites all the expressions into the canonic
 ```
 #### Complex expressions
 
-So far, we defined expressions only on a few number of sites, which can be indeed easily written explicitly. Here we will consider more complicated expressions, which require other techniques.
-One of the ways is to use the function `replace_indices`. For example, we can construct the sum of $\sigma^z$s defined on different sites: 
+So far, we have defined expressions only for a system with a few sites, which can be easily written explicitly. Here, we will consider more complicated expressions, which require other techniques.
+One of the ways to do this is by using the function `replace_indices`. For example, we can construct the sum of $\sigma^z$s defined on different sites: 
 
 ```pycon
 import operator #Import relevant methods
@@ -229,17 +229,17 @@ It is also possible to use `iGraph` to define an underlying (hyper)graph:
 ```pycon
 import igraph as ig
 
-e=Expr("σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀") #Here we have two-sites interaction, so we need an actual graph with edges consist of two vertices
+e=Expr("σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀") #Here, we have two-site interaction, so we need an actual graph with edges consisting of two vertices
 expr=e.on(ig.Graph.Lattice(dim=[2, 2])) # The graph is square 2x2
 >>>
 σ⁺₀ σ⁻₁ + σ⁺₀ σ⁻₂ + σ⁻₀ σ⁺₁ + σ⁻₀ σ⁺₂ + σ⁺₁ σ⁻₃ + σ⁻₁ σ⁺₃ + σ⁺₂ σ⁻₃ + σ⁻₂ σ⁺₃
 ```
 However, the method of defining a Hamiltonian on a graph works correctly only if the elementary expression (defined on one site) is homogeneous.
-It is also important that edges of a (hyper)graph should have the same number of vertices as the degree of the elementary expression.
+It is also important that the edges of a (hyper)graph should have the same number of vertices as the degree of the elementary expression.
 
 #### Properties
 
-One can make standard operations on expressions, such make an adjoint, as well
+One can make standard operations on expressions, such as making an adjoint, as well:
 
 ```pycon
 a = Expr("c†₀ c₁")
@@ -261,8 +261,9 @@ False
 
 ### Basis
 
-In this section we will consider generation of basises of various types.
-We will not consider symmetry adapted basises yet, we will do it in the section [Symmetry adapted basis](#Symmetry-adapted-basis).
+In this section, we will consider the generation of basises of various types.
+We will not consider symmetry-adapted bases yet; we will do so in the section [Symmetry-adapted basis](#Symmetry-adapted-basis).
+
 
 #### Spin basis
 Let's look at simple examples; at first we will not consider additional symmetries.
@@ -276,12 +277,15 @@ basis = ls.SpinBasis(3) #We create basis consisting of three $\frac{1}{2}$-spins
 basis.build() #build basis
 print(basis.index(basis.states)) #Print array of indices of basis states
 print(basis.number_states) #Print total number of states in the basis
-for i in range(basis.number_states): #Here we print all basis states as they stored in memory. Without symmetries, basis states equal their index
-    print(basis.states[i], basis.state_to_string(basis.states[i]))  
+
+def show_basis(basis): #Here we print all basis states as they are stored in memory 
+	for i in range(basis.number_states): 
+		print(basis.states[i], basis.state_to_string(basis.states[i])) 
+show_basis(basis)
 >>>
 [0 1 2 3 4 5 6 7]
 8
-0 |000⟩
+0 |000⟩ #Without symmetries, basis states are equal to their index
 1 |001⟩
 2 |010⟩
 3 |011⟩
@@ -292,26 +296,24 @@ for i in range(basis.number_states): #Here we print all basis states as they sto
 ```
 The basis states are equal to their indices and binary representations, as they should be.
 
-We can also consider only a part of the whole Hilbert space with a given number of spin ups (i.e. hamming weight of binary representations).
+We can also consider only a part of the whole Hilbert space with a given number of spins up (i.e. the hamming weight of binary representations).
 We can specify it as follows:
 ```pycon
 basis = ls.SpinBasis(3, hamming_weight=2) #We want the subspace with only 2 spins up 
 basis.build()
-for i in range(basis.number_states): #Print the states in the basis
-    print(basis.states[i], basis.state_to_string(basis.states[i]))  
+show_basis(basis)
 >>>
 3 |011⟩
 5 |101⟩
 6 |110⟩
 ```
-We see that basis states include only spins with 2 spins up. It is also interesting to note that in this case a basis state is not equal to its index.
+We see that basis states include only spins with 2 spins up. It is also interesting to note that, in this case, a basis state is not equal to its index.
 
 Sometimes, our system has spin inversion symmetry, and we can additionally shorten our basis. In this case, we can specify it as follows:
 ```pycon
 basis = ls.SpinBasis(4, spin_inversion=-1) #Spin inversion is present. It is not nessecary to specify hamming weight here, since it is fixed by spin inversion symmetry. 
 basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i]))  
+show_basis(basis) 
 >>>
 3 |011⟩
 5 |101⟩
@@ -319,14 +321,13 @@ for i in range(basis.number_states):
 ```
 
 #### Spinless Fermionic basis
-We can also consider the basis of fermions without spins. The basis states are stored as integers as for spin basis. However the binary representation has a second quantization interpretation.
+We can also consider the basis of fermions without spins. The basis states are stored as integers as for spin basis. However, the binary representation has a second quantization interpretation.
 Each basis state is given by the sequence of 0s and 1s, where 1 means a fermion on the corresponding site, and 0 means that the site is vacant.
 Let's consider the simplest example of fermions on two sites:
 ```pycon
-basis = ls.SpinlessFermionBasis(2) #We create fermionic basis on 2 sites
+basis = ls.SpinlessFermionBasis(2) #We create a fermionic basis on 2 sites
 basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+show_basis(basis)
 >>>
 0 |00⟩
 1 |01⟩
@@ -337,10 +338,9 @@ as one would expect.
 
 We can specify the number of particles as well:
 ```pycon
-basis = ls.SpinlessFermionBasis(4, number_particles=3) #We create fermionic basis on 4 sites with only 3 fermions
+basis = ls.SpinlessFermionBasis(4, number_particles=3) #We create a fermionic basis on 4 sites with only 3 fermions
 basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+show_basis(basis)
 >>>
 7 |0111⟩
 11 |1011⟩
@@ -354,10 +354,9 @@ We can see that the basis consists of states with three fermions.
 The last case includes fermions with spin. The binary representations of basis states can be read as a pair of (fermions with spin up on a lattice, fermions with spin down on a lattice).
 We can create a basis of spinful fermions as follows:
 ```pycon
-basis = ls.SpinfulFermionBasis(2) #We create basis of spinful fermions on 2 sites
+basis = ls.SpinfulFermionBasis(2) #We create a basis of spinful fermions on 2 sites
 basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+show_basis(basis)
 >>>
 0 |00⟩|00⟩
 1 |00⟩|01⟩
@@ -383,8 +382,7 @@ As before, we can specify a sector of the total Hilbert space with a given numbe
 ```pycon
 basis = ls.SpinfulFermionBasis(2, number_particles=(2,1)) #We specify the numbers of fermions with spins up and down (N_up, N_down)=(2,1)
 basis.build()
-for i in range(basis.number_states):
-    print(basis.states[i], basis.state_to_string(basis.states[i])) 
+show_basis(basis)
 >>>
 7 |01⟩|11⟩
 11 |10⟩|11⟩
@@ -392,8 +390,8 @@ for i in range(basis.number_states):
 
 #### Basis from Expressions
 
-Before we created basises explicitly, however, there is a way to construct a basis directly from expressions.
-However, in this case, one can obtain only full Fock basis, without restriction on the number of particles.
+Before we created a basis explicitly, however, there was a way to construct a basis directly from expressions.
+However, in this case, one can obtain only a full Fock basis without restriction on the number of particles.
 
 ```pycon
 expr=ls.Expr("Sˣ₀ Sˣ₁ + Sʸ₀ Sʸ₁ + Sᶻ₀ Sᶻ₁")
@@ -405,14 +403,14 @@ basis.build()
 ### Operators
 
 Based on an expression and a basis, we can build the corresponding operator acting on the chosen basis. Let's at first consider the Hubbard model on two sites without spins.
-One of the ways to contruct an operator is to create an expression and then create the operator:
+One of the ways to construct an operator is to create an expression and then create the operator:
 
 ```pycon
 import lattice_symmetries as ls
 
 op_expr = ls.Expr("-c†₀ c₁-c†₁ c₀+2.0 n₀ n₀+2.0 n₁ n₁")
 basis=ls.SpinlessFermionBasis(2)
-basis.build() #It is nessecary to build the basis before creating an operator
+basis.build() #It is necessary  to build the basis before creating an operator
 opr=ls.Operator(op_expr, basis) #Here we create the operator
 ```
 Another way is to create the expression for each small expression and then add them to each other:
@@ -427,7 +425,7 @@ opr.shape() #We can check the shape of the operator
 (4,4)
 ```
 
-We can make a multiplication of the operator on a vector using standart python notation.
+We can multiply the operator on a vector using standard Python notation.
 
 ```pycon
 vec1=np.array([1.0,0,0,0]) #One should work with floats
@@ -450,8 +448,8 @@ opr@vec3
 
 #### Symmetry as Permutations
 
-The full power of `lattice_symmetries` manifests if one use symmetries when constructing 
-symmetry adapted basis and linear operators acting on the corresponding Hilbert space. 
+The full power of `lattice_symmetries` manifests if one uses symmetries when constructing 
+symmetry-adapted basis and linear operators acting on the corresponding Hilbert space. 
 The symmetries are constructed with the help of expressions, and are represented as a permutation group of indices;
 `lattice_symmetries` uses sympy to represent permutations, therefore one can take a look at sympy documentation for more details.
 
@@ -463,7 +461,7 @@ p = ls.Permutation([1,2,3,0]) #This permutation shifts the indices, so that 0->1
 ```pycon
 p = ls.Permutation([0,1,2,3]) #This is the identity permutation
 >>> (3) #No cycles, therefore identity. 
-#If the largest index doesn't move, it is shown in brackets, the format is used by sympy
+#If the largest index doesn't move, it is shown in brackets; the format is used by sympy
 ```
 ```pycon
 p=ls.Permutation([1,0,3,2]) #Exchange indices 0<->1 and 2<->3
@@ -472,21 +470,21 @@ p=ls.Permutation([1,0,3,2]) #Exchange indices 0<->1 and 2<->3
 
 #### Symmetry from Expressions
 
-In the previous section, we constructed the symmetries by hand, however, this can be done only for relatively easy and small systems.
-In most of the cases, the more powerful and convenient way is to ask `lattice_symmetries` to calculate the simmetries of a given expression. 
+In the previous section, we constructed the symmetries by hand, however, this can only be done for relatively small and simple systems.
+In most cases, the more powerful and convenient way is to ask `lattice_symmetries` to calculate the symmetries of a given expression. 
 
-Since we need the characters to construct symmetry adapted basises, there are two possibilities.
+Since we need the characters to construct symmetry-adapted basis, there are two possibilities.
 
 - One can find all permutation symmetries of an expression:
 
 ```pycon
- #let's use a simple chain of 3 spins, which has the symmetry group D_3
+#Let's use a simple chain of 3 spins, which has the symmetry group D_3
 e=ls.Expr("σ^x_0 σ^x_1")
 expr=e.on(ig.Graph.Lattice(dim=[3], circular=True)) # The periodic chain with 3 sites
 sym=expr.permutation_group()
 >>>
 
-```pycon
+```
 This option can be used to study specific sectors and does not cover the whole Hilbert space.
 
 - Another option is to find the maximum abelian subgroup of the symmetry group:
@@ -496,12 +494,12 @@ ab_sym=expr.abelian_permutation_group()
 
 ```
 
-In this case the sectors cover the whole Hilbert space.
+In this case, the sectors cover the whole Hilbert space.
 
-#### Symmetry adapted basis
+#### Symmetry-adapted basis
 
-In order to make calculations with the help of symmetries, we need to construct a symmetry adapted basis.
-The basis is constructed with the help of one dimensional representations of the symmetry group of the Hamiltonian,
+To make calculations with the help of symmetries, we need to construct a symmetry-adapted basis.
+The basis is constructed with the help of one-dimensional representations of the Hamiltonian's symetry group,
 or in other words, the symmetry (permutation) group of the corresponding expression.
 
 The simplest example would be:
@@ -517,11 +515,11 @@ An `Operator` for symmetry adapted basis can be build in the same way as without
 
 ## Examples
 
-Here we will take a look at different examples of `lattice_symmetries` applications.
+Here, we will take a look at different examples of `lattice_symmetries` applications.
 
 ### Simple ED
 
-We will strat with the simplest example of exact diagonalization. We will consider Heisenberg chain on 10 sites and diagonalization with the help of symmetries.
+We will start with the simplest example of exact diagonalization. We will consider Heisenberg chain on 10 sites and diagonalization with the help of symmetries.
 For that, we will combine methods described in the previous sections.
 
 ```pycon
@@ -544,11 +542,11 @@ basis = ls.SpinBasis(
 	# NOTE: we don't actually need to specify hamming_weight when spin_inversion
 	# is set. The library will guess that hamming_weight = number_spins // 2.
 	spin_inversion=-1,
-	symmetries=[T, P], #Here we use the characters of the whole symmetry group (it's worth noting that, in general, translations and parity don't commute)
+	symmetries=[T, P], #Here, we use the characters of the whole symmetry group (it's worth noting that, in general, translations and parity don't commute)
 )
-basis.build()  # Build the list of representatives, we need it since we're doing ED
+basis.build()  # Build the list of representatives. We need it since we're doing ED
 
-# Constructing the expression of the hamiltonian
+# Constructing the expression of the Hamiltonian
 edges = [(i, (i + 1) % number_spins) for i in range(number_spins)]
 expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁", sites=edges)
 print("Expression:", expr)
@@ -568,7 +566,7 @@ eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA
 print("Ground state energy is {}".format(eigenvalues[0]))
 assert np.isclose(eigenvalues[0], -18.06178542)
 
-#Here we used a symmetry adapted basis, and therefore we were restricted to a specific sector of the total Hilbert space.
+#Above, we used the symmetry-adapted basis, and therefore we were restricted to the specific sector of the total Hilbert space.
 #Since the considered system is relatively small, we can make the diagonalization in the full basis and check that we chose the right sector.
 
 new_basis = ls.SpinBasis(
@@ -596,8 +594,8 @@ import scipy
 
 ### Time evolution
 
-Another example of capabilities of `lattice_symmetries` is time evolution.
-In order to apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
+Another example of the capabilities of `lattice_symmetries` is time evolution.
+To apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
 e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum\limits^{\infty} 2(-i)^n J_n(at)T_n(H')]

From bdde5cdc65e0946266931c6ef6ea9921c4521444 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 11 Apr 2024 19:03:22 +0200
Subject: [PATCH 46/67] Update tutorial.md

---
 python/tutorial.md | 59 ++++++++++++++++++++++++++++++++--------------
 1 file changed, 41 insertions(+), 18 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 0e80ca6..ead3969 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -86,7 +86,7 @@ Each basis state is represented as a sequence of 0s and 1s (i.e. a sequence of b
 - `Symmetry` is a method to work with symmetry-adapted basis.
 If an operator has symmetries, it is useful to work with a symmetry-adapted basis since it can drastically decrease the dimension of the Hilbert space.
 
-Now, we will take a look at basic functions and methods for these objects. 
+Now, we will take a look at basic functions and methods for these entities. 
 
 ### Expressions
 
@@ -475,43 +475,66 @@ In most cases, the more powerful and convenient way is to ask `lattice_symmetrie
 
 Since we need the characters to construct symmetry-adapted basis, there are two possibilities.
 
-- One can find all permutation symmetries of an expression:
+- One can find all permutation symmetries of an expression.
+This option can be used to study specific sectors and does not cover the whole Hilbert space.
+
+As a simple test case, we can consider a chain of 3 spins with the symmetry group $D_3=S_3$. This group is already non-abelian.
 
 ```pycon
-#Let's use a simple chain of 3 spins, which has the symmetry group D_3
 e=ls.Expr("σ^x_0 σ^x_1")
 expr=e.on(ig.Graph.Lattice(dim=[3], circular=True)) # The periodic chain with 3 sites
 sym=expr.permutation_group()
 >>>
-
-```
-This option can be used to study specific sectors and does not cover the whole Hilbert space.
-
-- Another option is to find the maximum abelian subgroup of the symmetry group:
+PermutationGroup([
+    (1 2),
+    (2)(0 1),
+    (0 1 2),
+    (0 2 1),
+    (0 2)]) 
+#The identity element is always present, so it is not shown.
+#The group has 6 elements, as it should
+```
+
+- Another option is to find the maximum abelian subgroup of the symmetry group. In this case, the sectors cover the whole Hilbert space.
+In the case of the same 3-site chain, we have:
 ```pycon
 ab_sym=expr.abelian_permutation_group()
 >>>
-
+PermutationGroup([
+    (0 1 2), #translation by 1
+    (0 2 1)]) #translation by 2
+#The maximal abelian sugroup consists of 2 translations and identity
 ```
 
-In this case, the sectors cover the whole Hilbert space.
-
 #### Symmetry-adapted basis
 
 To make calculations with the help of symmetries, we need to construct a symmetry-adapted basis.
-The basis is constructed with the help of one-dimensional representations of the Hamiltonian's symetry group,
-or in other words, the symmetry (permutation) group of the corresponding expression.
+The basis is constructed with the help of one-dimensional representations of the operator symmetry group,
+or, in other words, the symmetry (permutation) group of the corresponding expression.
+
+We describe the characters of the symmetry group as a list of tuples:
+```pycon
+[(permutation, rational number)]
+```
+where a permutation is a symmetry of the expression, and a rational number is its phase, defining the Hilbert space sector.
+One can think about this list as a list of symmetry generators with corresponding characters.
+If one considers only abelian symmetries, then the generator characters can be any rational number of the form $k/n$, 
+where $n$ is the order of the generator and $k$ is a natural number.
+However, if the symmetry group is non-abelian, the phases should be chosen in a proper way,
+so that characters of group elements organize a one-dimension representation of the symmetry group.
 
 The simplest example would be:
 ```pycon
-n=3
-translation = ls.Permutation([(1 + i) % n for i in range(n)]) #we consider one-dimensional translations as before
-b = ls.SpinBasis(number_spins=n, symmetries=[(translation, ls.Rational(k, n))])
+p = ls.Permutation([1,2,0])
+b = ls.SpinBasis(3, symmetries=[(p, ls.Rational(0, 1))]) #We specify the phase to be zero. It is the trivial character
 b.build()
 ```
 
-An `Operator` for symmetry adapted basis can be build in the same way as without symmetries.
+An `Operator` for symmetry-adapted basis is built in the same way as without symmetries:
 
+```pycon
+opr=ls.Operator(expr,basis)
+```
 
 ## Examples
 
@@ -608,4 +631,4 @@ $$
 H'=\frac{H-b}{a}
 $$
 
-where $a=1/(2-\epsilon)$, $b=(E_g+E_s)/2$, and $\epsilon$ is a safety factor to make the spectrum of $H'$ liying within $[-1,1]$.
\ No newline at end of file
+where $a=1/(2-\epsilon)$, $b=(E_g+E_s)/2$, and $\epsilon$ is a safety factor to make the spectrum of $H'$ be within $[-1,1]$.
\ No newline at end of file

From 985fb11e2491e6dc0bafe59c6195644fab0feafd Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Fri, 12 Apr 2024 11:45:17 +0200
Subject: [PATCH 47/67] Update tutorial.md

---
 python/tutorial.md | 81 +++++++++++++++++++++++++++++++++-------------
 1 file changed, 59 insertions(+), 22 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index ead3969..e2e6079 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -444,6 +444,9 @@ opr@vec3
 [ 0.  4. -2.  0.]
 ```
 
+The operators for the symmetry-adapted basis are created and manipulated in the same way.
+We discuss this in more detail in the section [Symmetry-adapted basis](#Symmetry-adapted-basis).
+
 ### Symmetry
 
 #### Symmetry as Permutations
@@ -451,7 +454,7 @@ opr@vec3
 The full power of `lattice_symmetries` manifests if one uses symmetries when constructing 
 symmetry-adapted basis and linear operators acting on the corresponding Hilbert space. 
 The symmetries are constructed with the help of expressions, and are represented as a permutation group of indices;
-`lattice_symmetries` uses sympy to represent permutations, therefore one can take a look at sympy documentation for more details.
+`lattice_symmetries` uses sympy to represent permutations, therefore one can take a look at [sympy documentation](https://docs.sympy.org/latest/modules/combinatorics/permutations.html) for more details.
 
 Let's take a look at a couple of examples:
 ```pycon
@@ -506,6 +509,11 @@ PermutationGroup([
 #The maximal abelian sugroup consists of 2 translations and identity
 ```
 
+These functions return a sympy permutation group. This allows us to apply all the existing sympy functions and analyze the symmetries of expressions.
+For more details, one can look [at sympy documentation](https://docs.sympy.org/latest/modules/combinatorics/perm_groups.html).
+
+Here, we will consider a few methods that can be particularly useful.
+
 #### Symmetry-adapted basis
 
 To make calculations with the help of symmetries, we need to construct a symmetry-adapted basis.
@@ -517,8 +525,16 @@ We describe the characters of the symmetry group as a list of tuples:
 [(permutation, rational number)]
 ```
 where a permutation is a symmetry of the expression, and a rational number is its phase, defining the Hilbert space sector.
-One can think about this list as a list of symmetry generators with corresponding characters.
-If one considers only abelian symmetries, then the generator characters can be any rational number of the form $k/n$, 
+To be more precise, the phase is described as follows:
+
+$$
+\hat P |\psi\rangle=e^{2\pi i \phi}|\psi\rangle
+$$
+
+where $\hat P$ is the operator representing the permutation, $\psi$ is its eigenvector and $\phi$ is the rational number.
+
+One can think about this list of tuples a list of symmetry generators with corresponding phases.
+If one considers only abelian symmetries, then the generator phases can be any rational number of the form $k/n$, 
 where $n$ is the order of the generator and $k$ is a natural number.
 However, if the symmetry group is non-abelian, the phases should be chosen in a proper way,
 so that characters of group elements organize a one-dimension representation of the symmetry group.
@@ -526,10 +542,21 @@ so that characters of group elements organize a one-dimension representation of
 The simplest example would be:
 ```pycon
 p = ls.Permutation([1,2,0])
-b = ls.SpinBasis(3, symmetries=[(p, ls.Rational(0, 1))]) #We specify the phase to be zero. It is the trivial character
-b.build()
+b0 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(0, 1))]) #We specify the phase to be zero. It defines the trivial character
+b0.build()
+```
+
+The order of the permutation is 3, wo we can have 2 other possible sectors:
+```pycon
+b1 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(1, 3))]) #The phase is 1/3, the total phase is 2/3π
+b1.build()
+
+b2 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(2, 3))]) #The phase is 2/3, the total phase is 4/3π
+b2.build()
 ```
+When making calculations, one needs to choose which sectors are suitable for the given task.
 
+When creating the `Operator`, one needs to be sure that the basis symmetries are consistent with the symmetries of the correspondning expression.
 An `Operator` for symmetry-adapted basis is built in the same way as without symmetries:
 
 ```pycon
@@ -545,20 +572,38 @@ Here, we will take a look at different examples of `lattice_symmetries` applicat
 We will start with the simplest example of exact diagonalization. We will consider Heisenberg chain on 10 sites and diagonalization with the help of symmetries.
 For that, we will combine methods described in the previous sections.
 
+First, we generate the expression for the Hamiltonian.
+
 ```pycon
 import lattice_symmetries as ls
 import numpy as np
 import scipy
 
+# Constructing the expression of the Hamiltonian
+edges = [(i, (i + 1) % number_spins) for i in range(number_spins)]
+expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁", sites=edges)
+print("Expression:", expr)
+
+# Alternatively, we can create expr in an algebraic way:
+two_site_expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁")
+many_exprs = [two_site_expr.replace_indices({0: i, 1: j}) for (i, j) in edges]
+expr2 = reduce(operator.add, many_exprs)
+assert expr == expr2 #we check that the expressions are equal
 number_spins = 10  # System size
+```
 
+Since the lattice is simple, we can construct the symmetries by hand and use one-dimensional representation of the whole symmetry group.
+```pycon
 # Constructing symmetries
 sites = np.arange(number_spins)
 # Momentum in x direction with phase φ=1/2 (i.e., with the eigenvalue λ=exp(-2πiφ)=-1)
 T = (ls.Permutation((sites + 1) % number_spins), ls.Rational(1, 2))
 # Parity with eigenvalue λ=-1
 P = (ls.Permutation(sites[::-1]), ls.Rational(1, 2))
+```
 
+Now, we construct the spin basis with 0 average magnetization.
+```pycon
 # Constructing the basis
 basis = ls.SpinBasis(
 	number_spins=number_spins,
@@ -568,18 +613,9 @@ basis = ls.SpinBasis(
 	symmetries=[T, P], #Here, we use the characters of the whole symmetry group (it's worth noting that, in general, translations and parity don't commute)
 )
 basis.build()  # Build the list of representatives. We need it since we're doing ED
-
-# Constructing the expression of the Hamiltonian
-edges = [(i, (i + 1) % number_spins) for i in range(number_spins)]
-expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁", sites=edges)
-print("Expression:", expr)
-
-# Alternatively, we can create expr in an algebraic way:
-two_site_expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁")
-many_exprs = [two_site_expr.replace_indices({0: i, 1: j}) for (i, j) in edges]
-expr2 = reduce(operator.add, many_exprs)
-assert expr == expr2 #we check that the expressions are equal
-
+```
+ The last step is to construct the actual Hamiltonian matrix and diagonalize it.
+```pycon
 # Construct the Hamiltonian
 hamiltonian = ls.Operator(expr, basis)
 
@@ -588,10 +624,11 @@ hamiltonian = ls.Operator(expr, basis)
 eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA")
 print("Ground state energy is {}".format(eigenvalues[0]))
 assert np.isclose(eigenvalues[0], -18.06178542)
+```
+Above, we used the symmetry-adapted basis, and therefore we were restricted to the specific sector of the total Hilbert space.
+Since the considered system is relatively small, we can make the diagonalization in the full basis and check that we chose the right sector.
 
-#Above, we used the symmetry-adapted basis, and therefore we were restricted to the specific sector of the total Hilbert space.
-#Since the considered system is relatively small, we can make the diagonalization in the full basis and check that we chose the right sector.
-
+```pycon
 new_basis = ls.SpinBasis(
 	number_spins=number_spins,
 	spin_inversion=-1,
@@ -601,9 +638,9 @@ new_hamiltonian = ls.Operator(expr, new_basis)
 eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA")
 print("Ground state energy is {}".format(eigenvalues[0]))
 assert np.isclose(eigenvalues[0], -18.06178542)
-
-#Indeed, the outcome is the same within machine precision. 
 ```
+Indeed, the outcome is the same within machine precision. 
+
 
 ### More complicated ED
 

From b020112047eb7dab91143c60db014d3369b6f3d5 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Fri, 12 Apr 2024 12:10:05 +0200
Subject: [PATCH 48/67] Update tutorial.md

---
 python/tutorial.md | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index e2e6079..ce2371d 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -16,7 +16,7 @@
 		- [Basis from Expressions](#Basis-from-expressions)
     - [Operators](#Operators)
     - [Symmetry](#Symmetry)
-		- [Symmetry as Permutations](#Symmetries-as-permutations)
+		- [Symmetry as Permutations](#Symmetry-as-permutations)
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry-adapted basis](#Symmetry-adapted-basis)
 - [Examples](#Examples)
@@ -531,9 +531,9 @@ $$
 \hat P |\psi\rangle=e^{2\pi i \phi}|\psi\rangle
 $$
 
-where $\hat P$ is the operator representing the permutation, $\psi$ is its eigenvector and $\phi$ is the rational number.
+where $\hat P$ is the operator representing the permutation, $\psi$ is a vector from the chosen sector and $\phi$ is the rational number.
 
-One can think about this list of tuples a list of symmetry generators with corresponding phases.
+One can think about this list of tuples as a list of symmetry generators with corresponding phases.
 If one considers only abelian symmetries, then the generator phases can be any rational number of the form $k/n$, 
 where $n$ is the order of the generator and $k$ is a natural number.
 However, if the symmetry group is non-abelian, the phases should be chosen in a proper way,
@@ -556,7 +556,7 @@ b2.build()
 ```
 When making calculations, one needs to choose which sectors are suitable for the given task.
 
-When creating the `Operator`, one needs to be sure that the basis symmetries are consistent with the symmetries of the correspondning expression.
+When creating an `Operator`, one needs to be sure that the basis symmetries are consistent with the symmetries of the correspondning expression.
 An `Operator` for symmetry-adapted basis is built in the same way as without symmetries:
 
 ```pycon
@@ -593,6 +593,7 @@ number_spins = 10  # System size
 ```
 
 Since the lattice is simple, we can construct the symmetries by hand and use one-dimensional representation of the whole symmetry group.
+It is worth noting that it is only possible, because we choose the proper sector for parity and translations permutations, which don't commute.
 ```pycon
 # Constructing symmetries
 sites = np.arange(number_spins)
@@ -602,7 +603,7 @@ T = (ls.Permutation((sites + 1) % number_spins), ls.Rational(1, 2))
 P = (ls.Permutation(sites[::-1]), ls.Rational(1, 2))
 ```
 
-Now, we construct the spin basis with 0 average magnetization.
+Now, we construct the spin basis with zero magnetization.
 ```pycon
 # Constructing the basis
 basis = ls.SpinBasis(
@@ -610,7 +611,7 @@ basis = ls.SpinBasis(
 	# NOTE: we don't actually need to specify hamming_weight when spin_inversion
 	# is set. The library will guess that hamming_weight = number_spins // 2.
 	spin_inversion=-1,
-	symmetries=[T, P], #Here, we use the characters of the whole symmetry group (it's worth noting that, in general, translations and parity don't commute)
+	symmetries=[T, P], #Here, we use the characters of the whole symmetry group 	
 )
 basis.build()  # Build the list of representatives. We need it since we're doing ED
 ```

From f8367e491d8cb8d76faac711e878d87a98f45456 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Fri, 12 Apr 2024 16:29:15 +0200
Subject: [PATCH 49/67] almost finished draft

---
 python/tutorial.md | 108 +++++++++++++++++++++++++++++++++++----------
 1 file changed, 85 insertions(+), 23 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index ce2371d..dd56d6c 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -20,8 +20,8 @@
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry-adapted basis](#Symmetry-adapted-basis)
 - [Examples](#Examples)
-    - [Simple ED](#Simple-ED)
-    - [More complicated ED](#More-complicated-ED)
+    - [ED, 1D chain](#Ed,-1d-chain)
+    - [ED, Star graph](#Ed,-star-graph)
     - [Time evolution](#Time-evolution)
 
 ## Installing
@@ -454,7 +454,7 @@ We discuss this in more detail in the section [Symmetry-adapted basis](#Symmetry
 The full power of `lattice_symmetries` manifests if one uses symmetries when constructing 
 symmetry-adapted basis and linear operators acting on the corresponding Hilbert space. 
 The symmetries are constructed with the help of expressions, and are represented as a permutation group of indices;
-`lattice_symmetries` uses sympy to represent permutations, therefore one can take a look at [sympy documentation](https://docs.sympy.org/latest/modules/combinatorics/permutations.html) for more details.
+`lattice_symmetries` uses sympy to represent permutations, therefore one can take a look at [`sympy` documentation](https://docs.sympy.org/latest/modules/combinatorics/permutations.html) for more details.
 
 Let's take a look at a couple of examples:
 ```pycon
@@ -471,6 +471,19 @@ p=ls.Permutation([1,0,3,2]) #Exchange indices 0<->1 and 2<->3
 >>> (0 1)(2 3) #Two cycles
 ```
 
+We can use some of the `sympy` functions to analyze permutations.
+For example, if one wants to find the order of a permutation:
+```pycon
+p = ls.Permutation([1,2,3,0])
+p.order()
+>>>
+4
+>>>
+p**(p.order()) #check that p^4==1
+>>>
+Permutation(3)
+```
+
 #### Symmetry from Expressions
 
 In the previous section, we constructed the symmetries by hand, however, this can only be done for relatively small and simple systems.
@@ -484,7 +497,9 @@ This option can be used to study specific sectors and does not cover the whole H
 As a simple test case, we can consider a chain of 3 spins with the symmetry group $D_3=S_3$. This group is already non-abelian.
 
 ```pycon
-e=ls.Expr("σ^x_0 σ^x_1")
+import igraph as ig
+
+e=ls.Expr("σ^x_0 σ^x_1") #A simple expression defined on an edge
 expr=e.on(ig.Graph.Lattice(dim=[3], circular=True)) # The periodic chain with 3 sites
 sym=expr.permutation_group()
 >>>
@@ -498,8 +513,10 @@ PermutationGroup([
 #The group has 6 elements, as it should
 ```
 
-- Another option is to find the maximum abelian subgroup of the symmetry group. In this case, the sectors cover the whole Hilbert space.
-In the case of the same 3-site chain, we have:
+- Another option is to find the maximum abelian subgroup of the symmetry group (usually it coincides with the translation subgroup).
+In this case, the sectors cover the whole Hilbert space.
+
+For the same 3-site chain, we have:
 ```pycon
 ab_sym=expr.abelian_permutation_group()
 >>>
@@ -509,10 +526,8 @@ PermutationGroup([
 #The maximal abelian sugroup consists of 2 translations and identity
 ```
 
-These functions return a sympy permutation group. This allows us to apply all the existing sympy functions and analyze the symmetries of expressions.
-For more details, one can look [at sympy documentation](https://docs.sympy.org/latest/modules/combinatorics/perm_groups.html).
-
-Here, we will consider a few methods that can be particularly useful.
+These functions return a sympy permutation group. This allows us to apply all the existing `sympy` functions and analyze the symmetries of expressions.
+For more details, one can look [at `sympy` documentation](https://docs.sympy.org/latest/modules/combinatorics/perm_groups.html).
 
 #### Symmetry-adapted basis
 
@@ -534,12 +549,13 @@ $$
 where $\hat P$ is the operator representing the permutation, $\psi$ is a vector from the chosen sector and $\phi$ is the rational number.
 
 One can think about this list of tuples as a list of symmetry generators with corresponding phases.
-If one considers only abelian symmetries, then the generator phases can be any rational number of the form $k/n$, 
-where $n$ is the order of the generator and $k$ is a natural number.
+However, the generators are not necessarily independent, and one can pass the whole symmetry group with the appropriate phases.
+If one considers only abelian symmetries and independent generators, then the generator phases can be any rational number of the form $k/n$, 
+where $n$ is the order of the generator, and $k$ is a natural number.
 However, if the symmetry group is non-abelian, the phases should be chosen in a proper way,
 so that characters of group elements organize a one-dimension representation of the symmetry group.
 
-The simplest example would be:
+Now, we are ready to build a symmetry-adapted basis. Th simplest example would be:
 ```pycon
 p = ls.Permutation([1,2,0])
 b0 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(0, 1))]) #We specify the phase to be zero. It defines the trivial character
@@ -556,18 +572,34 @@ b2.build()
 ```
 When making calculations, one needs to choose which sectors are suitable for the given task.
 
-When creating an `Operator`, one needs to be sure that the basis symmetries are consistent with the symmetries of the correspondning expression.
-An `Operator` for symmetry-adapted basis is built in the same way as without symmetries:
+Let's consider another example, where we will work with the symmetries of expressions:
+```pycon
+e=ls.Expr("σ^x_0 σ^x_1")
+expr=e.on(ig.Graph.Full(5)) # We define our lattice to be a complete graph with 5 vertices
+sym=expr.permutation_group()
+rep=[(p,ls.Rational(0,1)) for p in sym] #The trivial one-dimensional representation 
+#Notice that we can use the whole symmetry group for constructing 1D-representation
+basis = ls.SpinBasis(5, symmetries=rep) 
+basis.build()
+```
+
+The last step is to make calculations in symmetry adapted basis. For that we need to create an `Operator` object.
+When making an operator, one needs to be sure that the basis symmetries are consistent with the symmetries of the correspondning expression.
+An operator for symmetry-adapted basis is built in the same way as without symmetries. Schematically, the code would look as follows:
 
 ```pycon
-opr=ls.Operator(expr,basis)
+Make an expression
+Calculate symmetries
+Build symmetry-adapted basis
+
+opr=ls.Operator(expression,basis) # Make the operator in symmetry-adapted basis
 ```
 
 ## Examples
 
 Here, we will take a look at different examples of `lattice_symmetries` applications.
 
-### Simple ED
+### ED, 1D chain
 
 We will start with the simplest example of exact diagonalization. We will consider Heisenberg chain on 10 sites and diagonalization with the help of symmetries.
 For that, we will combine methods described in the previous sections.
@@ -579,6 +611,8 @@ import lattice_symmetries as ls
 import numpy as np
 import scipy
 
+number_spins = 10  # System size
+
 # Constructing the expression of the Hamiltonian
 edges = [(i, (i + 1) % number_spins) for i in range(number_spins)]
 expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁", sites=edges)
@@ -589,7 +623,6 @@ two_site_expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σ
 many_exprs = [two_site_expr.replace_indices({0: i, 1: j}) for (i, j) in edges]
 expr2 = reduce(operator.add, many_exprs)
 assert expr == expr2 #we check that the expressions are equal
-number_spins = 10  # System size
 ```
 
 Since the lattice is simple, we can construct the symmetries by hand and use one-dimensional representation of the whole symmetry group.
@@ -643,15 +676,42 @@ assert np.isclose(eigenvalues[0], -18.06178542)
 Indeed, the outcome is the same within machine precision. 
 
 
-### More complicated ED
-
-Now let's consider a more complicated example of ED.
+### ED, Star graph
 
+In the previous example, we considered a 1D chain, where we knew the graph's symmetries and could code them by hand.
+Sometimes, that task requires much time or is even practically impossible.
+Nevertheless, symmetries can still help, especially when the required sector lies in the one-dimensional representation of the whole symmetry group.
+As an example, we will consider the Heisenberg model on a star graph, where every vertex connects only to the center:
+ 
 ```pycon
 import lattice_symmetries as ls
 import numpy as np
 import scipy
+import igraph as ig
+
+number_spins=7
+expr = ls.Expr("2 (σ⁺₀ σ⁻₁ + σ⁺₁ σ⁻₀) + σᶻ₀ σᶻ₁").on(ig.Graph.Star(number_spins))
+sym=expr.permutation_group()
+rep=[(p,ls.Rational(0,1)) for p in sym] #The trivial one-dimensional representation 
+basis = ls.SpinBasis(number_spins, symmetries=rep) #We don't choose specific magnetization
+basis.build()
+hamiltonian = ls.Operator(expr, basis)
+eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA")
+print("Ground state energy is {}".format(eigenvalues[0]))
+assert np.isclose(eigenvalues[0], -8)
 ```
+Since the system is small, we can check that we chose the right sector:
+
+```pycon
+basis = ls.SpinBasis(number_spins) #No symmetries
+basis.build()
+hamiltonian = ls.Operator(expr, basis)
+eigenvalues, eigenstates = scipy.sparse.linalg.eigsh(hamiltonian, k=1, which="SA")
+print("Ground state energy is {}".format(eigenvalues[0]))
+assert np.isclose(eigenvalues[0], -8)
+```
+
+Of course, the considered example is somewhat artificial, nevertheless, it nicely demonstrates one of the possible applications of `lattice_symmetries`.
 
 ### Time evolution
 
@@ -659,7 +719,7 @@ Another example of the capabilities of `lattice_symmetries` is time evolution.
 To apply time evolution, we will use the Chebyshev expansion of the matrix exponent:
 
 $$
-e^{iHt}=e^{-i(E^{*}_g+aW')t}[J_0(at)+\sum\limits^{\infty} 2(-i)^n J_n(at)T_n(H')]
+e^{iHt}=e^{-i\frac{bt}{a}}[J_0(at)+\sum\limits^{\infty} 2(-i)^n J_n(at)T_n(H')]
 $$
 
 where we rescale the original Hamiltonian $H$ with bandwidth $[E_g, E_s]$ to the Hamiltonian $H'$ with bandwidth $[-1+\epsilon,1-\epsilon]$, so that the series converges.
@@ -669,4 +729,6 @@ $$
 H'=\frac{H-b}{a}
 $$
 
-where $a=1/(2-\epsilon)$, $b=(E_g+E_s)/2$, and $\epsilon$ is a safety factor to make the spectrum of $H'$ be within $[-1,1]$.
\ No newline at end of file
+where $a=1/(2-\epsilon)$, $b=(E_g+E_s)/2$, and $\epsilon$ is a safety factor to make the spectrum of $H'$ be within $[-1,1]$.
+
+We will consider the basic version of the algorithm, where instead of usual `numpy` arrays we will use `lattice_symmetries` operators.
\ No newline at end of file

From b5b49123d679d84e875bea534930473775136cae Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Fri, 12 Apr 2024 16:37:49 +0200
Subject: [PATCH 50/67] Update tutorial.md

---
 python/tutorial.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index dd56d6c..1572dc2 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -20,8 +20,8 @@
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry-adapted basis](#Symmetry-adapted-basis)
 - [Examples](#Examples)
-    - [ED, 1D chain](#Ed,-1d-chain)
-    - [ED, Star graph](#Ed,-star-graph)
+    - [ED of 1D chain](#Ed-of-1d-chain)
+    - [ED of Star graph](#Ed-of-star-graph)
     - [Time evolution](#Time-evolution)
 
 ## Installing
@@ -599,7 +599,7 @@ opr=ls.Operator(expression,basis) # Make the operator in symmetry-adapted basis
 
 Here, we will take a look at different examples of `lattice_symmetries` applications.
 
-### ED, 1D chain
+### ED of 1D chain
 
 We will start with the simplest example of exact diagonalization. We will consider Heisenberg chain on 10 sites and diagonalization with the help of symmetries.
 For that, we will combine methods described in the previous sections.
@@ -676,7 +676,7 @@ assert np.isclose(eigenvalues[0], -18.06178542)
 Indeed, the outcome is the same within machine precision. 
 
 
-### ED, Star graph
+### ED of Star graph
 
 In the previous example, we considered a 1D chain, where we knew the graph's symmetries and could code them by hand.
 Sometimes, that task requires much time or is even practically impossible.

From b0a20156b47386d0bd57cd95d563b2b0c500782e Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Tue, 30 Apr 2024 18:23:10 +0200
Subject: [PATCH 51/67] Update tutorial.md

---
 python/tutorial.md | 66 +++++++++++++++++++++++++++-------------------
 1 file changed, 39 insertions(+), 27 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 1572dc2..3ff5827 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -555,33 +555,45 @@ where $n$ is the order of the generator, and $k$ is a natural number.
 However, if the symmetry group is non-abelian, the phases should be chosen in a proper way,
 so that characters of group elements organize a one-dimension representation of the symmetry group.
 
-Now, we are ready to build a symmetry-adapted basis. Th simplest example would be:
-```pycon
-p = ls.Permutation([1,2,0])
-b0 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(0, 1))]) #We specify the phase to be zero. It defines the trivial character
-b0.build()
-```
-
-The order of the permutation is 3, wo we can have 2 other possible sectors:
-```pycon
-b1 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(1, 3))]) #The phase is 1/3, the total phase is 2/3π
-b1.build()
-
-b2 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(2, 3))]) #The phase is 2/3, the total phase is 4/3π
-b2.build()
-```
-When making calculations, one needs to choose which sectors are suitable for the given task.
-
-Let's consider another example, where we will work with the symmetries of expressions:
-```pycon
-e=ls.Expr("σ^x_0 σ^x_1")
-expr=e.on(ig.Graph.Full(5)) # We define our lattice to be a complete graph with 5 vertices
-sym=expr.permutation_group()
-rep=[(p,ls.Rational(0,1)) for p in sym] #The trivial one-dimensional representation 
-#Notice that we can use the whole symmetry group for constructing 1D-representation
-basis = ls.SpinBasis(5, symmetries=rep) 
-basis.build()
-```
+Now, we are ready to build a symmetry-adapted basis. 
+There are two main ways to do that. The first way is to generate the representations by hand, and the second is to use the built-in functions of `lattice-symmetries`.
+
+- 	Manually generated representations:
+
+	The simplest example would be:
+	```pycon
+	p = ls.Permutation([1,2,0]) #translation by 1
+	b0 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(0, 1))]) #We specify the phase to be zero. It defines the trivial character
+	b0.build()
+	```
+
+	The order of the permutation is 3, wo we can have 2 other possible sectors:
+	```pycon
+	b1 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(1, 3))]) #The phase is 1/3, the total phase is 2/3π
+	b1.build()
+
+	b2 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(2, 3))]) #The phase is 2/3, the total phase is 4/3π
+	b2.build()
+	```
+	When making calculations, one needs to choose which sectors are suitable for the given task.
+
+	Let's consider another example, where we will work with the symmetries of expressions:
+	```pycon
+	e=ls.Expr("σ^x_0 σ^x_1")
+	expr=e.on(ig.Graph.Full(5)) # We define our lattice to be a complete graph with 5 vertices
+	sym=expr.permutation_group()
+	rep=[(p,ls.Rational(0,1)) for p in sym] #The trivial one-dimensional representation 
+	#Notice that we can use the whole symmetry group for constructing 1D-representation
+	basis = ls.SpinBasis(5, symmetries=rep) 
+	basis.build()
+	```
+
+	It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
+	Of course, the generators should be consistent.
+
+- 	Functions `hilbert_space_sectors` and `ground_state_sectors`:
+
+	
 
 The last step is to make calculations in symmetry adapted basis. For that we need to create an `Operator` object.
 When making an operator, one needs to be sure that the basis symmetries are consistent with the symmetries of the correspondning expression.

From f7d2a7e04b58a3e14b699a835a77fbb110c0ace0 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <A.Iliasov@science.ru.nl>
Date: Wed, 1 May 2024 17:37:21 +0200
Subject: [PATCH 52/67] Update tutorial.md

---
 python/tutorial.md | 49 +++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 46 insertions(+), 3 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 3ff5827..3185ccb 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -261,7 +261,7 @@ False
 
 ### Basis
 
-In this section, we will consider the generation of basises of various types.
+In this section, we will consider the generation of bases of various types.
 We will not consider symmetry-adapted bases yet; we will do so in the section [Symmetry-adapted basis](#Symmetry-adapted-basis).
 
 
@@ -589,11 +589,54 @@ There are two main ways to do that. The first way is to generate the representat
 	```
 
 	It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
-	Of course, the generators should be consistent.
+	Of course, the generators should be consistent with each other.
 
 - 	Functions `hilbert_space_sectors` and `ground_state_sectors`:
 
+	There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
+	or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`. Let's consider an example of the most straightforward example of 3-sites chain:
+	```pycon
+	import igraph as ig
+
+	e=ls.Expr("σ^x_0 σ^x_1") #A simple expression defined on an edge
+	expr=e.on(ig.Graph.Lattice(dim=[3], circular=True))
+	```
+	
+	We can find the sectors (symmetry adapted bases) using an `Expression`, as we did before.
+	
+	```pycon
+	basis_list=expr.hilbert_space_sectors()
+	```
+	
+	The result is the array of symmetrized bases for each of the possible representations of the maximal Abelian subgroup.
+	It should be noted that this list covers the whole Hilbert space and can be used for the exhaustive search of the relevant sectors.
+	One can check the corresponding representations as follows:
+	```
+	for basis in basis_list:
+		print(basis.symmetries)
+	```
+	
+	Another useful function is `ground_state_sectors`:
+	
+	```pycon
+	basis_list=expr.ground_state_sectors()
+	```
+	The result is the array of symmetrized bases for every one-dimensional representation of the whole symmetry group.
+	The bases do not cover the whole Hilbert space, but they can be used for specific purposes. For example, if the ground state is unique, it for sure lies in one of the given sectors.
+	
+	We can check, that the result is expected by printing the corresponding representations:
+	```
+	for basis in basis_list:
+		print(basis.symmetries)
+	>>>
+	[(Permutation(2), 0), (Permutation(1, 2), 1/2), (Permutation(2)(0, 1), 1/2), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 1/2)]
+	[(Permutation(2), 0), (Permutation(1, 2), 0), (Permutation(2)(0, 1), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 0)]
+	[(Permutation(2), 0), (Permutation(1, 2), 1/2), (Permutation(2)(0, 1), 1/2), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 1/2)]
+	[(Permutation(2), 0), (Permutation(1, 2), 0), (Permutation(2)(0, 1), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 0)]
+	```
+	We have 4 representations in total with commuting parity and translations.
 	
+	It is important to notice that bases are not built yet. Therefore after choosing one of the bases, one should build it using `basis.build()` to make calculations.	
 
 The last step is to make calculations in symmetry adapted basis. For that we need to create an `Operator` object.
 When making an operator, one needs to be sure that the basis symmetries are consistent with the symmetries of the correspondning expression.
@@ -601,7 +644,7 @@ An operator for symmetry-adapted basis is built in the same way as without symme
 
 ```pycon
 Make an expression
-Calculate symmetries
+Calculate symmetries/choose symmetry-adapted basis
 Build symmetry-adapted basis
 
 opr=ls.Operator(expression,basis) # Make the operator in symmetry-adapted basis

From 8d3406fd0e3c25e495d8910a9151a97ace3b4980 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:13:22 +0200
Subject: [PATCH 53/67] Update tutorial.md

---
 python/tutorial.md | 189 ++++++++++++++++++++++++++-------------------
 1 file changed, 109 insertions(+), 80 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 3185ccb..fca31bb 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -19,10 +19,13 @@
 		- [Symmetry as Permutations](#Symmetry-as-permutations)
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry-adapted basis](#Symmetry-adapted-basis)
+			- [Manually generated representations](#Manually-generated-representations)
+			- [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Functions-hilbert-space-sectors-and-ground-state-sectors)
 - [Examples](#Examples)
     - [ED of 1D chain](#Ed-of-1d-chain)
     - [ED of Star graph](#Ed-of-star-graph)
     - [Time evolution](#Time-evolution)
+	- [Spin structure factor](#Spin-structure-factor)
 
 ## Installing
 
@@ -558,85 +561,101 @@ so that characters of group elements organize a one-dimension representation of
 Now, we are ready to build a symmetry-adapted basis. 
 There are two main ways to do that. The first way is to generate the representations by hand, and the second is to use the built-in functions of `lattice-symmetries`.
 
-- 	Manually generated representations:
-
-	The simplest example would be:
-	```pycon
-	p = ls.Permutation([1,2,0]) #translation by 1
-	b0 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(0, 1))]) #We specify the phase to be zero. It defines the trivial character
-	b0.build()
-	```
-
-	The order of the permutation is 3, wo we can have 2 other possible sectors:
-	```pycon
-	b1 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(1, 3))]) #The phase is 1/3, the total phase is 2/3π
-	b1.build()
-
-	b2 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(2, 3))]) #The phase is 2/3, the total phase is 4/3π
-	b2.build()
-	```
-	When making calculations, one needs to choose which sectors are suitable for the given task.
-
-	Let's consider another example, where we will work with the symmetries of expressions:
-	```pycon
-	e=ls.Expr("σ^x_0 σ^x_1")
-	expr=e.on(ig.Graph.Full(5)) # We define our lattice to be a complete graph with 5 vertices
-	sym=expr.permutation_group()
-	rep=[(p,ls.Rational(0,1)) for p in sym] #The trivial one-dimensional representation 
-	#Notice that we can use the whole symmetry group for constructing 1D-representation
-	basis = ls.SpinBasis(5, symmetries=rep) 
-	basis.build()
-	```
-
-	It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
-	Of course, the generators should be consistent with each other.
-
-- 	Functions `hilbert_space_sectors` and `ground_state_sectors`:
-
-	There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
-	or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`. Let's consider an example of the most straightforward example of 3-sites chain:
-	```pycon
-	import igraph as ig
-
-	e=ls.Expr("σ^x_0 σ^x_1") #A simple expression defined on an edge
-	expr=e.on(ig.Graph.Lattice(dim=[3], circular=True))
-	```
-	
-	We can find the sectors (symmetry adapted bases) using an `Expression`, as we did before.
-	
-	```pycon
-	basis_list=expr.hilbert_space_sectors()
-	```
-	
-	The result is the array of symmetrized bases for each of the possible representations of the maximal Abelian subgroup.
-	It should be noted that this list covers the whole Hilbert space and can be used for the exhaustive search of the relevant sectors.
-	One can check the corresponding representations as follows:
-	```
-	for basis in basis_list:
-		print(basis.symmetries)
-	```
-	
-	Another useful function is `ground_state_sectors`:
-	
-	```pycon
-	basis_list=expr.ground_state_sectors()
-	```
-	The result is the array of symmetrized bases for every one-dimensional representation of the whole symmetry group.
-	The bases do not cover the whole Hilbert space, but they can be used for specific purposes. For example, if the ground state is unique, it for sure lies in one of the given sectors.
-	
-	We can check, that the result is expected by printing the corresponding representations:
-	```
-	for basis in basis_list:
-		print(basis.symmetries)
-	>>>
-	[(Permutation(2), 0), (Permutation(1, 2), 1/2), (Permutation(2)(0, 1), 1/2), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 1/2)]
-	[(Permutation(2), 0), (Permutation(1, 2), 0), (Permutation(2)(0, 1), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 0)]
-	[(Permutation(2), 0), (Permutation(1, 2), 1/2), (Permutation(2)(0, 1), 1/2), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 1/2)]
-	[(Permutation(2), 0), (Permutation(1, 2), 0), (Permutation(2)(0, 1), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 0)]
-	```
-	We have 4 representations in total with commuting parity and translations.
-	
-	It is important to notice that bases are not built yet. Therefore after choosing one of the bases, one should build it using `basis.build()` to make calculations.	
+##### Manually generated representations:
+
+The simplest example would be:
+```pycon
+p = ls.Permutation([1,2,0]) #translation by 1
+b0 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(0, 1))]) #We specify the phase to be zero. It defines the trivial character
+b0.build()
+```
+
+The order of the permutation is 3, wo we can have 2 other possible sectors:
+```pycon
+b1 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(1, 3))]) #The phase is 1/3, the total phase is 2/3π
+b1.build()
+
+b2 = ls.SpinBasis(3, symmetries=[(p, ls.Rational(2, 3))]) #The phase is 2/3, the total phase is 4/3π
+b2.build()
+```
+When making calculations, one needs to choose which sectors are suitable for the given task.
+
+Let's consider another example, where we will work with the symmetries of expressions:
+```pycon
+e=ls.Expr("σ^x_0 σ^x_1")
+expr=e.on(ig.Graph.Full(5)) # We define our lattice to be a complete graph with 5 vertices
+sym=expr.permutation_group()
+rep=[(p,ls.Rational(0,1)) for p in sym] #The trivial one-dimensional representation 
+#Notice that we can use the whole symmetry group for constructing 1D-representation
+basis = ls.SpinBasis(5, symmetries=rep) 
+basis.build()
+```
+
+It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
+Of course, the generators should be consistent with each other.
+
+##### Functions `hilbert_space_sectors` and `ground_state_sectors`:
+
+There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
+or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.
+It is important to note that the function also calculate bases for different number of particles/hamming weight and spin inversion, so that the whole Hilbert space is covered. 
+Let's consider an example of the most straightforward example of 3-sites chain:
+```pycon
+import igraph as ig
+
+e=ls.Expr("σ^x_0 σ^x_1") #A simple expression defined on an edge
+expr=e.on(ig.Graph.Lattice(dim=[3], circular=True))
+```
+
+We can find the sectors (symmetry adapted bases) using an `Expression`, as we did before.
+
+```pycon
+basis_list=expr.hilbert_space_sectors()
+```
+
+The result is the array of symmetrized bases for each of the possible representations of the maximal Abelian subgroup, hamming weight and spin inversion.
+It should be noted that this list covers the whole Hilbert space and can be used for the exhaustive search of the relevant sectors.
+One can check the corresponding representations as follows:
+```
+for basis in basis_list:
+	print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible sectors depending on translation sector, hamming weight and spin inversion
+>>>
+[(Permutation(2), 0), (Permutation(0, 1, 2), 1/3), (Permutation(0, 2, 1), 2/3)] None -1
+[(Permutation(2), 0), (Permutation(0, 1, 2), 2/3), (Permutation(0, 2, 1), 1/3)] None -1
+[(Permutation(2), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0)] None -1
+[(Permutation(2), 0), (Permutation(0, 1, 2), 1/3), (Permutation(0, 2, 1), 2/3)] None 1
+[(Permutation(2), 0), (Permutation(0, 1, 2), 2/3), (Permutation(0, 2, 1), 1/3)] None 1
+[(Permutation(2), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0)] None 1
+```
+We see that the hamming weight is not defined, however, the spin inversion is possible, and for each value of spin inversion we have three one-dimensional representations of the translation group.
+We can also dcheck the expression indeed does not conserve the number of particles, so it is impossible to define basis for the chosen hamming weight:
+
+```pycon
+e.conserves_number_particles
+>>> False
+```
+
+Another useful function is `ground_state_sectors`:
+
+```pycon
+basis_list=expr.ground_state_sectors()
+```
+The result is the array of symmetrized bases for every one-dimensional representation of the whole symmetry group and possible values of particles/hamming weight and spin inversion.
+The bases do not cover the whole Hilbert space, but they can be used for specific purposes. For example, if the ground state is unique, it for sure lies in one of the given sectors.
+
+We can check, that the result is expected by printing the corresponding representations:
+```
+for basis in basis_list:
+	print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible sectors depending on translation sector, hamming weight and spin inversion
+>>>
+[(Permutation(2), 0), (Permutation(1, 2), 1/2), (Permutation(2)(0, 1), 1/2), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 1/2)] None -1
+[(Permutation(2), 0), (Permutation(1, 2), 0), (Permutation(2)(0, 1), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 0)] None -1
+[(Permutation(2), 0), (Permutation(1, 2), 1/2), (Permutation(2)(0, 1), 1/2), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 1/2)] None 1
+[(Permutation(2), 0), (Permutation(1, 2), 0), (Permutation(2)(0, 1), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 0)] None 1
+```
+We have 2 representations with commuting parity and translations for each values of spin inversion.
+
+It is important to notice that bases are not built yet. Therefore after choosing one of the bases, one should build it using `basis.build()` to make calculations.	
 
 The last step is to make calculations in symmetry adapted basis. For that we need to create an `Operator` object.
 When making an operator, one needs to be sure that the basis symmetries are consistent with the symmetries of the correspondning expression.
@@ -786,4 +805,14 @@ $$
 
 where $a=1/(2-\epsilon)$, $b=(E_g+E_s)/2$, and $\epsilon$ is a safety factor to make the spectrum of $H'$ be within $[-1,1]$.
 
-We will consider the basic version of the algorithm, where instead of usual `numpy` arrays we will use `lattice_symmetries` operators.
\ No newline at end of file
+We will consider the basic version of the algorithm, where instead of usual `numpy` arrays we will use `lattice_symmetries` operators.
+
+### Spin structure factor
+
+We can use `lattice-symmetries` to calculate static spin structure factor:
+
+$$
+S({\bf q})=\frac{1}{N}\sum e^{i{\bf q}\cdot ({\bf r_i}-{\bf r_j})}\langle S_i S_j \rangle
+$$
+
+The idea is to use the symmetries of the system to make calculations faster.
\ No newline at end of file

From 1076dfee32fc2bc04fb1387f99a55e538c7bf05b Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:16:18 +0200
Subject: [PATCH 54/67] Update tutorial.md

---
 python/tutorial.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index fca31bb..a3f269e 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -20,7 +20,7 @@
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry-adapted basis](#Symmetry-adapted-basis)
 			- [Manually generated representations](#Manually-generated-representations)
-			- [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Functions-hilbert-space-sectors-and-ground-state-sectors)
+			- [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Functions-`hilbert-space-sectors`-and-`ground-state-sectors`)
 - [Examples](#Examples)
     - [ED of 1D chain](#Ed-of-1d-chain)
     - [ED of Star graph](#Ed-of-star-graph)
@@ -561,7 +561,7 @@ so that characters of group elements organize a one-dimension representation of
 Now, we are ready to build a symmetry-adapted basis. 
 There are two main ways to do that. The first way is to generate the representations by hand, and the second is to use the built-in functions of `lattice-symmetries`.
 
-##### Manually generated representations:
+##### Manually generated representations
 
 The simplest example would be:
 ```pycon
@@ -594,7 +594,7 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### Functions `hilbert_space_sectors` and `ground_state_sectors`:
+##### Functions `hilbert_space_sectors` and `ground_state_sectors`
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.
@@ -628,7 +628,7 @@ for basis in basis_list:
 [(Permutation(2), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0)] None 1
 ```
 We see that the hamming weight is not defined, however, the spin inversion is possible, and for each value of spin inversion we have three one-dimensional representations of the translation group.
-We can also dcheck the expression indeed does not conserve the number of particles, so it is impossible to define basis for the chosen hamming weight:
+We can also check the expression indeed does not conserve the number of particles, so it is impossible to define basis for the chosen hamming weight:
 
 ```pycon
 e.conserves_number_particles

From 83fdcf50209796bf1c4a6b96ead469743614f01f Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:22:24 +0200
Subject: [PATCH 55/67] Update tutorial.md

---
 python/tutorial.md | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index a3f269e..e56db77 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -20,7 +20,7 @@
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry-adapted basis](#Symmetry-adapted-basis)
 			- [Manually generated representations](#Manually-generated-representations)
-			- [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Functions-`hilbert-space-sectors`-and-`ground-state-sectors`)
+			- [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Functions-`hilbert_space_sectors`-and-`ground_state_sectors`)
 - [Examples](#Examples)
     - [ED of 1D chain](#Ed-of-1d-chain)
     - [ED of Star graph](#Ed-of-star-graph)
@@ -618,7 +618,8 @@ It should be noted that this list covers the whole Hilbert space and can be used
 One can check the corresponding representations as follows:
 ```
 for basis in basis_list:
-	print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible sectors depending on translation sector, hamming weight and spin inversion
+	#Print all possible sectors depending on translation sector, hamming weight and spin inversion
+	print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) 
 >>>
 [(Permutation(2), 0), (Permutation(0, 1, 2), 1/3), (Permutation(0, 2, 1), 2/3)] None -1
 [(Permutation(2), 0), (Permutation(0, 1, 2), 2/3), (Permutation(0, 2, 1), 1/3)] None -1
@@ -627,7 +628,7 @@ for basis in basis_list:
 [(Permutation(2), 0), (Permutation(0, 1, 2), 2/3), (Permutation(0, 2, 1), 1/3)] None 1
 [(Permutation(2), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0)] None 1
 ```
-We see that the hamming weight is not defined, however, the spin inversion is possible, and for each value of spin inversion we have three one-dimensional representations of the translation group.
+We see that the hamming weight is not defined. However, the spin inversion is possible, and for each value of spin inversion we have three one-dimensional representations of the translation group.
 We can also check the expression indeed does not conserve the number of particles, so it is impossible to define basis for the chosen hamming weight:
 
 ```pycon
@@ -646,7 +647,8 @@ The bases do not cover the whole Hilbert space, but they can be used for specifi
 We can check, that the result is expected by printing the corresponding representations:
 ```
 for basis in basis_list:
-	print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible sectors depending on translation sector, hamming weight and spin inversion
+	#Print all possible sectors depending on translation sector, hamming weight and spin inversion
+	print(basis.symmetries, basis.hamming_weight, basis.spin_inversion)
 >>>
 [(Permutation(2), 0), (Permutation(1, 2), 1/2), (Permutation(2)(0, 1), 1/2), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 1/2)] None -1
 [(Permutation(2), 0), (Permutation(1, 2), 0), (Permutation(2)(0, 1), 0), (Permutation(0, 1, 2), 0), (Permutation(0, 2, 1), 0), (Permutation(0, 2), 0)] None -1

From 251264b39da5746cb91177a6a345b32e847394d5 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:25:11 +0200
Subject: [PATCH 56/67] Update tutorial.md

---
 python/tutorial.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index e56db77..cb5f1c8 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -20,7 +20,7 @@
 		- [Symmetry from Expressions](#Symmetry-from-expressions)
 		- [Symmetry-adapted basis](#Symmetry-adapted-basis)
 			- [Manually generated representations](#Manually-generated-representations)
-			- [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Functions-`hilbert_space_sectors`-and-`ground_state_sectors`)
+			- [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Sectors-section)
 - [Examples](#Examples)
     - [ED of 1D chain](#Ed-of-1d-chain)
     - [ED of Star graph](#Ed-of-star-graph)
@@ -594,7 +594,7 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### Functions `hilbert_space_sectors` and `ground_state_sectors`
+##### Functions `hilbert_space_sectors` and `ground_state_sectors` <a name="Sectors-section"></a>
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.

From 6a807079943e4482139fd776addeb8ed26b9702a Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:26:03 +0200
Subject: [PATCH 57/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index cb5f1c8..cb0fc87 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -594,7 +594,7 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### Functions `hilbert_space_sectors` and `ground_state_sectors` <a name="Sectors-section"></a>
+##### Functions `hilbert_space_sectors` and `ground_state_sectors` <a id="Sectors-section"></a>
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.

From 5815d6d9124f5b215fefaea1a0ecaa598dc3dc41 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:28:32 +0200
Subject: [PATCH 58/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index cb0fc87..8c55b57 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -594,7 +594,7 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### Functions `hilbert_space_sectors` and `ground_state_sectors` <a id="Sectors-section"></a>
+##### [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Sectors-section")
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.

From 5eb53bf6196da34ebf24f19588ecf19daea27469 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:29:14 +0200
Subject: [PATCH 59/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 8c55b57..90a61fe 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -594,7 +594,7 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Sectors-section")
+##### [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Sectors-section)
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.

From c93ecc86cd53c32e657fcb1eedebf57ce44571b6 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:35:25 +0200
Subject: [PATCH 60/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 90a61fe..5224580 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -594,7 +594,7 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### [Functions `hilbert_space_sectors` and `ground_state_sectors`](#Sectors-section)
+##### <a name="Sectors-section"></a> Functions `hilbert_space_sectors` and `ground_state_sectors`
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.

From 18bcf2b9ec44e93ec4696ecfe748a90f563d36e7 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:35:58 +0200
Subject: [PATCH 61/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 5224580..80fcb54 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -594,7 +594,7 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### <a name="Sectors-section"></a> Functions `hilbert_space_sectors` and `ground_state_sectors`
+##### <a name="#Sectors-section"></a> Functions `hilbert_space_sectors` and `ground_state_sectors`
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.

From 9351895207865b03ca3603d2d367dc02eb0d287b Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:36:23 +0200
Subject: [PATCH 62/67] Update tutorial.md

---
 python/tutorial.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index 80fcb54..aa4229e 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -594,7 +594,7 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### <a name="#Sectors-section"></a> Functions `hilbert_space_sectors` and `ground_state_sectors`
+##### <a id="#Sectors-section"></a> Functions `hilbert_space_sectors` and `ground_state_sectors`
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.

From 60b964276d32a6a1dc466d71bfbe18e3d5f92ff7 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Thu, 16 May 2024 15:38:11 +0200
Subject: [PATCH 63/67] Update tutorial.md

---
 python/tutorial.md | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/python/tutorial.md b/python/tutorial.md
index aa4229e..9be4941 100644
--- a/python/tutorial.md
+++ b/python/tutorial.md
@@ -594,7 +594,8 @@ basis.build()
 It should be noticed, that one can define a representation by any of its subset, and `lattice_symmetry` will generate the whole representation by the given generators.
 Of course, the generators should be consistent with each other.
 
-##### <a id="#Sectors-section"></a> Functions `hilbert_space_sectors` and `ground_state_sectors`
+<a name="Sectors-section"></a>
+##### Functions `hilbert_space_sectors` and `ground_state_sectors`
 
 There is a possibility to generate the symmetry-adapted basis in `lattice_symmetries automatically.` One can find all sectors of the maximal abelian subgroup using `hilbert_space_sectors`,
 or one can find all one-dimensional representations of the whole symmetry group by `ground_state_sectors`.

From fd65606874b2a0cdca5a020764b1d2d4db54efd4 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Wed, 17 Jul 2024 17:28:44 +0200
Subject: [PATCH 64/67] Create structure_factor.py

---
 python/example/structure_factor.py | 183 +++++++++++++++++++++++++++++
 1 file changed, 183 insertions(+)
 create mode 100644 python/example/structure_factor.py

diff --git a/python/example/structure_factor.py b/python/example/structure_factor.py
new file mode 100644
index 0000000..414b428
--- /dev/null
+++ b/python/example/structure_factor.py
@@ -0,0 +1,183 @@
+import lattice_symmetries as ls
+import numpy as np
+import operator
+import scipy.sparse.linalg
+import scipy.linalg as linalg
+import igraph as ig
+#import matplotlib.pyplot as plt
+
+def thermal_average(observable, H, T):
+    exp_H=linalg.expm(-H/T) #calculate matrix exponent
+    density_matrix=exp_H/np.trace(exp_H)
+    return np.trace(observable@density_matrix)
+
+def square_distance(r, L): #here we assume that L is odd
+    if r>L/2:
+        return r-L
+    else:
+        return r
+    
+def structure_factor_nosym(qx,qy,cor_array, coords, L, N_sites): #calculate structure factor without symmetries for specific values of q
+    sf=0
+    for i in range(N_sites):
+        rx=square_distance(coords[i][0],L)
+        ry=square_distance(coords[i][1],L)
+        sf += np.exp(1j*(qx*rx+qy*ry))*cor_array[i]/N_sites #here we use that the coordinates of the site with index 0 are (0,0)
+        #print(i, rx,ry, np.exp(1j*(qx*coords[i][0]+qy*coords[i][1])), sf)
+        #print(coords[i][0],coords[i][1])
+    return np.real(sf)
+
+def structure_factor_sym(qx,qy,sym_cor_array,coords, orbit_array, L,N_sites):
+#calculate structure factor using symmetries for specific values of q
+    sf=0
+    for orb_ind in range(len(orbit_array)):
+        for i in orbit_array[orb_ind]:
+            rx=square_distance(coords[i][0],L)
+            ry=square_distance(coords[i][1],L)
+            sf += np.exp(1j*(qx*rx+qy*ry))*sym_cor_array[orb_ind]/N_sites #here we also use that the coordinates of the site with index 0 are (0,0)
+        #print(i, rx,ry, np.exp(1j*(qx*coords[i][0]+qy*coords[i][1])), sf)
+        #print(coords[i][0],coords[i][1])
+    return np.real(sf)
+
+
+def main():
+    
+    T=1 #temperature
+    Nq=60 #discretization of q-space
+    qxs = np.linspace(-3, 3, Nq)
+    qys = np.linspace(-3, 3, Nq)
+    L=3
+
+    #define hamiltonian
+    basic_expr = ls.Expr("Sx0 Sx1 + Sy0 Sy1 + Sz0 Sz1") # this basic expression is the building block of our Hamiltonian and correlation operators
+    lattice=ig.Graph.Lattice(dim=[L, L])
+    H_expr=basic_expr.on(lattice) # Hamiltonian expression, on square 3x3
+    N_sites=H_expr.number_sites
+    #For calculations we will also need coordinates of the vertices. Let's use ig.layout
+    coords=lattice.layout_grid()
+
+    #At first let's calculate structure factor without using any symmetry.
+    #For that, we need to calculate correlations functions.
+    full_basis=H_expr.full_basis()
+    full_basis.build()
+    H_op=ls.Operator(H_expr,full_basis)
+    H_matrix=H_op.to_dense()#Make numpy matrix for matrix exponentiation
+    cor_array=[]
+    for i in range(N_sites):
+        cor_expr=basic_expr.replace_indices({1: i}) #make operator S_0S_i
+        cor_op=ls.Operator(cor_expr,full_basis)
+        cor_matrix=cor_op.to_dense() #make numpy matrix
+        cor_array.append(thermal_average(cor_matrix,H_matrix,T)) #calculate observable and add to array
+    print("corrs", cor_array)
+
+    SF=np.zeros((Nq,Nq))
+    for i in range(Nq):
+        for j in range(Nq):
+            SF[i,j]=structure_factor_nosym(qxs[i],qys[j],cor_array,coords, L,N_sites)
+    #print(SF)
+    np.save("nosym_cors.npy", SF)
+    
+    #plot structure factor
+    #h = plt.contourf(x, y, zs)
+    #plt.axis('scaled')
+    #plt.colorbar()
+    #plt.savefig("structure_factor_no_symmetries.png")
+    #plt.close()
+
+    #Now, let's exploit the power of symmetries! 
+    #We will divide our space on symmetric sectors, and use the fact that structure factor will be the same for some of the sectors.
+    #In particular, we are interested in the symmetries generated by the point group.
+    #Let's find it. We are looking for the subgroup, which fixes index 0.
+    
+    symmetries=H_expr.permutation_group()
+    point_group=[]
+    for symmetry in symmetries:
+        list_form=symmetry.list()
+        if list_form[0]==0:
+            point_group.append(symmetry)
+    #print("point group", point_group)
+    #The next step is to find orbits of non-zero indices. 
+    #If two indices can be transformed to each other by an element of a point group, they correspond to the same correlation function C_{0,i}.
+    index_array=[i for i in range(1,N_sites)]
+    #print("index array", index_array)
+    orbit_array=[[0]] #zero index is always here and unique
+    while index_array!=[]:
+        index=index_array[0]
+        index_array.remove(index)
+        orbit=[index]
+        #print("orbit", orbit)
+        for symmetry in point_group:
+            list_form=symmetry.list()
+            new_index=list_form[index]
+            #print("symmetry", symmetry, "new index", new_index)
+            if new_index not in orbit:
+                orbit.append(new_index)
+                index_array.remove(new_index)
+        orbit_array.append(orbit)
+    #print("orbit array", orbit_array)
+    #We are ready to calculate structure factor!
+    #At first we calculate correlation function for an element of an orbit, then add the correct phase factor
+    sym_cor_array=[]
+    for orbit in orbit_array:
+        cor_expr=basic_expr.replace_indices({1: orbit[0]}) #make operator S_0S_i
+        cor_op=ls.Operator(cor_expr,full_basis)
+        cor_matrix=cor_op.to_dense() #make numpy matrix
+        sym_cor_array.append(thermal_average(cor_matrix,H_matrix,T)) #calculate observable and add to array
+    print("sym corrs", sym_cor_array)
+
+
+    SF_sym=np.zeros((Nq,Nq))
+    for i in range(Nq):
+        for j in range(Nq):
+            SF_sym[i,j]=structure_factor_sym(qxs[i],qys[j],sym_cor_array,coords, orbit_array, L,N_sites)
+    #print(SF)
+    np.save("sym_cors.npy", SF_sym)
+
+    #Now, let's reverse the order of the symmetry calculations.
+    #At first we will make a summation over different translation sectors using the translations and the action of point group.
+    #The idea is to construct the observables invariant under the action of the point group so that there is the exact correspondence between different sectors of translation group.
+    #The symmetric observables are made from already calculated orbits of C_{0,i}
+
+    hf_bases=[] #array of bases with specific hamming weight
+    translation_bases=H_expr.hilbert_space_sectors()
+    for basis in translation_bases:
+        #print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible symmetries. We have various hamming weights as well
+        if basis.hamming_weight==5:
+            hf_bases.append(basis)
+            print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible symmetries. We have various hamming weights as well
+    translation_group=H_expr.abelian_permutation_group()
+
+    symmetric_exprs=[]
+    for orbit in orbit_array:
+        if len(orbit)>1:
+            print("orbit", orbit, "basic_expr", basic_expr)
+            point_expr=ls.Expr("0.0 Sx0 Sx1")
+            for j in orbit:
+                point_expr+=basic_expr.replace_indices({1: j})#make operator S_0S_i
+            print("point_expr", point_expr)
+            cor_expr=point_expr
+            for translation in translation_group:
+                print(T)
+                cor_expr+=point_expr.replace_indices({n: n^translation for n in range(9)})
+            print(cor_expr)
+            symmetric_exprs.append(cor_expr)
+    print("sym exprs", symmetric_exprs)
+
+    #Now, let's make a check, which representations will give the same answer. Let's at first focus on a specifiv value of hamming weight.
+    for basis in hf_bases:
+        print(basis.symmetries)
+        basis.build()
+        cor_op_sym=ls.Operator(symmetric_exprs[0],basis)
+        cor_matrix_sym=cor_op.to_dense() #make numpy matrix
+        H_op_sym=ls.Operator(H_expr,basis)
+        H_matrix_sym=H_op_sym.to_dense()
+        print(cor_matrix_sym)
+        print(H_matrix_sym)
+        corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
+        print("corr", corr)
+
+
+    
+
+
+main()

From ef966358b335bf585bde470fdce2f0dbb2fd63e2 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Wed, 2 Oct 2024 18:44:08 +0200
Subject: [PATCH 65/67] Update structure_factor.py

---
 python/example/structure_factor.py | 86 +++++++++++++++++++++++-------
 1 file changed, 66 insertions(+), 20 deletions(-)

diff --git a/python/example/structure_factor.py b/python/example/structure_factor.py
index 414b428..2baa641 100644
--- a/python/example/structure_factor.py
+++ b/python/example/structure_factor.py
@@ -23,8 +23,6 @@ def structure_factor_nosym(qx,qy,cor_array, coords, L, N_sites): #calculate stru
         rx=square_distance(coords[i][0],L)
         ry=square_distance(coords[i][1],L)
         sf += np.exp(1j*(qx*rx+qy*ry))*cor_array[i]/N_sites #here we use that the coordinates of the site with index 0 are (0,0)
-        #print(i, rx,ry, np.exp(1j*(qx*coords[i][0]+qy*coords[i][1])), sf)
-        #print(coords[i][0],coords[i][1])
     return np.real(sf)
 
 def structure_factor_sym(qx,qy,sym_cor_array,coords, orbit_array, L,N_sites):
@@ -35,8 +33,6 @@ def structure_factor_sym(qx,qy,sym_cor_array,coords, orbit_array, L,N_sites):
             rx=square_distance(coords[i][0],L)
             ry=square_distance(coords[i][1],L)
             sf += np.exp(1j*(qx*rx+qy*ry))*sym_cor_array[orb_ind]/N_sites #here we also use that the coordinates of the site with index 0 are (0,0)
-        #print(i, rx,ry, np.exp(1j*(qx*coords[i][0]+qy*coords[i][1])), sf)
-        #print(coords[i][0],coords[i][1])
     return np.real(sf)
 
 
@@ -95,26 +91,21 @@ def main():
         list_form=symmetry.list()
         if list_form[0]==0:
             point_group.append(symmetry)
-    #print("point group", point_group)
     #The next step is to find orbits of non-zero indices. 
     #If two indices can be transformed to each other by an element of a point group, they correspond to the same correlation function C_{0,i}.
     index_array=[i for i in range(1,N_sites)]
-    #print("index array", index_array)
     orbit_array=[[0]] #zero index is always here and unique
     while index_array!=[]:
         index=index_array[0]
         index_array.remove(index)
         orbit=[index]
-        #print("orbit", orbit)
         for symmetry in point_group:
             list_form=symmetry.list()
             new_index=list_form[index]
-            #print("symmetry", symmetry, "new index", new_index)
             if new_index not in orbit:
                 orbit.append(new_index)
                 index_array.remove(new_index)
         orbit_array.append(orbit)
-    #print("orbit array", orbit_array)
     #We are ready to calculate structure factor!
     #At first we calculate correlation function for an element of an orbit, then add the correct phase factor
     sym_cor_array=[]
@@ -130,7 +121,6 @@ def main():
     for i in range(Nq):
         for j in range(Nq):
             SF_sym[i,j]=structure_factor_sym(qxs[i],qys[j],sym_cor_array,coords, orbit_array, L,N_sites)
-    #print(SF)
     np.save("sym_cors.npy", SF_sym)
 
     #Now, let's reverse the order of the symmetry calculations.
@@ -144,38 +134,94 @@ def main():
         #print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible symmetries. We have various hamming weights as well
         if basis.hamming_weight==5:
             hf_bases.append(basis)
-            print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible symmetries. We have various hamming weights as well
+            #print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible symmetries. We have various hamming weights as well
     translation_group=H_expr.abelian_permutation_group()
 
     symmetric_exprs=[]
     for orbit in orbit_array:
         if len(orbit)>1:
-            print("orbit", orbit, "basic_expr", basic_expr)
+            #print("orbit", orbit, "basic_expr", basic_expr)
             point_expr=ls.Expr("0.0 Sx0 Sx1")
             for j in orbit:
-                point_expr+=basic_expr.replace_indices({1: j})#make operator S_0S_i
-            print("point_expr", point_expr)
+                point_expr+=basic_expr.replace_indices({1: j})#make operator S_0S_i as it occurs in spin factor
+            #print("point_expr", point_expr)
             cor_expr=point_expr
             for translation in translation_group:
-                print(T)
+                #print("translation", translation)
                 cor_expr+=point_expr.replace_indices({n: n^translation for n in range(9)})
-            print(cor_expr)
+            #print("cor expr", cor_expr)
             symmetric_exprs.append(cor_expr)
     print("sym exprs", symmetric_exprs)
 
-    #Now, let's make a check, which representations will give the same answer. Let's at first focus on a specifiv value of hamming weight.
+    #Now, let's make a check, which representations will give the same answer. Let's at first focus on a specific value of hamming weight.
     for basis in hf_bases:
         print(basis.symmetries)
         basis.build()
         cor_op_sym=ls.Operator(symmetric_exprs[0],basis)
-        cor_matrix_sym=cor_op.to_dense() #make numpy matrix
+        cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
         H_op_sym=ls.Operator(H_expr,basis)
         H_matrix_sym=H_op_sym.to_dense()
-        print(cor_matrix_sym)
-        print(H_matrix_sym)
+        #print(cor_matrix_sym.shape)
+        #print(H_matrix_sym.shape)
         corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
         print("corr", corr)
 
+    #We see that we have only 3 independent values of the correlator. We can do calculate them in a more automatic way.
+    #The representations of translation group related to each other by an action of point group (lie in the same conjugacy class) give the same correlators.
+    #Let's select the representatives of the bases lying in the same conjugacy class.
+    bases_classes=[]
+    #At first we will ort all symmetries, so that we will be able to compare them
+    #for basis in hf_bases:
+    #    basis.symmetries=sorted(basis.symmetries, key=lambda x: list(x[0].array_form))
+    bases_count=hf_bases.copy()         
+    #print("point group", point_group)
+    
+    for basis in hf_basis:
+        #print("basis_symmetry", basis.symmetries) #Print the representation corresponding to the given basis
+        tb=basis.symmetries.copy()
+        conj_class=[]
+        for g in point_group:
+            for i in range(len(tb)):
+                tb[i]=(g*tb[i][0]*(g**(-1)),tb[i][1])
+            tb=sorted(tb, key=lambda x: list(x[0].array_form))
+            #print("tb", tb)
+            for bs in bases_count:
+                symsort=sorted(bs.symmetries, key=lambda x: list(x[0].array_form))
+                if symsort==tb:
+                    check=True #We check that we do not overcount bases since some elements of the point group can give the same representations
+                    for bs_check in basis_count:
+                        symsort_check=sorted(bs_check.symmetries, key=lambda x: list(x[0].array_form))
+                        if symsort_check==tb:
+                            check=False
+                            #print("check")
+                            continue
+                    if check==True:
+                        conj_class.append(bs)
+                        #print("symsort", symsort)
+                        #print("tb", tb)
+                        #print("symmetries", bs.symmetries)
+        #print(conj_class)
+        if len(conj_class)>0:
+            bases_classes.append(conj_class)
+    
+    #Now, let's check if these conjugacy classes indeed give the same results
+    for base_class in bases_classes:
+        print("base_class", base_class)
+        for bs in base_class:
+            print(bs.symmetries)
+            bs.build()
+            cor_op_sym=ls.Operator(symmetric_exprs[0],bs)
+            cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
+            H_op_sym=ls.Operator(H_expr,bs)
+            H_matrix_sym=H_op_sym.to_dense()
+            #print(cor_matrix_sym.shape)
+            #print(H_matrix_sym.shape)
+            corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
+            print("corr", corr)
+
+
+
+
 
     
 

From 8e87630e647c9691986419182067f7a710132a21 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Tue, 22 Oct 2024 18:33:12 +0200
Subject: [PATCH 66/67] Update structure_factor.py

---
 python/example/structure_factor.py | 203 +++++++++++++++++++----------
 1 file changed, 134 insertions(+), 69 deletions(-)

diff --git a/python/example/structure_factor.py b/python/example/structure_factor.py
index 2baa641..cd7d0b7 100644
--- a/python/example/structure_factor.py
+++ b/python/example/structure_factor.py
@@ -6,6 +6,10 @@
 import igraph as ig
 #import matplotlib.pyplot as plt
 
+def weighted_trace(H,T):
+    exp_H=linalg.expm(-H/T) #calculate matrix exponent
+    return np.trace(exp_H)
+
 def thermal_average(observable, H, T):
     exp_H=linalg.expm(-H/T) #calculate matrix exponent
     density_matrix=exp_H/np.trace(exp_H)
@@ -36,9 +40,31 @@ def structure_factor_sym(qx,qy,sym_cor_array,coords, orbit_array, L,N_sites):
     return np.real(sf)
 
 
+def conjugate_classes(bases, point_group):
+    bases_classes=[]
+    bases_count=bases.copy()          
+    while bases_count!=[]:
+        basis = bases_count[0]
+        conj_class=[basis]
+        bases_count.remove(basis)
+        del_array=[]
+        for g in point_group:
+            tb=basis.symmetries.copy()
+            for i in range(len(tb)):
+                tb[i]=(g*tb[i][0]*(g**(-1)),tb[i][1])
+            tb=sorted(tb, key=lambda x: list(x[0].array_form))
+            for bs in bases_count:
+                if bs.symmetries==tb and (bs not in conj_class):
+                    conj_class.append(bs)
+                    del_array.append(bs)
+        for del_bs in del_array:
+            bases_count.remove(del_bs)
+        bases_classes.append(conj_class)
+        return bases_classes
+
 def main():
     
-    T=1 #temperature
+    T= 1 #temperature
     Nq=60 #discretization of q-space
     qxs = np.linspace(-3, 3, Nq)
     qys = np.linspace(-3, 3, Nq)
@@ -91,6 +117,7 @@ def main():
         list_form=symmetry.list()
         if list_form[0]==0:
             point_group.append(symmetry)
+    #print("point group", point_group)
     #The next step is to find orbits of non-zero indices. 
     #If two indices can be transformed to each other by an element of a point group, they correspond to the same correlation function C_{0,i}.
     index_array=[i for i in range(1,N_sites)]
@@ -126,98 +153,136 @@ def main():
     #Now, let's reverse the order of the symmetry calculations.
     #At first we will make a summation over different translation sectors using the translations and the action of point group.
     #The idea is to construct the observables invariant under the action of the point group so that there is the exact correspondence between different sectors of translation group.
-    #The symmetric observables are made from already calculated orbits of C_{0,i}
-
+    #We start with specific hamming weight and generate symmetric bases
     hf_bases=[] #array of bases with specific hamming weight
     translation_bases=H_expr.hilbert_space_sectors()
     for basis in translation_bases:
-        #print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible symmetries. We have various hamming weights as well
         if basis.hamming_weight==5:
             hf_bases.append(basis)
-            #print(basis.symmetries, basis.hamming_weight, basis.spin_inversion) #Print all possible symmetries. We have various hamming weights as well
     translation_group=H_expr.abelian_permutation_group()
-
+    
+    #Make symmetrized expressions for the correlators
+    #The symmetric observables are made from already calculated orbits of C_{0,i}
     symmetric_exprs=[]
     for orbit in orbit_array:
-        if len(orbit)>1:
-            #print("orbit", orbit, "basic_expr", basic_expr)
-            point_expr=ls.Expr("0.0 Sx0 Sx1")
-            for j in orbit:
-                point_expr+=basic_expr.replace_indices({1: j})#make operator S_0S_i as it occurs in spin factor
-            #print("point_expr", point_expr)
-            cor_expr=point_expr
-            for translation in translation_group:
-                #print("translation", translation)
-                cor_expr+=point_expr.replace_indices({n: n^translation for n in range(9)})
-            #print("cor expr", cor_expr)
-            symmetric_exprs.append(cor_expr)
-    print("sym exprs", symmetric_exprs)
-
-    #Now, let's make a check, which representations will give the same answer. Let's at first focus on a specific value of hamming weight.
-    for basis in hf_bases:
-        print(basis.symmetries)
-        basis.build()
-        cor_op_sym=ls.Operator(symmetric_exprs[0],basis)
-        cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
-        H_op_sym=ls.Operator(H_expr,basis)
-        H_matrix_sym=H_op_sym.to_dense()
-        #print(cor_matrix_sym.shape)
-        #print(H_matrix_sym.shape)
-        corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
-        print("corr", corr)
+        point_expr=ls.Expr("0.0 Sx0 Sx1")
+        for j in orbit:
+            point_expr+=basic_expr.replace_indices({1: j})#make operator S_0S_i as it occurs in spin factor
+        
+        cor_expr=point_expr
+        for translation in translation_group:
+            cor_expr+=point_expr.replace_indices({n: n^translation for n in range(9)})
+
+        norm=1/(len(orbit)*(len(translation_group)+1))#plus 1 is necessary, because in the symmetry array we don't count the identity
+        cor_expr=norm*cor_expr #here we also apply the normalization: divide by number of primitives in symmetric expression
+        symmetric_exprs.append(cor_expr)
+    
+    #Now, let's make a check, that some representations indeed give the same answer. Let's at first focus on a specific value of hamming weight.
+    #for basis in hf_bases:
+    #    print(basis.symmetries)
+    #    basis.build()
+    #    cor_op_sym=ls.Operator(symmetric_exprs[1],basis) # we pick a particular expression
+    #    cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
+    #    H_op_sym=ls.Operator(H_expr,basis)
+    #    H_matrix_sym=H_op_sym.to_dense()
+    #    corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
+    #    print("corr", corr)
 
     #We see that we have only 3 independent values of the correlator. We can do calculate them in a more automatic way.
     #The representations of translation group related to each other by an action of point group (lie in the same conjugacy class) give the same correlators.
-    #Let's select the representatives of the bases lying in the same conjugacy class.
+    #Therefore, we need to divide the bases into different conjugacy classes, and we expect the structure 4x4x1, 
+    #where one quadruple consists of representations where (k_x,k_y) is diagonal, 
+    #another quadruple consists of representations where one of momenta is zero, and the trivial representation
+
     bases_classes=[]
-    #At first we will ort all symmetries, so that we will be able to compare them
-    #for basis in hf_bases:
-    #    basis.symmetries=sorted(basis.symmetries, key=lambda x: list(x[0].array_form))
-    bases_count=hf_bases.copy()         
-    #print("point group", point_group)
-    
-    for basis in hf_basis:
-        #print("basis_symmetry", basis.symmetries) #Print the representation corresponding to the given basis
-        tb=basis.symmetries.copy()
-        conj_class=[]
+    bases_count=hf_bases.copy()          
+    while bases_count!=[]:
+        basis = bases_count[0]
+        #print("#######basis_symmetry#####", basis.symmetries) #Print the representation corresponding to the given basis
+        conj_class=[basis]
+        bases_count.remove(basis)
+        #print("BASES_COUNT", bases_count)
+        del_array=[]
         for g in point_group:
+            tb=basis.symmetries.copy()
+            #print("element of point group", g)
             for i in range(len(tb)):
                 tb[i]=(g*tb[i][0]*(g**(-1)),tb[i][1])
             tb=sorted(tb, key=lambda x: list(x[0].array_form))
-            #print("tb", tb)
+            #print("transformed basis", tb)
             for bs in bases_count:
-                symsort=sorted(bs.symmetries, key=lambda x: list(x[0].array_form))
-                if symsort==tb:
-                    check=True #We check that we do not overcount bases since some elements of the point group can give the same representations
-                    for bs_check in basis_count:
-                        symsort_check=sorted(bs_check.symmetries, key=lambda x: list(x[0].array_form))
-                        if symsort_check==tb:
-                            check=False
-                            #print("check")
-                            continue
-                    if check==True:
-                        conj_class.append(bs)
-                        #print("symsort", symsort)
-                        #print("tb", tb)
-                        #print("symmetries", bs.symmetries)
-        #print(conj_class)
-        if len(conj_class)>0:
-            bases_classes.append(conj_class)
+                if bs.symmetries==tb and (bs not in conj_class):
+                    #print("symmetries are equal")
+                    conj_class.append(bs)
+                    #print("bs_symmetries", bs.symmetries)
+                    del_array.append(bs)
+                    #print("bases_count", bases_count)
+        for del_bs in del_array:
+            bases_count.remove(del_bs)
+        bases_classes.append(conj_class)
     
+
+    #print("BASES CLASSES", bases_classes) #print conjugacy classes, so that we see that we have the correct  
     #Now, let's check if these conjugacy classes indeed give the same results
-    for base_class in bases_classes:
-        print("base_class", base_class)
-        for bs in base_class:
-            print(bs.symmetries)
+    #for base_class in bases_classes:
+    #    print("base_class", base_class)
+    #    for bs in base_class:
+    #        print(bs.symmetries)
+    #        bs.build()
+    #        cor_op_sym=ls.Operator(symmetric_exprs[1],bs)
+    #        cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
+    #        H_op_sym=ls.Operator(H_expr,bs)
+    #        H_matrix_sym=H_op_sym.to_dense()
+    #        corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
+    #        print("corr", corr)
+    
+    #It works! We still need to make a few more steps. 
+    #First, we need to obtain the correct multiplier factor to reproduce the result for the full basis with given magnetization.
+    #Next, we need to sum over all magnetizations to obtain the final result (of course, we need to add proper phases as well, but we did it before)
+    
+    #Let's proceed. We can expect that the weight of a correlator calculated in a given symmetric basis with fixed magnetization
+    #should be divided by total number of primitive expression in the symmetrized expression and multipled by the number of representations in the conjugacy lass.
+    #We also need to take into account the normalization due to the trace calculation which for high temperatures gives (dimension of subspace)/(dimension of the total space)
+    #For low temperatures, the normalilzation should take into account exponents of Hamiltonians (exp(-H_{subspace}T))/(exp(-H_{total}T))
+    
+    #Correlators in spin basis with magnetization 5 are
+    basis_m5=ls.SpinBasis(number_spins=9,hamming_weight=5)
+    basis_m5.build()
+    H_op_m5=ls.Operator(H_expr,basis_m5)
+    H_matrix_m5=H_op_m5.to_dense()#Make numpy matrix for matrix exponentiation
+    cor_array_m5=[]
+    cor_array_m5_se=[]
+    for i in range(N_sites):
+        cor_expr=basic_expr.replace_indices({1: i}) #make operator S_0S_i
+        cor_op_m5=ls.Operator(cor_expr,basis_m5)
+        cor_matrix_m5=cor_op_m5.to_dense() #make numpy matrix
+        cor_array_m5.append(thermal_average(cor_matrix_m5,H_matrix_m5,T)) #calculate observable and add to array
+    print("corrs magnetisation=5", cor_array_m5)
+    
+    #Correlators calculated from representatives
+    sym_corr_m5=[]
+    for sym_expr in symmetric_exprs:
+        sym_corr=0
+        for base_class in bases_classes:
+            bs=base_class[0]
             bs.build()
-            cor_op_sym=ls.Operator(symmetric_exprs[0],bs)
+            cor_op_sym=ls.Operator(sym_expr,bs)
             cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
             H_op_sym=ls.Operator(H_expr,bs)
             H_matrix_sym=H_op_sym.to_dense()
-            #print(cor_matrix_sym.shape)
-            #print(H_matrix_sym.shape)
             corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
-            print("corr", corr)
+            #sym_corr+=(len(base_class)*(bs.number_states)/basis_m5.number_states)*corr #for high temperatures
+            sym_corr+=(len(base_class)*weighted_trace(H_matrix_sym,T)/weighted_trace(H_matrix_m5,T))*corr #general expression
+        sym_corr_m5.append(sym_corr)
+    print("correlators calculated using translation symmetry, magnetisation=5", sym_corr_m5)
+
+    #It works! Now let's write the general cycle.
+    for hamming_weight in range(10):
+
+
+
+
+
 
 
 

From ecf7f69f014e0c17e25cdb237c98dc9f81c98756 Mon Sep 17 00:00:00 2001
From: Askar Iliasov <66430227+asjosik1991@users.noreply.github.com>
Date: Wed, 23 Oct 2024 16:14:44 +0200
Subject: [PATCH 67/67] Update structure_factor.py

---
 python/example/structure_factor.py | 179 +++++++++++++++++------------
 1 file changed, 107 insertions(+), 72 deletions(-)

diff --git a/python/example/structure_factor.py b/python/example/structure_factor.py
index cd7d0b7..7654f80 100644
--- a/python/example/structure_factor.py
+++ b/python/example/structure_factor.py
@@ -39,7 +39,6 @@ def structure_factor_sym(qx,qy,sym_cor_array,coords, orbit_array, L,N_sites):
             sf += np.exp(1j*(qx*rx+qy*ry))*sym_cor_array[orb_ind]/N_sites #here we also use that the coordinates of the site with index 0 are (0,0)
     return np.real(sf)
 
-
 def conjugate_classes(bases, point_group):
     bases_classes=[]
     bases_count=bases.copy()          
@@ -60,7 +59,14 @@ def conjugate_classes(bases, point_group):
         for del_bs in del_array:
             bases_count.remove(del_bs)
         bases_classes.append(conj_class)
-        return bases_classes
+    return bases_classes
+
+def pick_hamming_weight(sectors,n):
+    picked_bases=[] #array of bases with specific hamming weight
+    for basis in sectors:
+        if basis.hamming_weight==n:
+            picked_bases.append(basis)
+    return picked_bases
 
 def main():
     
@@ -154,15 +160,16 @@ def main():
     #At first we will make a summation over different translation sectors using the translations and the action of point group.
     #The idea is to construct the observables invariant under the action of the point group so that there is the exact correspondence between different sectors of translation group.
     #We start with specific hamming weight and generate symmetric bases
-    hf_bases=[] #array of bases with specific hamming weight
-    translation_bases=H_expr.hilbert_space_sectors()
-    for basis in translation_bases:
-        if basis.hamming_weight==5:
-            hf_bases.append(basis)
-    translation_group=H_expr.abelian_permutation_group()
+    #hf_bases=[] #array of bases with specific hamming weight
+    #translation_bases=H_expr.hilbert_space_sectors()
+    #for basis in translation_bases:
+    #    if basis.hamming_weight==5:
+    #        hf_bases.append(basis)
     
     #Make symmetrized expressions for the correlators
     #The symmetric observables are made from already calculated orbits of C_{0,i}
+    
+    translation_group=H_expr.abelian_permutation_group()
     symmetric_exprs=[]
     for orbit in orbit_array:
         point_expr=ls.Expr("0.0 Sx0 Sx1")
@@ -194,32 +201,32 @@ def main():
     #where one quadruple consists of representations where (k_x,k_y) is diagonal, 
     #another quadruple consists of representations where one of momenta is zero, and the trivial representation
 
-    bases_classes=[]
-    bases_count=hf_bases.copy()          
-    while bases_count!=[]:
-        basis = bases_count[0]
-        #print("#######basis_symmetry#####", basis.symmetries) #Print the representation corresponding to the given basis
-        conj_class=[basis]
-        bases_count.remove(basis)
-        #print("BASES_COUNT", bases_count)
-        del_array=[]
-        for g in point_group:
-            tb=basis.symmetries.copy()
-            #print("element of point group", g)
-            for i in range(len(tb)):
-                tb[i]=(g*tb[i][0]*(g**(-1)),tb[i][1])
-            tb=sorted(tb, key=lambda x: list(x[0].array_form))
-            #print("transformed basis", tb)
-            for bs in bases_count:
-                if bs.symmetries==tb and (bs not in conj_class):
-                    #print("symmetries are equal")
-                    conj_class.append(bs)
-                    #print("bs_symmetries", bs.symmetries)
-                    del_array.append(bs)
-                    #print("bases_count", bases_count)
-        for del_bs in del_array:
-            bases_count.remove(del_bs)
-        bases_classes.append(conj_class)
+    #bases_classes=[]
+    #bases_count=hf_bases.copy()          
+    #while bases_count!=[]:
+    #    basis = bases_count[0]
+    #    #print("#######basis_symmetry#####", basis.symmetries) #Print the representation corresponding to the given basis
+    #    conj_class=[basis]
+    #    bases_count.remove(basis)
+    #    #print("BASES_COUNT", bases_count)
+    #    del_array=[]
+    #    for g in point_group:
+    #        tb=basis.symmetries.copy()
+    #        #print("element of point group", g)
+    #        for i in range(len(tb)):
+    #            tb[i]=(g*tb[i][0]*(g**(-1)),tb[i][1])
+    #        tb=sorted(tb, key=lambda x: list(x[0].array_form))
+    #        #print("transformed basis", tb)
+    #        for bs in bases_count:
+    #            if bs.symmetries==tb and (bs not in conj_class):
+    #                #print("symmetries are equal")
+    #                conj_class.append(bs)
+    #                #print("bs_symmetries", bs.symmetries)
+    #                del_array.append(bs)
+    #                #print("bases_count", bases_count)
+    #    for del_bs in del_array:
+    #        bases_count.remove(del_bs)
+    #    bases_classes.append(conj_class)
     
 
     #print("BASES CLASSES", bases_classes) #print conjugacy classes, so that we see that we have the correct  
@@ -246,49 +253,77 @@ def main():
     #For low temperatures, the normalilzation should take into account exponents of Hamiltonians (exp(-H_{subspace}T))/(exp(-H_{total}T))
     
     #Correlators in spin basis with magnetization 5 are
-    basis_m5=ls.SpinBasis(number_spins=9,hamming_weight=5)
-    basis_m5.build()
-    H_op_m5=ls.Operator(H_expr,basis_m5)
-    H_matrix_m5=H_op_m5.to_dense()#Make numpy matrix for matrix exponentiation
-    cor_array_m5=[]
-    cor_array_m5_se=[]
-    for i in range(N_sites):
-        cor_expr=basic_expr.replace_indices({1: i}) #make operator S_0S_i
-        cor_op_m5=ls.Operator(cor_expr,basis_m5)
-        cor_matrix_m5=cor_op_m5.to_dense() #make numpy matrix
-        cor_array_m5.append(thermal_average(cor_matrix_m5,H_matrix_m5,T)) #calculate observable and add to array
-    print("corrs magnetisation=5", cor_array_m5)
+    #basis_m5=ls.SpinBasis(number_spins=9,hamming_weight=5)
+    #basis_m5.build()
+    #H_op_m5=ls.Operator(H_expr,basis_m5)
+    #H_matrix_m5=H_op_m5.to_dense()#Make numpy matrix for matrix exponentiation
+    #cor_array_m5=[]
+    #cor_array_m5_se=[]
+    #for i in range(N_sites):
+    #    cor_expr=basic_expr.replace_indices({1: i}) #make operator S_0S_i
+    #    cor_op_m5=ls.Operator(cor_expr,basis_m5)
+    #    cor_matrix_m5=cor_op_m5.to_dense() #make numpy matrix
+    #    cor_array_m5.append(thermal_average(cor_matrix_m5,H_matrix_m5,T)) #calculate observable and add to array
+    #print("corrs magnetisation=5", cor_array_m5)
     
     #Correlators calculated from representatives
-    sym_corr_m5=[]
-    for sym_expr in symmetric_exprs:
-        sym_corr=0
-        for base_class in bases_classes:
-            bs=base_class[0]
-            bs.build()
-            cor_op_sym=ls.Operator(sym_expr,bs)
-            cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
-            H_op_sym=ls.Operator(H_expr,bs)
-            H_matrix_sym=H_op_sym.to_dense()
-            corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
-            #sym_corr+=(len(base_class)*(bs.number_states)/basis_m5.number_states)*corr #for high temperatures
-            sym_corr+=(len(base_class)*weighted_trace(H_matrix_sym,T)/weighted_trace(H_matrix_m5,T))*corr #general expression
-        sym_corr_m5.append(sym_corr)
-    print("correlators calculated using translation symmetry, magnetisation=5", sym_corr_m5)
+    #sym_corr_m5=[]
+    #for sym_expr in symmetric_exprs:
+    #    sym_corr=0
+    #    for base_class in bases_classes:
+    #        bs=base_class[0]
+    #        bs.build()
+    #        cor_op_sym=ls.Operator(sym_expr,bs)
+    #        cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
+    #        H_op_sym=ls.Operator(H_expr,bs)
+    #        H_matrix_sym=H_op_sym.to_dense()
+    #        corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
+    #        #sym_corr+=(len(base_class)*(bs.number_states)/basis_m5.number_states)*corr #for high temperatures
+    #        sym_corr+=(len(base_class)*weighted_trace(H_matrix_sym,T)/weighted_trace(H_matrix_m5,T))*corr #general expression
+    #    sym_corr_m5.append(sym_corr)
+    #print("correlators calculated using translation symmetry, magnetisation=5", sym_corr_m5)
 
     #It works! Now let's write the general cycle.
-    for hamming_weight in range(10):
-
-
-
-
-
-
-
-
-
-
     
+    
+    translation_bases=H_expr.hilbert_space_sectors()
+    sym_corr_full=[]
+    for sym_expr in symmetric_exprs:
+        print("the correlator for this expression is being calculated", sym_expr)
+        sym_corr=0
+        full_trace=weighted_trace(H_matrix,T)
+        for hamming_weight in range(10):
+            print("hamming weight", hamming_weight)
+            hf_bases=pick_hamming_weight(translation_bases, hamming_weight)
+            #print("bases", hf_bases)
+            bases_classes=conjugate_classes(hf_bases, point_group)
+            #print("bases classes", bases_classes)
+            for base_class in bases_classes:
+                bs=base_class[0]
+                bs.build()
+                cor_op_sym=ls.Operator(sym_expr,bs)
+                cor_matrix_sym=cor_op_sym.to_dense() #make numpy matrix
+                H_op_sym=ls.Operator(H_expr,bs)
+                H_matrix_sym=H_op_sym.to_dense()
+                corr=thermal_average(cor_matrix_sym,H_matrix_sym,T) #calculate observable
+                #sym_corr+=(len(base_class)*(bs.number_states)/basis_m5.number_states)*corr #for high temperatures
+                sym_corr+=(len(base_class)*weighted_trace(H_matrix_sym,T)/full_trace)*corr #general expression
+        sym_corr_full.append(sym_corr)
+    print("correlators calculated from all hilbert sectors", sym_corr_full)
 
+    #Everything works!
+    #Let's plot some figures!
+    SF_sym_full=np.zeros((Nq,Nq))
+    for i in range(Nq):
+        for j in range(Nq):
+            SF_sym_full[i,j]=structure_factor_sym(qxs[i],qys[j],sym_corr_full,coords, orbit_array, L,N_sites)
+    np.save("sym_corr_full.npy", SF_sym_full)
+    
+    #plot structure factor
+    #h = plt.contourf(x, y, SF_sym_full)
+    #plt.axis('scaled')
+    #plt.colorbar()
+    #plt.savefig("structure_factor_using_translations.png")
+    #plt.close()
 
 main()