Repository: tenpy/tenpy (language: python)

Issue:
npc.tensordot drops labels when using axes=0
While playing around I ran into some unexpected (for me at least) behavior for npc.tensordot() when called with the axes=0 option - from what I understand it should be the same as npc.outer(), but somehow when doing it eg. with two simple vectors, the axes labels of the first get lost along the way, so for example:

```
import tenpy.linalg.np_conserved as npc

v1 = npc.Array.from_ndarray_trivial([1,2,3],labels=['one'])
v2 = npc.Array.from_ndarray_trivial([2,3,4],labels=['two'])

tdot = npc.tensordot(v1,v2,axes=0)

tdot
```

will give 
```
<npc.Array shape=(3, 3) labels=['two']>
```

I had a quick look and a possible fix might be as simple as pulling the `elif axes=0` check up in the tensordot function, so that outer() gets called before other possible messups occur inbetween, something along the lines of https://github.com/tenpy/tenpy/commit/af86845bb3cea9f06ba050a5f9cdb5437874c068

For some simple cases it seems to be working fine, I didn't check whether this breaks anything else though.

[Zheng, Y.-G., et al. (October 2022). Efficiently Extracting Multi-Point Correlations of a Floquet Thermalized System. ]
I've used TeNPy in a recent work:  https://arxiv.org/abs/2210.08556  

Please include it into the list of [papers using TeNPy](https://tenpy.readthedocs.io/en/latest/papers_using_tenpy.html)



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/doc-build.yml b/.github/workflows/doc-build.yml
index 501e62ae7..da71cf759 100644
--- a/.github/workflows/doc-build.yml
+++ b/.github/workflows/doc-build.yml
@@ -1,6 +1,11 @@
 name: sphinx-doctest
 
-on: [push]
+on:
+  push: 
+    branches:
+      - main
+  pull_request:
+
 
 jobs:
   build:
diff --git a/doc/conf.py b/doc/conf.py
index 95b8cfd1f..f9a53ea49 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 #
 import sys
 import os
diff --git a/doc/literature.bib b/doc/literature.bib
index bef8c0416..941c81cb8 100644
--- a/doc/literature.bib
+++ b/doc/literature.bib
@@ -1,4 +1,3 @@
-
 @article{affleck1987,
   title = {Rigorous Results on Valence-Bond Ground States in Antiferromagnets},
   author = {Affleck, Ian and Kennedy, Tom and Lieb, Elliott H. and Tasaki, Hal},
@@ -659,5 +658,3 @@ @article{zauner2015
   doi = {10.1088/0953-8984/27/42/425602},
   archiveprefix = {arXiv}
 }
-
-
diff --git a/doc/papers_using_tenpy.bib b/doc/papers_using_tenpy.bib
index a015aab6c..8d1cce0cf 100644
--- a/doc/papers_using_tenpy.bib
+++ b/doc/papers_using_tenpy.bib
@@ -1,4 +1,3 @@
-
 @article{alsallom2021,
   title = {Fate of the Non-{{Hermitian}} Skin Effect in Many-Body Fermionic Systems},
   author = {Alsallom, Faisal and Herviou, Lo{\"i}c and Yazyev, Oleg V. and Brzezi{\'n}ska, Marta},
@@ -67,7 +66,7 @@ @article{andrews2021
   archiveprefix = {arXiv}
 }
 
-@misc{argüello-luengo2022,
+@misc{arguello-luengo2022,
   title = {Tuning Long-Range Fermion-Mediated Interactions in Cold-Atom Quantum Simulators},
   author = {{Arg{\"u}ello-Luengo}, Javier and {Gonz{\'a}lez-Tudela}, Alejandro and {Gonz{\'a}lez-Cuadra}, Daniel},
   year = {2022},
@@ -96,6 +95,20 @@ @article{banuls2020
   archiveprefix = {arXiv}
 }
 
+@misc{bauer2023,
+  title = {Non-{{Abelian Anyons}} with {{Rydberg Atoms}}},
+  author = {Bauer, Nora M. and Kokkas, Elias and Ale, Victor and Siopsis, George},
+  year = {2023},
+  month = jan,
+  number = {arXiv:2301.04250},
+  eprint = {2301.04250},
+  eprinttype = {arxiv},
+  primaryclass = {quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2301.04250},
+  archiveprefix = {arXiv}
+}
+
 @misc{bhakuni2022,
   title = {Interaction-Driven Phase Transition in One Dimensional Mirror-Symmetry Protected Topological Insulator},
   author = {Bhakuni, Devendra Singh and Ghosh, Amrita and Grosfeld, Eytan},
@@ -310,6 +323,20 @@ @article{cervia2022
   archiveprefix = {arXiv}
 }
 
+@misc{chan2022,
+  title = {Plaquette Valence Bond State in Spin-1/2 {{J1-J2 XY}} Model on Square Lattice},
+  author = {Chan, Y.-H. and Jiang, Hong-Chen and Chen, Y.-C.},
+  year = {2022},
+  month = dec,
+  number = {arXiv:2212.04863},
+  eprint = {2212.04863},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2212.04863},
+  archiveprefix = {arXiv}
+}
+
 @misc{chen2022,
   title = {Continuous {{Symmetry Breaking}} in a {{Two-dimensional Rydberg Array}}},
   author = {Chen, Cheng and Bornet, Guillaume and Bintz, Marcus and Emperauger, Gabriel and Leclerc, Lucas and Liu, Vincent S. and Scholl, Pascal and Barredo, Daniel and Hauschild, Johannes and Chatterjee, Shubhayu and Schuler, Michael and Laeuchli, Andreas M. and Zaletel, Michael P. and Lahaye, Thierry and Yao, Norman Y. and Browaeys, Antoine},
@@ -393,6 +420,37 @@ @misc{cookmeyer2022
   archiveprefix = {arXiv}
 }
 
+@article{cubela2023,
+  title = {Particle Zoo in a Doped Spin Chain: {{Correlated}} States of Mesons and Magnons},
+  shorttitle = {Particle Zoo in a Doped Spin Chain},
+  author = {{\v C}ubela, Petar and Bohrdt, Annabelle and Greiner, Markus and Grusdt, Fabian},
+  year = {2023},
+  month = jan,
+  journal = {Physical Review B},
+  volume = {107},
+  number = {3},
+  eprint = {2210.02320},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  pages = {035105},
+  doi = {10.1103/PhysRevB.107.035105},
+  archiveprefix = {arXiv}
+}
+
+@misc{dai2022,
+  title = {Fermionic {{Isometric Tensor Network States}} in {{Two Dimensions}}},
+  author = {Dai, Zhehao and Wu, Yantao and Wang, Taige and Zaletel, Michael P.},
+  year = {2022},
+  month = dec,
+  number = {arXiv:2211.00043},
+  eprint = {2211.00043},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2211.00043},
+  archiveprefix = {arXiv}
+}
+
 @article{demidio2021,
   title = {Diagnosing Weakly First-Order Phase Transitions by Coupling to Order Parameters},
   author = {D'Emidio, Jonathan and Eberharter, Alexander A. and L{\"a}uchli, Andreas M.},
@@ -420,6 +478,20 @@ @article{denner2020
   archiveprefix = {arXiv}
 }
 
+@misc{desaules2023,
+  title = {Ergodicity {{Breaking Under Confinement}} in {{Cold-Atom Quantum Simulators}}},
+  author = {Desaules, Jean-Yves and Su, Guo-Xian and McCulloch, Ian P. and Yang, Bing and Papi{\'c}, Zlatko and Halimeh, Jad C.},
+  year = {2023},
+  month = jan,
+  number = {arXiv:2301.07717},
+  eprint = {2301.07717},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:hep-lat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2301.07717},
+  archiveprefix = {arXiv}
+}
+
 @article{devakul2021,
   ids = {devakul2021a},
   title = {Quantum Anomalous {{Hall}} Effect from Inverted Charge Transfer Gap},
@@ -505,6 +577,20 @@ @article{dong2020
   archiveprefix = {arXiv}
 }
 
+@misc{drescher2022,
+  title = {Dynamical {{Signatures}} of {{Symmetry Broken}} and {{Liquid Phases}} in an \${{S}}=1/2\$ {{Heisenberg Antiferromagnet}} on the {{Triangular Lattice}}},
+  author = {Drescher, Markus and Vanderstraeten, Laurens and Moessner, Roderich and Pollmann, Frank},
+  year = {2022},
+  month = sep,
+  number = {arXiv:2209.03344},
+  eprint = {2209.03344},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2209.03344},
+  archiveprefix = {arXiv}
+}
+
 @misc{fromholz2022,
   title = {Phase Diagram of {{Rydberg-dressed}} Atoms on Two-Leg Triangular Ladders},
   author = {Fromholz, Pierre and Tsitsishvili, Mikheil and Votto, Matteo and Dalmonte, Marcello and Nersesyan, Alexander and Chanda, Titas},
@@ -602,6 +688,20 @@ @article{hartmann2021
   archiveprefix = {arXiv}
 }
 
+@misc{hatsugai2022,
+  title = {Topological Pump of \${{SU}}({{Q}})\$ Quantum Chain and {{Diophantine}} Equation},
+  author = {Hatsugai, Yasuhiro and Kuno, Yoshihito},
+  year = {2022},
+  month = oct,
+  number = {arXiv:2210.11646},
+  eprint = {2210.11646},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:hep-th},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2210.11646},
+  archiveprefix = {arXiv}
+}
+
 @article{hauschild2018_dup,
   title = {Finding Purifications with Minimal Entanglement},
   author = {Hauschild, Johannes and Leviatan, Eyal and Bardarson, Jens H. and Altman, Ehud and Zaletel, Michael P. and Pollmann, Frank},
@@ -808,6 +908,20 @@ @article{jin2021
   archiveprefix = {arXiv}
 }
 
+@misc{jin2023,
+  title = {Bosonization of the Interacting {{Su-Schrieffer-Heeger}} Model},
+  author = {Jin, Tony and Ruggiero, Paola and Giamarchi, Thierry},
+  year = {2023},
+  month = jan,
+  number = {arXiv:2301.00472},
+  eprint = {2301.00472},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2301.00472},
+  archiveprefix = {arXiv}
+}
+
 @misc{jobst2022,
   title = {Finite-Depth Scaling of Infinite Quantum Circuits for Quantum Critical Points},
   author = {Jobst, Bernhard and Smith, Adam and Pollmann, Frank},
@@ -849,6 +963,20 @@ @article{julia-farre2021
   archiveprefix = {arXiv}
 }
 
+@misc{julia-farre2023,
+  title = {Topological Stripe State in an Extended {{Fermi-Hubbard}} Model},
+  author = {{Juli{\`a}-Farr{\'e}}, Sergi and Cardarelli, Lorenzo and Lewenstein, Maciej and M{\"u}ller, Markus and Dauphin, Alexandre},
+  year = {2023},
+  month = jan,
+  number = {arXiv:2301.03312},
+  eprint = {2301.03312},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2301.03312},
+  archiveprefix = {arXiv}
+}
+
 @misc{kadow2022,
   title = {Hole {{Spectral Function}} of a {{Chiral Spin Liquid}} in the {{Triangular Lattice Hubbard Model}}},
   author = {Kadow, Wilhelm and Vanderstraeten, Laurens and Knap, Michael},
@@ -863,7 +991,7 @@ @misc{kadow2022
   archiveprefix = {arXiv}
 }
 
-@misc{kebrič2022,
+@misc{kebric2022,
   title = {Confinement {{Induced Frustration}} in a {{One-Dimensional Z}}\_2 {{Lattice Gauge Theory}}},
   author = {Kebri{\v c}, Matja{\v z} and Borla, Umberto and Schollw{\"o}ck, Ulrich and Moroz, Sergej and Barbiero, Luca and Grusdt, Fabian},
   year = {2022},
@@ -978,6 +1106,20 @@ @article{kottmann2021
   archiveprefix = {arXiv}
 }
 
+@misc{kottmann2022,
+  title = {Investigating {{Quantum Many-Body Systems}} with {{Tensor Networks}}, {{Machine Learning}} and {{Quantum Computers}}},
+  author = {Kottmann, Korbinian},
+  year = {2022},
+  month = oct,
+  number = {arXiv:2210.11130},
+  eprint = {2210.11130},
+  eprinttype = {arxiv},
+  primaryclass = {quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2210.11130},
+  archiveprefix = {arXiv}
+}
+
 @article{kraus2021,
   title = {Quantum Phases of Dipolar Bosons in One-Dimensional Optical Lattices},
   author = {Kraus, Rebecca and Chanda, Titas and Zakrzewski, Jakub and Morigi, Giovanna},
@@ -990,6 +1132,34 @@ @article{kraus2021
   archiveprefix = {arXiv}
 }
 
+@misc{kuhlenkamp2022,
+  title = {Tunable Topological Order of Pseudo Spins in Semiconductor Heterostructures},
+  author = {Kuhlenkamp, Clemens and Kadow, Wilhelm and Imamoglu, Atac and Knap, Michael},
+  year = {2022},
+  month = sep,
+  number = {arXiv:2209.05506},
+  eprint = {2209.05506},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2209.05506},
+  archiveprefix = {arXiv}
+}
+
+@misc{kumar2022,
+  title = {Thermal {{Hall}} Conductivity near Field-Suppressed Magnetic Order in a {{Kitaev-Heisenberg}} Model},
+  author = {Kumar, Aman and Tripathi, Vikram},
+  year = {2022},
+  month = oct,
+  number = {arXiv:2208.13708},
+  eprint = {2208.13708},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2208.13708},
+  archiveprefix = {arXiv}
+}
+
 @article{kuno2020,
   title = {Interaction {{Induced Topological Charge Pump}}},
   author = {Kuno, Yoshihito and Hatsugai, Yasuhiro},
@@ -1098,6 +1268,48 @@ @article{li2021
   archiveprefix = {arXiv}
 }
 
+@misc{liu2022,
+  title = {Model-{{Independent Learning}} of {{Quantum Phases}} of {{Matter}} with {{Quantum Convolutional Neural Networks}}},
+  author = {Liu, Yu-Jie and Smith, Adam and Knap, Michael and Pollmann, Frank},
+  year = {2022},
+  month = dec,
+  number = {arXiv:2211.11786},
+  eprint = {2211.11786},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2211.11786},
+  archiveprefix = {arXiv}
+}
+
+@misc{luo2022,
+  title = {Global Quantum Phase Diagram and Non-{{Abelian}} Chiral Spin Liquid in a Spin-3/2 Square Lattice Antiferromagnet},
+  author = {Luo, Wei-Wei and Huang, Yixuan and Sheng, D. N. and Zhu, W.},
+  year = {2022},
+  month = dec,
+  number = {arXiv:2212.14223},
+  eprint = {2212.14223},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2212.14223},
+  archiveprefix = {arXiv}
+}
+
+@misc{ma2022,
+  title = {Average {{Symmetry-Protected Topological Phases}}},
+  author = {Ma, Ruochen and Wang, Chong},
+  year = {2022},
+  month = sep,
+  number = {arXiv:2209.02723},
+  eprint = {2209.02723},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2209.02723},
+  archiveprefix = {arXiv}
+}
+
 @article{mahyaeh2020,
   title = {Phase Diagram of the {{Z}}\_3-{{Fock}} Parafermion Chain with Pair Hopping},
   author = {Mahyaeh, Iman and Wouters, Jurriaan and Schuricht, Dirk},
@@ -1110,6 +1322,21 @@ @article{mahyaeh2020
   archiveprefix = {arXiv}
 }
 
+@misc{mendez-cordoba2022,
+  title = {From Edge to Bulk: {{Cavity}} Induced Displacement of Topological Non-Local Qubits},
+  shorttitle = {From Edge to Bulk},
+  author = {{M{\'e}ndez-C{\'o}rdoba}, F. P. M. and Rodr{\'i}guez, F. J. and Tejedor, C. and Quiroga, L.},
+  year = {2022},
+  month = dec,
+  number = {arXiv:2211.14145},
+  eprint = {2211.14145},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2211.14145},
+  archiveprefix = {arXiv}
+}
+
 @misc{metz2022,
   title = {Self-{{Correcting Quantum Many-Body Control}} Using {{Reinforcement Learning}} with {{Tensor Networks}}},
   author = {Metz, Friederike and Bukov, Marin},
@@ -1163,6 +1390,20 @@ @article{morera2021
   archiveprefix = {arXiv}
 }
 
+@misc{morera2022,
+  title = {High-Temperature Kinetic Magnetism in Triangular Lattices},
+  author = {Morera, Ivan and {Kan{\'a}sz-Nagy}, M{\'a}rton and Smolenski, Tomasz and Ciorciaro, Livio and Imamo{\u g}lu, Ata{\c c} and Demler, Eugene},
+  year = {2022},
+  month = sep,
+  number = {arXiv:2209.05398},
+  eprint = {2209.05398},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2209.05398},
+  archiveprefix = {arXiv}
+}
+
 @article{motruk2020,
   title = {Detecting Fractional {{Chern}} Insulators in Optical Lattices through Quantized Displacement},
   author = {Motruk, Johannes and Na, Ilyoun},
@@ -1175,6 +1416,20 @@ @article{motruk2020
   archiveprefix = {arXiv}
 }
 
+@misc{motruk2022,
+  title = {Kagome {{Chiral Spin Liquid}} in {{Transition Metal Dichalcogenide Moir}}\textbackslash '\{e\} {{Bilayers}}},
+  author = {Motruk, Johannes and Rossi, Dario and Abanin, Dmitry A. and Rademaker, Louk},
+  year = {2022},
+  month = nov,
+  number = {arXiv:2211.15696},
+  eprint = {2211.15696},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2211.15696},
+  archiveprefix = {arXiv}
+}
+
 @misc{nietner2022,
   title = {A Route towards Engineering Many-Body Localization in Real Materials},
   author = {Nietner, A. and Kshetrimayum, A. and Eisert, J. and Lake, B.},
@@ -1257,6 +1512,20 @@ @article{parker2021a
   archiveprefix = {arXiv}
 }
 
+@misc{passetti2022,
+  title = {Cavity {{Light-Matter Entanglement}} through {{Quantum Fluctuations}}},
+  author = {Passetti, Giacomo and Eckhardt, Christian J. and Sentef, Michael A. and Kennes, Dante M.},
+  year = {2022},
+  month = dec,
+  number = {arXiv:2212.03011},
+  eprint = {2212.03011},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:physics, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2212.03011},
+  archiveprefix = {arXiv}
+}
+
 @article{pena2021,
   title = {Fractional Resonances and Prethermal States in {{Floquet}} Systems},
   author = {Pe{\~n}a, R. and Bastidas, V. M. and Torres, F. and Munro, W. J. and Romero, G.},
@@ -1269,6 +1538,36 @@ @article{pena2021
   archiveprefix = {arXiv}
 }
 
+@article{pena2022,
+  title = {Stable Many-Body Resonances in Open Quantum Systems},
+  author = {Pe{\~n}a, Ruben and Kyaw, Thi Ha and Romero, Guillermo},
+  year = {2022},
+  month = dec,
+  journal = {Symmetry},
+  volume = {14},
+  number = {12},
+  eprint = {2209.07307},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  pages = {2562},
+  doi = {10.3390/sym14122562},
+  archiveprefix = {arXiv}
+}
+
+@misc{prembabu2022,
+  title = {Boundary Deconfined Quantum Criticality at Transitions between Symmetry-Protected Topological Chains},
+  author = {Prembabu, Saranesh and Thorngren, Ryan and Verresen, Ruben},
+  year = {2022},
+  month = aug,
+  number = {arXiv:2208.12258},
+  eprint = {2208.12258},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2208.12258},
+  archiveprefix = {arXiv}
+}
+
 @article{quella2020,
   title = {Symmetry Protected Topological Phases beyond Groups: {{The}} q-Deformed Bilinear-Biquadratic Spin Chain},
   shorttitle = {Symmetry Protected Topological Phases beyond Groups},
@@ -1453,6 +1752,34 @@ @article{semeghini2021
   archiveprefix = {arXiv}
 }
 
+@misc{sim2022,
+  title = {Kekule Spin-Orbit Dimer Phase and Triplon Dynamics},
+  author = {Sim, GiBaik},
+  year = {2022},
+  month = nov,
+  number = {arXiv:2211.16094},
+  eprint = {2211.16094},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2211.16094},
+  archiveprefix = {arXiv}
+}
+
+@misc{sim2022a,
+  title = {Nonlinear Spectroscopy of Bound States in Perturbed {{Ising}} Spin Chains},
+  author = {Sim, GiBaik and Knolle, Johannes and Pollmann, Frank},
+  year = {2022},
+  month = sep,
+  number = {arXiv:2209.00720},
+  eprint = {2209.00720},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2209.00720},
+  archiveprefix = {arXiv}
+}
+
 @article{soejima2020,
   title = {Efficient Simulation of Moire Materials Using the Density Matrix Renormalization Group},
   author = {Soejima, Tomohiro and Parker, Daniel E. and Bultinck, Nick and Hauschild, Johannes and Zaletel, Michael P.},
@@ -1531,7 +1858,7 @@ @article{strkalj2021
   archiveprefix = {arXiv}
 }
 
-@misc{štrkalj2022,
+@misc{strkalj2022,
   title = {Coexistence of Localization and Transport in Many-Body Two-Dimensional {{Aubry-Andr}}\textbackslash 'e Models},
   author = {{\v S}trkalj, Antonio and Doggen, Elmer V. H. and Castelnovo, Claudio},
   year = {2022},
@@ -1946,6 +2273,20 @@ @misc{wybo2022
   archiveprefix = {arXiv}
 }
 
+@misc{xie2022,
+  title = {Complete Field-Induced Spectral Response of the Spin-1/2 Triangular-Lattice Antiferromagnet {{CsYbSe}}\$\_2\$},
+  author = {Xie, Tao and Eberharter, A. A. and Xing, Jie and Nishimoto, S. and Brando, M. and Khanenko, P. and Sichelschmidt, J. and Turrini, A. A. and Mazzone, D. G. and Naumov, P. G. and Sanjeewa, L. D. and Harrison, N. and Sefat, A. S. and Normand, B. and Lauchli, A. M. and Podlesnyak, A. and Nikitin, S. E.},
+  year = {2022},
+  month = oct,
+  number = {arXiv:2210.04928},
+  eprint = {2210.04928},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2210.04928},
+  archiveprefix = {arXiv}
+}
+
 @article{you2021,
   title = {Fractonic Critical Point Proximate to a Higher-Order Topological Insulator: {{How}} Does {{UV}} Blend with {{IR}}?},
   shorttitle = {Fractonic Critical Point Proximate to a Higher-Order Topological Insulator},
@@ -1973,6 +2314,20 @@ @misc{zechmann2022
   archiveprefix = {arXiv}
 }
 
+@misc{zechmann2022a,
+  title = {Fractonic {{Luttinger Liquids}} and {{Supersolids}} in a {{Constrained Bose-Hubbard Model}}},
+  author = {Zechmann, Philip and Altman, Ehud and Knap, Michael and Feldmeier, Johannes},
+  year = {2022},
+  month = oct,
+  number = {arXiv:2210.11072},
+  eprint = {2210.11072},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2210.11072},
+  archiveprefix = {arXiv}
+}
+
 @article{zhang2020,
   title = {Symmetric Pseudogap Metal in a Generalized T-{{J}} Model},
   author = {Zhang, Ya-Hui and Zhu, Zheng},
@@ -2012,6 +2367,21 @@ @misc{zhang2022
   archiveprefix = {arXiv}
 }
 
+@misc{zhang2022a,
+  title = {Transformer {{Quantum State}}: {{A Multi-Purpose Model}} for {{Quantum Many-Body Problems}}},
+  shorttitle = {Transformer {{Quantum State}}},
+  author = {Zhang, Yuan-Hang and Di Ventra, Massimiliano},
+  year = {2022},
+  month = dec,
+  number = {arXiv:2208.01758},
+  eprint = {2208.01758},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:physics, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2208.01758},
+  archiveprefix = {arXiv}
+}
+
 @misc{zhao2022,
   title = {The Fractional {{Chern}} Insulator with {{Rydberg-dressed}} Neutral Atoms},
   author = {Zhao, Yang and Shi, Xiao-Feng},
@@ -2026,6 +2396,20 @@ @misc{zhao2022
   archiveprefix = {arXiv}
 }
 
+@misc{zheng2022,
+  title = {Efficiently {{Extracting Multi-Point Correlations}} of a {{Floquet Thermalized System}}},
+  author = {Zheng, Yong-Guang and Zhang, Wei-Yong and Shen, Ying-Chao and Luo, An and Liu, Ying and He, Ming-Gen and Zhang, Hao-Ran and Lin, Wan and Wang, Han-Yi and Zhu, Zi-Hang and Chen, Ming-Cheng and Lu, Chao-Yang and Thanasilp, Supanut and Angelakis, Dimitris G. and Yuan, Zhen-Sheng and Pan, Jian-Wei},
+  year = {2022},
+  month = oct,
+  number = {arXiv:2210.08556},
+  eprint = {2210.08556},
+  eprinttype = {arxiv},
+  primaryclass = {cond-mat, physics:quant-ph},
+  publisher = {{arXiv}},
+  doi = {10.48550/arXiv.2210.08556},
+  archiveprefix = {arXiv}
+}
+
 @misc{zhou2022,
   title = {Observation of {{Universal Hall Response}} in {{Strongly Interacting Fermions}}},
   author = {Zhou, T.-W. and Cappellini, G. and Tusi, D. and Franchi, L. and Parravicini, J. and Repellin, C. and Greschner, S. and Inguscio, M. and Giamarchi, T. and Filippone, M. and Catani, J. and Fallani, L.},
@@ -2039,5 +2423,3 @@ @misc{zhou2022
   doi = {10.48550/arXiv.2205.13567},
   archiveprefix = {arXiv}
 }
-
-
diff --git a/examples/a_np_conserved.py b/examples/a_np_conserved.py
index fcf4e8d31..e9783ade8 100644
--- a/examples/a_np_conserved.py
+++ b/examples/a_np_conserved.py
@@ -13,7 +13,7 @@
 Compare it to the example `b_mps.py`,
 which does the same steps using a few predefined classes like MPS and MPO.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy.linalg.np_conserved as npc
 import numpy as np
diff --git a/examples/advanced/central_charge_ising.py b/examples/advanced/central_charge_ising.py
index e645dd4ab..0e632830d 100644
--- a/examples/advanced/central_charge_ising.py
+++ b/examples/advanced/central_charge_ising.py
@@ -6,7 +6,7 @@
 
 For the theoretical background why :math:`S = c/6 log(xi)`, see :cite:`pollmann2009`.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import tenpy
diff --git a/examples/advanced/mpo_exponential_decay.py b/examples/advanced/mpo_exponential_decay.py
index 8307a4b97..bdead5827 100644
--- a/examples/advanced/mpo_exponential_decay.py
+++ b/examples/advanced/mpo_exponential_decay.py
@@ -10,7 +10,7 @@
 We run the iDMRG algorithm to find the ground state and energy density of the
 system in the thermodynamic limit.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import tenpy.linalg.np_conserved as npc
diff --git a/examples/advanced/tfi_phase_transition.py b/examples/advanced/tfi_phase_transition.py
index 5c7ace45f..893862b37 100644
--- a/examples/advanced/tfi_phase_transition.py
+++ b/examples/advanced/tfi_phase_transition.py
@@ -4,7 +4,7 @@
 through the phase transition by changing the field `g`.
 It plots a few observables in the end.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/examples/advanced/xxz_corr_length.py b/examples/advanced/xxz_corr_length.py
index 4955f6835..a978f6ace 100644
--- a/examples/advanced/xxz_corr_length.py
+++ b/examples/advanced/xxz_corr_length.py
@@ -5,7 +5,7 @@
 :meth:`~tenpy.networks.mps.MPS.correlation_length` to extract the correlation length of the ground
 state, and plots it vs. hz in the end.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/examples/b_mps.py b/examples/b_mps.py
index d8f899fb0..143e9199d 100644
--- a/examples/b_mps.py
+++ b/examples/b_mps.py
@@ -12,7 +12,7 @@
 Note that this example performs the same steps as `a_np_conserved.py`,
 but makes use of other predefined classes except npc.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy.linalg.np_conserved as npc
 import numpy as np
diff --git a/examples/c_tebd.py b/examples/c_tebd.py
index 859698509..f23fca454 100644
--- a/examples/c_tebd.py
+++ b/examples/c_tebd.py
@@ -3,7 +3,7 @@
 The example functions in this class do the same as the ones in `toycodes/c_tebd.py`, but make use
 of the classes defined in tenpy.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/examples/chern_insulators/chiral_pi_flux.py b/examples/chern_insulators/chiral_pi_flux.py
index 55f4dc56b..8833b493e 100644
--- a/examples/chern_insulators/chiral_pi_flux.py
+++ b/examples/chern_insulators/chiral_pi_flux.py
@@ -3,7 +3,7 @@
 Based on the model in [Neupert2011]_
 """
 
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/examples/chern_insulators/haldane.py b/examples/chern_insulators/haldane.py
index 448859869..5bbcf9c32 100644
--- a/examples/chern_insulators/haldane.py
+++ b/examples/chern_insulators/haldane.py
@@ -3,7 +3,7 @@
 Reproduces Fig. 2.a,b) in [Grushin2015]_
 """
 
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/examples/chern_insulators/haldane_C3.py b/examples/chern_insulators/haldane_C3.py
index 9025250f5..219766cac 100644
--- a/examples/chern_insulators/haldane_C3.py
+++ b/examples/chern_insulators/haldane_C3.py
@@ -3,7 +3,7 @@
 Based on the model in [Yang2012]_
 """
 
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/examples/chern_insulators/haldane_FCI.py b/examples/chern_insulators/haldane_FCI.py
index fba86cd50..6d885622f 100644
--- a/examples/chern_insulators/haldane_FCI.py
+++ b/examples/chern_insulators/haldane_FCI.py
@@ -6,7 +6,7 @@
 - 1/2 filling of the lowest band (i.e. 1/4 total filling)
 """
 
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import warnings
diff --git a/examples/d_dmrg.py b/examples/d_dmrg.py
index d60749436..ecaf162f4 100644
--- a/examples/d_dmrg.py
+++ b/examples/d_dmrg.py
@@ -3,7 +3,7 @@
 The example functions in this class do the same as the ones in `toycodes/d_dmrg.py`,
 but make use of the classes defined in tenpy.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/examples/e_tdvp.py b/examples/e_tdvp.py
index b6af46c0e..242f30bc1 100644
--- a/examples/e_tdvp.py
+++ b/examples/e_tdvp.py
@@ -4,7 +4,7 @@
 difference is that we can run one-site TDVP or two-site TDVP. In the former, the bond dimension can
 not grow; the latter allows to grow the bond dimension and hence requires a truncation.
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 import numpy as np
 import tenpy.linalg.np_conserved as npc
 import tenpy
diff --git a/examples/tfi_exact.py b/examples/tfi_exact.py
index 622ac3ff0..30fc14d2b 100644
--- a/examples/tfi_exact.py
+++ b/examples/tfi_exact.py
@@ -4,7 +4,7 @@
 .. math ::
     H = - J \\sum_{i} \\sigma^x_i \\sigma^x_{i+1} - g \\sum_{i} \\sigma^z_i
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import scipy.sparse as sparse
diff --git a/examples/z_exact_diag.py b/examples/z_exact_diag.py
index bfcce5233..d094c6746 100644
--- a/examples/z_exact_diag.py
+++ b/examples/z_exact_diag.py
@@ -2,7 +2,7 @@
 
 Sorry that this is not well documented! ED is meant to be used for debugging only ;)
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy.linalg.np_conserved as npc
 from tenpy.models.xxz_chain import XXZChain
diff --git a/setup.py b/setup.py
index 6a0a12052..97fa9f1e3 100644
--- a/setup.py
+++ b/setup.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 from setuptools import setup, find_packages
 from setuptools import Extension
 import numpy
@@ -77,7 +77,7 @@ def write_version_py(full_version, git_rev, filename='tenpy/_version.py'):
     content = """\
 # THIS FILE IS GENERATED FROM setup.py
 # thus, it contains the version during compilation
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 version = '{version!s}'
 short_version = 'v' + version
 released = {released!s}
diff --git a/tenpy/__init__.py b/tenpy/__init__.py
index 89374ac21..2e8fcff96 100644
--- a/tenpy/__init__.py
+++ b/tenpy/__init__.py
@@ -6,7 +6,7 @@
 The code is intended to be accessible for newcommers
 and yet powerful enough for day-to-day research.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 # This file marks this directory as a python package.
 
 import warnings
diff --git a/tenpy/__main__.py b/tenpy/__main__.py
index e1b1f21ab..00b713934 100644
--- a/tenpy/__main__.py
+++ b/tenpy/__main__.py
@@ -1,5 +1,5 @@
 """The tenpy entry point, called by ``python -m tenpy``."""
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy
 
diff --git a/tenpy/algorithms/__init__.py b/tenpy/algorithms/__init__.py
index df22f731f..1fc326355 100644
--- a/tenpy/algorithms/__init__.py
+++ b/tenpy/algorithms/__init__.py
@@ -18,7 +18,7 @@
     exact_diag
     disentangler
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from . import algorithm, truncation, dmrg, dmrg_parallel, disentangler, mps_common, tebd, tdvp, \
     exact_diag, purification, network_contractor, mpo_evolution
diff --git a/tenpy/algorithms/algorithm.py b/tenpy/algorithms/algorithm.py
index 7c11f7ea2..c421be504 100644
--- a/tenpy/algorithms/algorithm.py
+++ b/tenpy/algorithms/algorithm.py
@@ -1,5 +1,5 @@
 """This module contains some base classes for algorithms."""
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import warnings
 import time
diff --git a/tenpy/algorithms/disentangler.py b/tenpy/algorithms/disentangler.py
index 250762445..71a065d40 100644
--- a/tenpy/algorithms/disentangler.py
+++ b/tenpy/algorithms/disentangler.py
@@ -7,7 +7,7 @@
 
 .. autodata:: disentanglers_atom_parse_dict
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import logging
diff --git a/tenpy/algorithms/dmrg.py b/tenpy/algorithms/dmrg.py
index 4b7b64228..922651659 100644
--- a/tenpy/algorithms/dmrg.py
+++ b/tenpy/algorithms/dmrg.py
@@ -30,7 +30,7 @@
 A generic protocol for approaching a physics question using DMRG is given in
 :doc:`/intro/dmrg-protocol`.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import time
diff --git a/tenpy/algorithms/dmrg_parallel.py b/tenpy/algorithms/dmrg_parallel.py
index 019628fb8..f7a0a1e43 100644
--- a/tenpy/algorithms/dmrg_parallel.py
+++ b/tenpy/algorithms/dmrg_parallel.py
@@ -3,7 +3,7 @@
 .. warning ::
     This module is still under active development. Use with care!
 """
-# Copyright 2021 TeNPy Developers, GNU GPLv3
+# Copyright 2021-2023 TeNPy Developers, GNU GPLv3
 
 from ..tools.thread import Worker
 
diff --git a/tenpy/algorithms/exact_diag.py b/tenpy/algorithms/exact_diag.py
index 7f6ccc13f..82bcc353d 100644
--- a/tenpy/algorithms/exact_diag.py
+++ b/tenpy/algorithms/exact_diag.py
@@ -13,7 +13,7 @@
     but just the ability to diagonalize the defined models for small system sizes
     without addional extra work.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import warnings
diff --git a/tenpy/algorithms/mpo_evolution.py b/tenpy/algorithms/mpo_evolution.py
index 6058c4484..ab89190b5 100644
--- a/tenpy/algorithms/mpo_evolution.py
+++ b/tenpy/algorithms/mpo_evolution.py
@@ -1,6 +1,6 @@
 """Time evolution using the WI or WII approximation of the time evolution operator."""
 
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import time
diff --git a/tenpy/algorithms/mps_common.py b/tenpy/algorithms/mps_common.py
index d1d7dafa0..f5ad5bc3d 100644
--- a/tenpy/algorithms/mps_common.py
+++ b/tenpy/algorithms/mps_common.py
@@ -18,7 +18,7 @@
 The :class:`VariationalCompression` and :class:`VariationalApplyMPO`
 implemented here also directly use the :class:`Sweep` class.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from .algorithm import Algorithm
 from ..linalg.sparse import NpcLinearOperator, SumNpcLinearOperator, OrthogonalNpcLinearOperator
diff --git a/tenpy/algorithms/mps_sweeps.py b/tenpy/algorithms/mps_sweeps.py
index b1c5e1bb8..8053de421 100644
--- a/tenpy/algorithms/mps_sweeps.py
+++ b/tenpy/algorithms/mps_sweeps.py
@@ -5,7 +5,7 @@
 .. deprecated :: 0.7.0
     This module is just around for backwards compatibility.
 """
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import warnings
 
diff --git a/tenpy/algorithms/network_contractor.py b/tenpy/algorithms/network_contractor.py
index d710c39b5..0ccd59d47 100644
--- a/tenpy/algorithms/network_contractor.py
+++ b/tenpy/algorithms/network_contractor.py
@@ -13,7 +13,7 @@
     - improve helpfulness of Warnings
     - _do_trace: trace over all pairs of legs at once. need the corresponding npc function first.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import warnings
@@ -25,6 +25,46 @@
 outer_product = -66666666  #: a constant that represents an outer product in the sequence of ncon
 
 
+def ncon(tensor_list, leg_links, sequence=None):
+    """Implementation of ``ncon.m`` for TeNPy Arrays.
+
+    This function is a python implementation of ``ncon.m`` (:arxiv:`1304.6112`) for tenpy
+    :class:`~tenpy.linalg.np_conserved.Array`.
+    :func:`contract` is a wrapper that translates from a more python/tenpy input style
+
+    Parameters
+    ----------
+    tensor_list : iterable of :class:'Array'
+        Tensors to be contracted.
+    leg_links : iterable of iterable of int
+        Each entry of leg_links describes the connectivity of the corresponding tensor in
+        `tensor_list`.
+        Each entry is a sequence (e.g. a list) that has an integer value for each leg of the corresponding tensor.
+        Positive values ``1,2,...`` are labels of contracted legs and should appear
+        exactly twice in `leg_links`.
+        Negative values ``-1,-2,-3,...`` are labels of uncontracted legs and indicate the final ordering
+        (``-1`` is the first axis).
+    sequence : list of int, optional
+        The order in which the contractions (indicated by positive values in `leg_links` are to be performed.
+        Ascending order is used by default
+        An entry of network_contractor.outer_product indicates performing an outer product.
+        This corresponds to the zero-in-sequence convention of :arxiv:`1304.6112`
+
+    Returns
+    -------
+    result : :class:`Array` | complex
+        The number or tensor resulting from the contraction.
+    """
+    tensor_list, leg_links, sequence = _ncon_input_checks(tensor_list, leg_links, sequence)
+    tensor_list, leg_links, sequence = _ncon_do_traces(tensor_list, leg_links, sequence)
+    tensor_list, leg_links, sequence = _ncon_do_binary_contractions(tensor_list, leg_links, sequence)
+    tensor_list, leg_links = _ncon_do_outer_products(tensor_list, leg_links)
+    result, = tensor_list
+    if len(leg_links[0]) > 0:
+        result.itranspose(np.argsort(-leg_links[0]))
+    return result
+
+
 def contract(tensor_list, tensor_names=None, leg_contractions=None, open_legs=None, sequence=None):
     """Contract a network of tensors.
 
@@ -74,7 +114,7 @@ def contract(tensor_list, tensor_names=None, leg_contractions=None, open_legs=No
 
     # default sequence
     if sequence is None:
-        sequence = list(range(len(leg_contractions)))
+        sequence = list(range(1, len(leg_contractions) + 1))
 
     # translate leg_contractions and open_legs to a leg_links list as used by ncon
     # initialise leg_links
@@ -83,14 +123,14 @@ def contract(tensor_list, tensor_names=None, leg_contractions=None, open_legs=No
         leg_links.append([None] * len(tensor.legs))
 
     # fill in the contractions
-    contraction_counter = 0
+    contraction_counter = 1
     new_sequence = []
     for n in sequence:
         if n == outer_product:
             new_sequence.append(outer_product)
             continue
 
-        con = leg_contractions[n]
+        con = leg_contractions[n - 1]
         leg_idx1 = tensor_list[con[0]].get_leg_index(con[1])
         leg_idx2 = tensor_list[con[2]].get_leg_index(con[3])
 
@@ -121,395 +161,139 @@ def contract(tensor_list, tensor_names=None, leg_contractions=None, open_legs=No
         res.iset_leg_labels(final_labels)
 
     return res
+    
 
-
-def ncon(tensor_list, leg_links, sequence):
-    """Implementation of ``ncon.m`` for TeNPy Arrays.
-
-    This function is a python implementation of ``ncon.m`` (:arxiv:`1304.6112`) for tenpy
-    :class:`~tenpy.linalg.np_conserved.Array`.
-    :func:`contract` is a wrapper that translates from a more python/tenpy input style
-
-    Parameters
-    ----------
-    tensor_list : list of :class:'Array'
-        Tensors to be contracted.
-    leg_links : list of list of int
-        Each entry of leg_links describes the connectivity of the corresponding tensor in
-        `tensor_list`.
-        Each entry is a list that has an entry for each leg of the corresponding tensor.
-        Values ``0,1,2,...`` are labels of contracted legs and should appear
-        exactly twice in `leg_links`.
-        Values ``-1,-2,-3,...`` are labels of uncontracted legs and indicate the final ordering
-        (``-1`` is first axis).
-    sequence : list of int
-        The order in which the contractions are to be performed.
-        An entry of network_contractor.outer_product indicates performing an outer product.
-        This corresponds to the zero-in-sequence convention of :arxiv:`1304.6112`
-
-    Returns
-    -------
-    result : :class:`Array` | complex
-        The number or tensor resulting from the contraction.
+def _ncon_input_checks(tensor_list, leg_links, sequence):
+    """Check inputs for consistency and convert to the following format
+    leg_links: list of np.ndarray
+    sequence: np.ndarray
     """
-
     tensor_list = list(tensor_list)
-
-    # check contractibility of legs?
-
-    def flatten(l):
-        return [item for sublist in l for item in sublist]
-
-    while len(tensor_list) > 1 or any(i >= 0 for i in flatten(leg_links)):
-        if len(sequence) == 0:
-            sequence = [outer_product] * (len(tensor_list) - 1)
-        if sequence[0] == outer_product:
-            # outer product
-            tensor_list, leg_links, sequence = _outer_product(tensor_list, leg_links, sequence)
-        else:
-            # identify and perform tensor contraction
-            # find tensors that the index sequence[0] corresponds to
-            tensors = []
-            for a in range(len(leg_links)):
-                if sequence[0] in leg_links[a]:
-                    tensors.append(a)
-
-            if len(tensors) == 1:
-                # its a trace
-
-                # find all traced indices on that tensor
-                traced_indices = np.sort(leg_links[tensors[0]])
-                traced_indices = [
-                    item for item, count in collections.Counter(traced_indices).items()
-                    if count == 2
-                ]
-
-                # check if this is in accordance with sequence and update sequence
-                doing_traces, sequence = _find_in_sequence(traced_indices, sequence)
-                if not np.array_equal(np.sort(traced_indices), np.sort(doing_traces)):
-                    warnings.warn(
-                        'Suboptimal contraction sequence. When tracing legs ' + str(doing_traces) +
-                        ' the legs ' +
-                        str(list(filter(lambda n: n not in doing_traces, traced_indices))) +
-                        'should also be traced')
-                # TODO translate this back to human readable leg label?
-
-                # perform traces
-                tensor_list[tensors[0]], leg_links[tensors[0]] = \
-                    _do_trace(tensor_list[tensors[0]], leg_links[tensors[0]], doing_traces)
-
-                # update leg links
-                for idx in doing_traces:
-                    leg_links[tensors[0]] = list(filter(lambda b: b != idx, leg_links[tensors[0]]))
-
-            else:
-                # its a contraction
-
-                # find all other contracted legs between the two tensors
-                common_indices = []
-                for idx in leg_links[tensors[0]]:
-                    if idx in leg_links[tensors[1]]:
-                        common_indices.append(idx)
-
-                # check if this is in accordance with sequence
-                contraction_indices, sequence = _find_in_sequence(common_indices, sequence)
-                if not np.array_equal(np.sort(contraction_indices), np.sort(common_indices)):
-                    warnings.warn(
-                        'Suboptimal contraction sequence. When contracting legs ' +
-                        str(contraction_indices) + ' the legs ' +
-                        str(list(filter(lambda n: n not in contraction_indices, common_indices))) +
-                        'should also be traced')
-                    # TODO translate this back to human readable leg names
-
-                # are there traced indices on either of these tensors?
-                # noinspection PyArgumentList
-                traces0 = [
-                    item for item, count in collections.Counter(leg_links[tensors[0]]).items()
-                    if count == 2
-                ]
-                traces1 = [
-                    item for item, count in collections.Counter(leg_links[tensors[1]]).items()
-                    if count == 2
-                ]
-                if len(traces0) + len(traces1) != 0:
-                    warnings.warn(
-                        'Suboptimal contraction sequence. When processing ' + str(None) +
-                        ' one of the involved tensors have legs that could be traced over first')
-                    # TODO human readable identifier
-                # contract all contraction_indices and update leg_links
-                tensor_list[tensors[0]], leg_links[tensors[0]] = \
-                    _tcontract(tensor_list[tensors[0]], tensor_list[tensors[1]], leg_links[tensors[0]],
-                               leg_links[tensors[1]], contraction_indices)
-                tensor_list.pop(tensors[1])
-                leg_links.pop(tensors[1])
-    assert len(tensor_list) == 1
-    return tensor_list[0]
-
-
-def _find_in_sequence(indices, sequence):
-    """Helper function for ncon.
-
-    check if the supplied indices appear at the beginning of sequence
-
-    Parameters
-    ----------
-    indices : list of int
-    sequence : list of int
-
-    Returns
-    -------
-    idcs : list of int
-        All the given indices that appear consecutively at the front of sequence
-    """
-
-    ptr = 0
-    while ptr < len(sequence) and sequence[ptr] in indices:
-        ptr = ptr + 1
-    rtn_indices = sequence[:ptr]
-    sequence = sequence[ptr:]
-    return rtn_indices, sequence
+    
+    if len(leg_links) != len(tensor_list):
+        msg = f'Mismatching lengths: Got {len(tensor_list)} Arrays and {len(leg_links)} leg_links'
+        raise ValueError(msg)
+    leg_links = [np.asarray(tensor_links) for tensor_links in leg_links]
+    leg_links_flat = [value for tensor_links in leg_links for value in tensor_links]
+    num_contractions = max((value for value in leg_links_flat if value > 0), default=0)
+    missing_positive_values = [n for n in range(1, num_contractions + 1) if n not in leg_links_flat]
+    if missing_positive_values:
+        which = ", ".join(map(str, missing_positive_values))
+        msg = f'The following positive values are missing in leg_links: {which}'
+        raise ValueError(msg)
+    
+    if sequence is None:
+        sequence = np.arange(1, num_contractions + 1)
+    else:
+        sequence = np.asarray(sequence)
+        if len(sequence) != num_contractions or set(sequence) != set(range(1, num_contractions + 1)):
+            msg = f'Invalid sequence. Expected a permutation of [1, ..., {num_contractions}]. Got {sequence}.'
+            raise ValueError(msg)
+        if outer_product in sequence:
+            raise ValueError('Outer product in sequence is deprecated.')
+        
+    return tensor_list, leg_links, sequence
 
 
-def _do_trace(a, leg_link, traced_indices):
-    """Helper function for ncon.
+def _ncon_do_traces(tensor_list, leg_links, sequence):
+    for n in range(len(leg_links)):
+        if len(leg_links[n]) > len(np.unique(leg_links[n])):
+            tensor_list[n], leg_links[n], used_values = _partial_trace(tensor_list[n], leg_links[n], loc=n)
+            sequence = np.delete(sequence, np.intersect1d(sequence, used_values, return_indices=True)[1])
+    return tensor_list, leg_links, sequence
 
-    Trace over pair(s) of legs on a given tensor
-    Update the leg_link entry
-    .. todo :
-        perform all traces simultaneously
 
+def _partial_trace(tensor, tensor_links, loc):
+    """Perform all partial traces on a given tensor.
+    
     Parameters
     ----------
-    a : :class:'Array'
-        the tensor to perform traces on
-    leg_link : list of int
-        the leg_links entry of a
-    traced_indices list of int
-        the labels of the legs to be traced. each should appear twice in leg_link
-
+    tensor: :class:'Array'
+    tensor_links: np.ndarray
+        the corresponding entry of `tensor_links`
+    loc: int
+        the index of `tensor` in `tensor_list`
+        
     Returns
     -------
-    a : :class:'Array
+    result: :class:'Array'
         the traced tensor
-    leg_link : list of int
-        the updated entry for leg_links
-    """
-
-    untraced_indices = list(filter(lambda n: n not in traced_indices, leg_link))
-
-    # sort traced legs into two blocks blocks to trace over
-    block_a = []
-    block_b = []
-    block_untraced = []
-    for idx in traced_indices:
-        # position of first and last appearance of idx in leg_labels
-        block_a.append(leg_link.index(idx))
-        block_b.append(len(leg_link) - 1 - leg_link[::-1].index(idx))
-    for idx in untraced_indices:
-        block_untraced.append(leg_link.index(idx))
-
-    new_order = block_untraced + [l[n] for n in range(len(block_a)) for l in [block_a, block_b]]
-    a.itranspose(new_order)
-    for n in range(len(new_order) - 1, len(block_untraced), -2):
-        a = npc.trace(a, n, n - 1)
-    leg_link = untraced_indices
-
-    return a, leg_link
-
-
-def _tcontract(t1, t2, links1, links2, contract_legs):
-    """Helper function for ncon.
-
-    Contract two tensors along one or multiple axis
-
-    Parameters
-    ----------
-    t1 : :class:'Array'
-        The first tensor.
-    t2 : :class:'Array'
-        The second tensor.
-    links1 : list of int
-        The leg_links entry of the first tensor
-    links2 : list of int
-        The leg_links entry of the second tensor
-    contract_legs : list of int
-        The labels of the legs to be contracted. Each should appear exactly once in links1 and exactly once in links2.
-
-    Returns
-    -------
-    res : :class:'Array'
-        The result of the pairwise contraction
-    links : list of int
-        a leg_links entry for the res tensor
+    tensor_links: np.ndarray
+        the proper replacement entry for `tensor_links`
+    used_values: np.ndarray
+        the entries of `tensor_links` that indicated the trace and are missing in the returned `tensor_links`
     """
-
-    # if either tensor is not an :class:'Array' try converting
-    # this may occur if a closed disconnected diagram is part of the contraction.
-    # ncon will then try to process the resulting number with _tcontract of _outer_product
-    if type(t1) is not npc.Array:
-        t1 = npc.Array.from_ndarray_trivial(t1)
-    if type(t2) is not npc.Array:
-        t2 = npc.Array.from_ndarray_trivial(t2)
-
-    # find uncontracted legs
-    free_legs_1 = list(filter(lambda n: n not in contract_legs, links1))
-    free_legs_2 = list(filter(lambda n: n not in contract_legs, links2))
-
-    # find positions of legs
-    pos_cont_legs_1 = [links1.index(leg) for leg in contract_legs]
-    pos_cont_legs_2 = [links2.index(leg) for leg in contract_legs]
-
-    # contract
-    res = npc.tensordot(t1, t2, axes=(pos_cont_legs_1, pos_cont_legs_2))
-    # tensordot keeps order of uncontracted legs intact. first those of T1 then those of T2
-    links = free_legs_1 + free_legs_2
-
-    return res, links
-
-
-def _outer_product(tensor_list, leg_links, sequence):
-    """Helper function for ncon.
-
-    Perform an outer product of multiple tensors and optionally contract all their legs with one single tensor.
-    This can be caused a value OP in the sequence or if there are more than one tensors
-    left but no legs to be contracted
-
-    Details see :arxiv:`1304.6112`
-
-    Parameters
-    ----------
-    tensor_list : list of :class:'Array'
-        the whole list of tensors currently processed by ncon
-    leg_links : list of list of int
-        the whole leg_links currently processed by ncon
-    sequence :
-     the remaining sequence currently processed by ncon
-
-    Returns
-    -------
-    tensor_list : list of :class:'Array'
-        an updated list of tensors containing the result of the outer product and all untouched tensors
-    leg_links : list of list of int
-        the corresponding updated leg_links
-    sequence : list of int
-        the remaining sequence
-    """
-
-    if all(n == outer_product for n in sequence):
-        # final outer product of all remaining tensors
-        # ensure there are enough entries
-        if len(sequence) < len(tensor_list) - 1:
-            sequence = [outer_product] * len(tensor_list - 1)
-            warnings.warn('Not enough OP entries in sequence')
-
-    # determine number of pending outer products
-    num_op = len(
-        sequence
-    )  # default value: if no entry in sequence is find that is not outer_product, all of them are
-    for n in range(len(sequence)):
-        if sequence[n] != outer_product:
-            num_op = n
-            break
-
-    # determine tensors on which OPs are to be performed
-    # for num_op outer products we need num_op+1 tensors that are all contracted with one extra tensor
-    # see :arxiv:`1304.6112`
-    # find the next num_op+2 tensors coming up in sequence
-    # failure to find this many implies an invalid sequence.
-
-    if num_op == len(leg_links) - 1:
-        # OP of all remaining tensors
-        op_list = list(range(len(leg_links)))
-    else:
-        # flag relevant num_op + 2 tensors
-        flags = [False] * len(leg_links)
-        ptr = num_op  # sequence[num_op] is the first entry that is not _outer
-        while sum(flags) < num_op + 2:
-            if ptr >= len(sequence):
-                raise ValueError(
-                    'sequence contains outer products but ended before finding enough '
-                    'tensors for the outer product.')
-            if sequence[ptr] == outer_product:
-                raise ValueError(
-                    'sequence contains outer products but ncon encountered another OP before '
-                    'identifying all tensors involved in the first.')
-            count = 0
-            for a in range(len(leg_links)):
-                if sequence[ptr] in leg_links[a]:
-                    flags[a] = True
-                    count = count + 1
-            if count != 2:
-                raise ValueError(
-                    'sequence contains outer products. An index on one of the legs is not appearing twice.'
-                )
-            ptr = ptr + 1
-
-        # identify which of these tensors is *not* participating in the OP
-        # but is instead contracted with the result of the OP
-        # - identify the two tensors involved in the first contraction
-        # - examine following contractions until one of them is with a third tensor
-        #   (thus only occurs on one of the initial candidates)
-        first_tensors = [None, None]
-        ptr = num_op
-        for a in range(len(leg_links)):
-            if sequence[ptr] in leg_links[a]:
-                if first_tensors[0] is None:
-                    first_tensors[0] = a
-                else:
-                    first_tensors[1] = a
-                    break
-        done = False
-        next_tensors = [None, None]
-        while not done:
-            next_tensors = [None, None]
-            ptr = ptr + 1
-            for a in range(len(leg_links)):
-                if sequence[ptr] in leg_links[a]:
-                    if next_tensors[0] is None:
-                        next_tensors[0] = a
-                    else:
-                        next_tensors[1] = a
-                        break
-            if not first_tensors == next_tensors:
-                done = True
-        if next_tensors[0] in first_tensors:
-            post_op_tensor = next_tensors[0]
-        else:
-            post_op_tensor = next_tensors[1]
-        flags[post_op_tensor] = False
-        op_list = list(filter(lambda n: flags[n], range(len(leg_links))))
-
-        # check that all indices of op_list are contracted with post_op_tensor
-        op_indices = [leg_links[n] for n in op_list]
-        op_indices = [item for sublist in op_indices for item in sublist]
-        for a in range(len(op_indices)):
-            if not op_indices[a] in leg_links[post_op_tensor]:
-                raise ValueError(
-                    'Outer product failure. OP tensor has contraction with multiple other tensors.')
-
-    # if either tensor is not an :class:'Array' try converting
-    # this may occur if a closed disconnected diagram is part of the contraction.
-    # ncon will then try to process the resulting number with _tcontract of _outer_product
-    for n in op_list:
-        if type(tensor_list[n]) is not npc.Array:
-            tensor_list[n] = npc.Array.from_ndarray_trivial(tensor_list[n])
-
-    # perform OPs, starting with the smallest tensors
-    op_sizes = [tensor_list[n].size for n in op_list]
-    while len(op_sizes) > 1:
-        order = np.argsort(op_sizes)
-
-        # construct outer product of the two smallest tensors
-        tensor_list[op_list[order[0]]] = npc.outer(tensor_list[op_list[order[0]]],
-                                                   tensor_list[op_list[order[1]]])
-        tensor_list.pop(op_list[order[1]])
-        leg_links[op_list[order[0]]] = leg_links[op_list[order[0]]] + leg_links[op_list[order[1]]]
-        leg_links.pop(op_list[order[1]])
-
-        # re adjust op_sizes
-        op_sizes[order[0]] = op_sizes[order[0]] + op_sizes[order[1]]
-        op_sizes.pop(order[1])
-        op_list = [n - 1 if n > op_list[order[1]] else n for n in op_list]
-        op_list.pop(order[1])
-
+    num_ocurrences = np.sum(tensor_links[:, None] == tensor_links[None, :], axis=1)
+    trace_links = np.unique(tensor_links[np.where(num_ocurrences > 1)[0]])
+    res_links = tensor_links[np.where(num_ocurrences == 1)[0]]
+    num_traces = len(trace_links)
+    assert num_traces > 0
+
+    trace_axes = np.zeros((num_traces, 2), dtype=np.int_)
+    for n, value in enumerate(trace_links):
+        trace_axes[n, :] = np.where(tensor_links == value)[0]
+    other_axes = [ax for ax in range(len(tensor_links)) if ax not in trace_axes]
+    
+    try:
+        tensor = tensor.combine_legs(combine_legs=(trace_axes[:, 0], trace_axes[:, 1]), new_axes=(-2, -1))
+        tensor = npc.trace(tensor, leg1=-2, leg2=-1)
+    except Exception as e:
+        msg = (f'An error ocurred while performing the partial trace in tensor_list[{loc}]. '
+               f'Original stacktrace below.')
+        raise type(e)(msg) from e
+    
+    return tensor, res_links, trace_links
+    
+
+def _ncon_do_binary_contractions(tensor_list, leg_links, sequence):
+    while len(sequence) > 0:
+        to_contract = sequence[0]
+        which_tensors = [n for n in range(len(tensor_list)) if sum(leg_links[n] == to_contract) > 0]
+        try:
+            loc_a, loc_b = which_tensors
+        except ValueError:
+            msg = f'Invalid leg_links. Value {to_contract} appeared on {len(which_tensors)} different tensors!'
+            raise ValueError(msg) from None
+
+        # pop b first since lob_b > loc_a
+        tensor_b = tensor_list.pop(loc_b)
+        links_b = leg_links.pop(loc_b)
+        tensor_a = tensor_list.pop(loc_a)
+        links_a = leg_links.pop(loc_a)
+        
+        common_values, a_axes, b_axes = np.intersect1d(links_a, links_b, assume_unique=True, return_indices=True)
+        
+        try:
+            res = npc.tensordot(tensor_a, tensor_b, (a_axes, b_axes))
+        except Excpetion as e:
+            msg = (f'An error occured while performing the pairwise contraction indicated by '
+                   f'values {", ".join(common_values)} in leg_links. '
+                   f'Original stacktrace below.')
+            raise type(e)(msg) from e
+
+        res_links = np.append(np.delete(links_a, a_axes), np.delete(links_b, b_axes))
+        
+        tensor_list.append(res)
+        leg_links.append(res_links)
+        used_sequence_values = np.intersect1d(sequence, common_values, assume_unique=True, return_indices=True)[1]
+        sequence = np.delete(sequence, used_sequence_values)
+        
     return tensor_list, leg_links, sequence
+    
+    
+def _ncon_do_outer_products(tensor_list, leg_links):
+    while len(tensor_list) > 1:
+        tensor_b = tensor_list.pop(-1)
+        links_b = leg_links.pop(-1)
+        try:
+            tensor_list[-1] = npc.outer(tensor_list[-1], tensor_b)
+        except Exception as e:
+            msg = (f'An error occured while performing a final outer product between the last two '
+                   f'of {len(tensor_list) + 1} remaining tensors. '
+                   f'Original stacktrace below.')
+            raise type(e)(msg) from e
+        leg_links[-1] = np.append(leg_links[-1], links_b)
+    
+    return tensor_list, leg_links
+    
+        
\ No newline at end of file
diff --git a/tenpy/algorithms/purification.py b/tenpy/algorithms/purification.py
index dbaa237c7..18035e445 100644
--- a/tenpy/algorithms/purification.py
+++ b/tenpy/algorithms/purification.py
@@ -1,6 +1,6 @@
 """Algorithms for using Purification."""
 
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import logging
diff --git a/tenpy/algorithms/purification_tebd.py b/tenpy/algorithms/purification_tebd.py
index 8a4ad8a69..0b16a96ee 100644
--- a/tenpy/algorithms/purification_tebd.py
+++ b/tenpy/algorithms/purification_tebd.py
@@ -6,7 +6,7 @@
 .. deprecated :: 0.7.1
     This module is just around for backwards compatibility.
 """
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 from .purification import PurificationTEBD, PurificationTEBD2
 from .disentangler import *
diff --git a/tenpy/algorithms/tdvp.py b/tenpy/algorithms/tdvp.py
index a6e3eab20..6a934b688 100644
--- a/tenpy/algorithms/tdvp.py
+++ b/tenpy/algorithms/tdvp.py
@@ -29,7 +29,7 @@
 .. todo ::
     allow for increasing bond dimension in SingleSiteTDVPEngine, similar to DMRG Mixer
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 from tenpy.linalg.lanczos import LanczosEvolution
 from tenpy.algorithms.truncation import svd_theta, TruncationError
diff --git a/tenpy/algorithms/tebd.py b/tenpy/algorithms/tebd.py
index 5706f008b..8cbb05d9d 100644
--- a/tenpy/algorithms/tebd.py
+++ b/tenpy/algorithms/tebd.py
@@ -37,7 +37,7 @@
     Yet, imaginary TEBD might be usefull for cross-checks and testing.
 
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import time
diff --git a/tenpy/algorithms/truncation.py b/tenpy/algorithms/truncation.py
index c7bcf88d0..554ccb402 100644
--- a/tenpy/algorithms/truncation.py
+++ b/tenpy/algorithms/truncation.py
@@ -41,7 +41,7 @@
     There might be other sources of error as well, for example TEBD has also an discretisation
     error depending on the chosen time step.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 from ..linalg import np_conserved as npc
diff --git a/tenpy/linalg/__init__.py b/tenpy/linalg/__init__.py
index 9c7747544..2a58e78f4 100644
--- a/tenpy/linalg/__init__.py
+++ b/tenpy/linalg/__init__.py
@@ -21,7 +21,7 @@
     lanczos
 
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from . import charges, np_conserved, lanczos, random_matrix, sparse, svd_robust
 
diff --git a/tenpy/linalg/_npc_helper.pyx b/tenpy/linalg/_npc_helper.pyx
index e88586500..6b22bf7c2 100644
--- a/tenpy/linalg/_npc_helper.pyx
+++ b/tenpy/linalg/_npc_helper.pyx
@@ -8,7 +8,7 @@ functions/classes defined here to overwrite those written in pure Python wheneve
 decorator ``@use_cython`` is used in other python files of tenpy.
 If this module was not compiled and could not be imported, a warning is issued.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 DEF DEBUG_PRINT = 0  # set this to 1 for debug output (e.g. benchmark timings within the functions)
 DEF USE_MKL_GEMM_BATCH = 1 # whether to use ?gemm_batch function of MKL
diff --git a/tenpy/linalg/charges.py b/tenpy/linalg/charges.py
index d57a487f9..07b7a2403 100644
--- a/tenpy/linalg/charges.py
+++ b/tenpy/linalg/charges.py
@@ -19,7 +19,7 @@
 
 .. autodata:: QTYPE
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import copy
diff --git a/tenpy/linalg/lanczos.py b/tenpy/linalg/lanczos.py
index f38a245c5..84af92f50 100644
--- a/tenpy/linalg/lanczos.py
+++ b/tenpy/linalg/lanczos.py
@@ -1,5 +1,5 @@
 """Lanczos algorithm for np_conserved arrays."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import warnings
 import numpy as np
diff --git a/tenpy/linalg/np_conserved.py b/tenpy/linalg/np_conserved.py
index 887ea0121..a4f4c9c0a 100644
--- a/tenpy/linalg/np_conserved.py
+++ b/tenpy/linalg/np_conserved.py
@@ -81,7 +81,7 @@
     speigs
 
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import scipy.linalg
@@ -3522,7 +3522,7 @@ def tensordot(a, b, axes=2):
         # #### the main work
         res = _tensordot_worker(a, b, axes)
     # labels
-    res._labels = _drop_duplicate_labels(a._labels[:-axes], b._labels[axes:])
+    res._labels = _drop_duplicate_labels(a._labels[:a.rank-axes], b._labels[axes:])
     return res
 
 
diff --git a/tenpy/linalg/random_matrix.py b/tenpy/linalg/random_matrix.py
index c5532d03f..54c1491ae 100644
--- a/tenpy/linalg/random_matrix.py
+++ b/tenpy/linalg/random_matrix.py
@@ -29,7 +29,7 @@
     npc.Array.from_func_square(GOE, [leg, leg.conj()])
 
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/linalg/sparse.py b/tenpy/linalg/sparse.py
index 54e3cb918..5a2cf6b0c 100644
--- a/tenpy/linalg/sparse.py
+++ b/tenpy/linalg/sparse.py
@@ -7,7 +7,7 @@
 :class:`FlatLinearOperator` allows to use all the scipy sparse methods by providing functionality
 to convert flat numpy arrays to and from np_conserved arrays.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 from . import np_conserved as npc
diff --git a/tenpy/linalg/svd_robust.py b/tenpy/linalg/svd_robust.py
index 28fe55420..9e90ffd2c 100644
--- a/tenpy/linalg/svd_robust.py
+++ b/tenpy/linalg/svd_robust.py
@@ -34,7 +34,7 @@
 >>> from tenpy.linalg.svd_robust import svd
 >>> U, S, VT = svd([[1., 1.], [0., 1.]])
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import scipy
diff --git a/tenpy/models/__init__.py b/tenpy/models/__init__.py
index 9c5281e26..1e5b1c496 100644
--- a/tenpy/models/__init__.py
+++ b/tenpy/models/__init__.py
@@ -31,7 +31,7 @@
     toric_code
     mixed_xk
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from . import lattice, model
 from . import tf_ising, xxz_chain, spins, spins_nnn
diff --git a/tenpy/models/aklt.py b/tenpy/models/aklt.py
index 25cae6212..f6934285c 100644
--- a/tenpy/models/aklt.py
+++ b/tenpy/models/aklt.py
@@ -4,7 +4,7 @@
 Writing down the Hamiltonian is easiest done in terms of bond couplings.
 This class thus serves as an example how this can be done.
 """
-# Copyright 2021 TeNPy Developers, GNU GPLv3
+# Copyright 2021-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/models/bose_hubbard.py b/tenpy/models/bose_hubbard.py
index 80d0b7276..1ee019ea5 100644
--- a/tenpy/models/bose_hubbard.py
+++ b/tenpy/models/bose_hubbard.py
@@ -3,7 +3,7 @@
 .. deprecated :: 0.4.1
     This module is just around for backwards compatibility.
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 from .hubbard import BoseHubbardModel, BoseHubbardChain
 
diff --git a/tenpy/models/bose_hubbard_chain.py b/tenpy/models/bose_hubbard_chain.py
index f659ae083..d3db6c697 100644
--- a/tenpy/models/bose_hubbard_chain.py
+++ b/tenpy/models/bose_hubbard_chain.py
@@ -5,7 +5,7 @@
 .. deprecated :: 0.4.1
     This module is just around for backwards compatibility.
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 from .hubbard import BoseHubbardModel, BoseHubbardChain
 
diff --git a/tenpy/models/fermion_chain.py b/tenpy/models/fermion_chain.py
index 37adbeeb2..c765f40d7 100644
--- a/tenpy/models/fermion_chain.py
+++ b/tenpy/models/fermion_chain.py
@@ -3,7 +3,7 @@
 .. deprecated :: 0.4.1
     This module is just around for backwards compatibility.
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 from .fermions_spinless import FermionModel, FermionChain
 
diff --git a/tenpy/models/fermions_hubbard.py b/tenpy/models/fermions_hubbard.py
index 1773d6a29..ffd7840dd 100644
--- a/tenpy/models/fermions_hubbard.py
+++ b/tenpy/models/fermions_hubbard.py
@@ -3,7 +3,7 @@
 .. deprecated :: 0.4.1
     This module is just around for backwards compatibility.
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 from .hubbard import FermiHubbardModel as FermionicHubbardModel
 from .hubbard import FermiHubbardChain as FermionicHubbardChain
diff --git a/tenpy/models/fermions_spinless.py b/tenpy/models/fermions_spinless.py
index 6883b6a9c..083901f43 100644
--- a/tenpy/models/fermions_spinless.py
+++ b/tenpy/models/fermions_spinless.py
@@ -2,7 +2,7 @@
 
 .. todo ::     add further terms (e.g. c^dagger c^dagger + h.c.) to the Hamiltonian.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/models/haldane.py b/tenpy/models/haldane.py
index 2b12b7545..fccf24cef 100644
--- a/tenpy/models/haldane.py
+++ b/tenpy/models/haldane.py
@@ -1,5 +1,5 @@
 """Bosonic and fermionic Haldane models."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/models/hofstadter.py b/tenpy/models/hofstadter.py
index 1513ee73c..5e86c9be9 100644
--- a/tenpy/models/hofstadter.py
+++ b/tenpy/models/hofstadter.py
@@ -4,7 +4,7 @@
     Long term: implement different lattices.
     Long term: implement variable hopping strengths Jx, Jy.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import warnings
diff --git a/tenpy/models/hubbard.py b/tenpy/models/hubbard.py
index 1ca32a6b0..dacde4526 100644
--- a/tenpy/models/hubbard.py
+++ b/tenpy/models/hubbard.py
@@ -1,5 +1,5 @@
 """Bosonic and fermionic Hubbard models."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/models/lattice.py b/tenpy/models/lattice.py
index 1650cd389..40809d497 100644
--- a/tenpy/models/lattice.py
+++ b/tenpy/models/lattice.py
@@ -16,7 +16,7 @@
 :doc:`/notebooks/90_overview_predefined_lattices`
 
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import itertools
diff --git a/tenpy/models/mixed_xk.py b/tenpy/models/mixed_xk.py
index 85671af18..a7f890955 100644
--- a/tenpy/models/mixed_xk.py
+++ b/tenpy/models/mixed_xk.py
@@ -62,7 +62,7 @@
     exclusion principle implies a possibly large occupation on single k modes, i.e., hard-core
     bosons in x-y-space don't map to hard-core bosons in x-k-space!
 """
-# Copyright 2021 TeNPy Developers, GNU GPLv3
+# Copyright 2021-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import itertools as it
diff --git a/tenpy/models/model.py b/tenpy/models/model.py
index d96b3996d..5f672b272 100644
--- a/tenpy/models/model.py
+++ b/tenpy/models/model.py
@@ -27,7 +27,7 @@
 
 See also the introduction in :doc:`/intro/model`.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import warnings
@@ -231,8 +231,8 @@ def update_time_parameter(self, new_time):
         reads out the parameter ``self.options['time']``, and reinitialize/update the model
         calling this method.
 
-        Paramters
-        ---------
+        Parameters
+        ----------
         new_time : float
             Time at which the (time-dependent) Hamiltonian should be constructed.
 
diff --git a/tenpy/models/spins.py b/tenpy/models/spins.py
index 49f568c8f..8ff399ea2 100644
--- a/tenpy/models/spins.py
+++ b/tenpy/models/spins.py
@@ -2,7 +2,7 @@
 
 Uniform lattice of spin-S sites, coupled by nearest-neighbour interactions.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/models/spins_nnn.py b/tenpy/models/spins_nnn.py
index 57b72e38e..790965d2c 100644
--- a/tenpy/models/spins_nnn.py
+++ b/tenpy/models/spins_nnn.py
@@ -14,7 +14,7 @@
 :meth:`~tenpy.models.model.NearestNeighbormodel.from_MPOModel`.
 An example for such a case is given in the file ``examples/c_tebd.py``.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/models/tf_ising.py b/tenpy/models/tf_ising.py
index 5959a28ce..bef31c3eb 100644
--- a/tenpy/models/tf_ising.py
+++ b/tenpy/models/tf_ising.py
@@ -6,7 +6,7 @@
 
 We choose the field along z to allow to conserve the parity, if desired.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/models/toric_code.py b/tenpy/models/toric_code.py
index 65a13d31e..b7ec346be 100644
--- a/tenpy/models/toric_code.py
+++ b/tenpy/models/toric_code.py
@@ -3,7 +3,7 @@
 As we put the model on a cylinder, the name "toric code" is a bit misleading, but it is the
 established name for this model...
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/models/xxz_chain.py b/tenpy/models/xxz_chain.py
index 8d5256d8b..4f1ca59e7 100644
--- a/tenpy/models/xxz_chain.py
+++ b/tenpy/models/xxz_chain.py
@@ -3,7 +3,7 @@
 The XXZ chain is contained in the more general :class:`~tenpy.models.spins.SpinChain`; the idea of
 this module is more to serve as a pedagogical example for a model.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/networks/__init__.py b/tenpy/networks/__init__.py
index b1eb5fae7..388a797d9 100644
--- a/tenpy/networks/__init__.py
+++ b/tenpy/networks/__init__.py
@@ -14,7 +14,7 @@
     terms
     purification_mps
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from . import site, mps, mpo, purification_mps
 
diff --git a/tenpy/networks/mpo.py b/tenpy/networks/mpo.py
index 8924edf8d..56353726d 100644
--- a/tenpy/networks/mpo.py
+++ b/tenpy/networks/mpo.py
@@ -34,7 +34,7 @@
 Similar as for the MPS, a bond index ``i`` is *left* of site `i`,
 i.e. between sites ``i-1`` and ``i``.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 from scipy.linalg import expm
diff --git a/tenpy/networks/mps.py b/tenpy/networks/mps.py
index dc72f08b9..9cd0fe925 100644
--- a/tenpy/networks/mps.py
+++ b/tenpy/networks/mps.py
@@ -144,7 +144,7 @@
 iMPS, so this is the form *assumed* when calling MPO :meth:`~tenpy.networks.mpo.MPO.apply` on an
 MPS.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import warnings
@@ -1433,6 +1433,7 @@ def group_split(self, trunc_par=None):
         """
         if trunc_par is None:
             trunc_par = {}
+        trunc_par = asConfig(trunc_par, 'trunc_params')
         self.convert_form('B')
         if self.L > 1:
             trunc_par.setdefault('chi_max', max(self.chi))
@@ -1468,7 +1469,7 @@ def group_split(self, trunc_par=None):
                 U, S, V, err, _ = svd_theta(theta, trunc_par, inner_labels=['vR', 'vL'])
                 Ss_new.append(S)
                 trunc_err += err
-                theta = U.split_legs(0)  # vL p0 ... pj-1 vR
+                theta = U.scale_axis(S, 'vR').split_legs(0)  # vL p0 ... pj-1 vR
                 for _ in range(n_p_label):
                     combine[0].pop()
                 B = V.split_legs(1).iset_leg_labels(self._B_labels)  # vL p vR
@@ -1667,6 +1668,8 @@ def entanglement_entropy(self, n=1, bonds=None, for_matrix_S=False):
         if bonds is None:
             nt = self.nontrivial_bonds
             bonds = range(nt.start, nt.stop)
+        if isinstance(bonds, int):
+            bonds = [bonds]
         res = []
         for ib in bonds:
             s = self._S[ib]
@@ -3043,7 +3046,9 @@ def canonical_form_finite(self, renormalize=True, cutoff=0., envs_to_update=None
         Parameters
         ----------
         renormalize: bool
-            Whether a change in the norm should be discarded or used to update :attr:`norm`.
+            Whether a change in the norm should be discarded or used to *update* :attr:`norm`.
+            Note that even `renoramlize=True` *does not reset* the :attr:`norm` to 1.
+            To do that, you would rather have to set ``psi.norm = 1`` explicitly!
         cutoff : float | None
             Cutoff of singular values used in the SVDs.
         envs_to_update : None | list of :class:`MPSEnvironment`
@@ -3183,6 +3188,8 @@ def canonical_form_infinite1(self, renormalize=True, tol_xi=1.e6):
         ----------
         renormalize: bool
             Whether a change in the norm should be discarded or used to update :attr:`norm`.
+            Note that even `renoramlize=True` *does not reset* the :attr:`norm` to 1.
+            To do that, you would rather have to set ``psi.norm = 1`` explicitly!
         tol_xi : float
             Raise an error if the correlation length is larger than that
             (which indicates a degenerate "cat" state, e.g., for spontaneous symmetry breaking).
@@ -3256,6 +3263,8 @@ def canonical_form_infinite2(self, renormalize=True, tol=1.e-15, arnoldi_params=
         ----------
         renormalize: bool
             Whether a change in the norm should be discarded or used to update :attr:`norm`.
+            Note that even `renoramlize=True` *does not reset* the :attr:`norm` to 1.
+            To do that, you would rather have to set ``psi.norm = 1`` explicitly!
         tol : float
             Precision down to which the state should be in canonical form.
         arnoldi_params : dict
diff --git a/tenpy/networks/purification_mps.py b/tenpy/networks/purification_mps.py
index 79f1d3396..e643b6203 100644
--- a/tenpy/networks/purification_mps.py
+++ b/tenpy/networks/purification_mps.py
@@ -116,16 +116,18 @@
     Moreover, we don't split the physical and auxiliar space into separate sites, which makes
     TEBD as costly as :math:`O(d^6 \chi^3)`.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
+import copy
 import numpy as np
 import itertools
 
 from .mps import MPS
 from ..linalg import np_conserved as npc
 from ..tools.math import entropy
+from ..tools.misc import lexsort
 
-__all__ = ['PurificationMPS']
+__all__ = ['PurificationMPS', 'convert_model_purification_canonical_conserve_ancilla_charge']
 
 
 class PurificationMPS(MPS):
@@ -197,32 +199,37 @@ def from_infiniteT(cls, sites, bc='finite', form='B', dtype=np.float64):
         return res
 
     @classmethod
-    def from_infiniteT_canonical(cls, sites, charge_sector, form='B', dtype=np.float64):
+    def from_infiniteT_canonical(cls, sites, charge_sector, dtype=np.float64,
+                                 conserve_ancilla_charge=False):
         """Initial state corresponding to *canonical* infinite-temperature ensemble.
 
         Works only for finite boundary conditions, following the idea outlined in
         :cite:`barthel2016`.
-        However, we just put trivial charges on the ancilla legs,
-        and do *not* double the number of charges as suggested in that paper - there's no need to.
-
-        Note that the 'backwards' disentanglers doesn't work with the canonical ensemble.
 
         Parameters
         ----------
         sites : list of :class:`~tenpy.networks.site.Site`
-            The sites defining the local Hilbert space.
+            The sites defining the local Hilbert space (on the physical legs).
             For usual :class:`tenpy.models.model.Model` given by `model.lat.mps_sites()`.
         charge_sector : tuple of int
             The desired charge sector to be taken for the canonical ensemble.
-        form : (list of) {``'B' | 'A' | 'C' | 'G' | None`` | tuple(float, float)}
-            The canonical form of the stored 'matrices', see table in :mod:`~tenpy.networks.mps`.
-            A single choice holds for all of the entries.
+        dtype : type or string
+            The data type of the array entries.
+        conserve_ancilla_charge : bool
+            Whether to conserve the charges on the ancilla leg.
+            If False, do *not* double the number of conserved charges to get a separate charge
+            for the ancilla degrees of freedom.
+            If True, separately conserve charges on physical and ancilla spaces.
+            In that case, use the function
+            :func:`convert_model_purification_canonical_conserve_ancilla_charge`
+            to get a converted model before using algorithms like the `PurificationTEBD`.
 
         Returns
         -------
         infiniteT_MPS : :class:`PurificationMPS`
             Describes the infinite-temperature (grand canonical) ensemble,
             i.e. expectation values give a trace over all basis states.
+
         """
         sites = list(sites)
         L = len(sites)
@@ -234,54 +241,87 @@ def from_infiniteT_canonical(cls, sites, charge_sector, form='B', dtype=np.float
         charge_sector_right = chinfo.make_valid(charge_sector)
         assert charge_sector_right.ndim == 1
         # get bounds for the maximal and minimal charge values at each bond
-        qflats = [s.leg.to_qflat() for s in sites]
-        min_p_Q = [np.min(qflat, axis=0) for qflat in qflats]
-        max_p_Q = [np.max(qflat, axis=0) for qflat in qflats]
-        min_Q_left = charge_sector_left + np.cumsum(min_p_Q, axis=0)  # on bonds 1, 2, 3... L
-        max_Q_left = charge_sector_left + np.cumsum(max_p_Q, axis=0)
-        min_Q_right = charge_sector_right - np.cumsum(max_p_Q[::-1], axis=0)[::-1]  # 0, 1, ... L-1
-        max_Q_right = charge_sector_right - np.cumsum(min_p_Q[::-1], axis=0)[::-1]
-        min_Q = np.max([min_Q_left[:-1], min_Q_right[1:]], axis=0)  # on non-trivial bonds
-        max_Q = np.min([max_Q_left[:-1], max_Q_right[1:]], axis=0)  # on non-trivial bonds
-        min_Q = np.append(min_Q, [charge_sector_right], axis=0)  # bonds 1, 2, ... L
-        max_Q = np.append(max_Q, [charge_sector_right], axis=0)
-        assert np.all(max_Q >= min_Q)
-        chi = np.prod(max_Q - min_Q + 1, axis=1)  # assumes that charges can be incremented by 1
+        Q_from_right = [None] * L + [set([tuple(charge_sector_right)])]  # all bonds 0, ... L
+        for i in reversed(range(L)):
+            Q_R = np.array(list(Q_from_right[i+1]))
+            # find new charges possible on left of site i, coming from the right
+            Q_L = set()
+            for Q_p in sites[i].leg.charges:
+                Q_L_add = chinfo.make_valid(Q_R - Q_p[np.newaxis, :])
+                Q_L_add = set([tuple(q) for q in Q_L_add])
+                Q_L = Q_L.union(Q_L_add)
+            Q_from_right[i] = Q_L
+        if tuple(charge_sector_left) not in Q_from_right[0]:
+            raise ValueError("can't get desired charge sector {charge_sector!r} "
+                             "for the given charges on physical sites!")
+        Q_from_left = [set([tuple(charge_sector_left)])] + [None] * L
+        Q_L_arrays = []
+        for i in range(L):
+            Q_L = np.array(list(Q_from_left[i]))
+            Q_L = Q_L[lexsort(Q_L.T), :]
+            Q_L_arrays.append(Q_L)
+            Q_R = set()
+            for Q_p in sites[i].leg.charges:
+                Q_R_add = chinfo.make_valid(Q_L + Q_p[:, np.newaxis])
+                Q_R_add = set([tuple(q) for q in Q_R_add])
+                Q_R = Q_R.union(Q_R_add)
+            Q_from_left[i+1] = Q_R.intersection(Q_from_right[i+1])
+        assert Q_from_left[-1] == Q_from_right[-1]  # should match charge_sector on the right
+        Q_L_arrays.append(chinfo.make_valid(charge_sector_right[np.newaxis, :]))
 
-        Bs = []
-        Ss = [np.ones(1, dtype=np.float64)]
         # now we can define the tensors following section VI.C) of [barthel2016]_:
         # B[vL, vR, p, q] = delta_{p,q} delta_{Q(p) + Q(vL), Q(vR)}
         # the normalization will be ensured by a call to `canonical_form_finite()` in the end.
-        right_Q = chinfo.make_valid([charge_sector_left])
-        right_leg = npc.LegCharge.from_qflat(chinfo, right_Q, qconj=-1)
-        for s in range(L):
-            p_leg = sites[s].leg
-            p_Q = p_leg.to_qflat()
-            q_leg = p_leg.conj().copy()
-            q_leg.charges = np.zeros_like(q_leg.charges)
-            left_Q = right_Q
-            left_leg = right_leg.conj()
-            right_Q = np.array(
-                list(
-                    itertools.product(
-                        *[range(min_Q[s][c], max_Q[s][c] + 1) for c in range(chinfo.qnumber)])))
-            right_Q = right_Q[np.lexsort(right_Q.T), :]  # sort charges
-            right_leg = npc.LegCharge.from_qflat(chinfo, right_Q, qconj=-1)
-            B = npc.zeros([left_leg, right_leg, p_leg, q_leg],
+        Bs = []
+        Ss = [np.ones(1, dtype=np.float64)]
+        Q_R = Q_L_arrays[0]
+        if not conserve_ancilla_charge:  # cac := conserve ancilla charges
+            leg_R = npc.LegCharge.from_qflat(chinfo, Q_R, qconj=-1)
+        else:
+            chinfo_cac = npc.ChargeInfo(list(chinfo.mod) * 2,
+                                        chinfo.names + [n + ' ancilla' for n in chinfo.names])
+            Q_R_cac = chinfo_cac.make_valid(np.hstack([Q_R, -Q_R]))
+            leg_R = npc.LegCharge.from_qflat(chinfo_cac, Q_R_cac, qconj=-1)
+            sites_cac = []
+        for i in range(L):
+            leg_p = sites[i].leg
+            Q_p = leg_p.to_qflat()
+            Q_L = Q_L_arrays[i]
+            Q_R = Q_L_arrays[i+1]
+            Q_R_map = dict((tuple(q),i) for i, q in enumerate(Q_R))
+
+            leg_L = leg_R.conj()
+            if not conserve_ancilla_charge:
+                leg_q = leg_p.conj().copy()
+                leg_q.charges = np.zeros_like(leg_q.charges)
+                leg_R = npc.LegCharge.from_qflat(chinfo, Q_R, qconj=-1)
+            else:
+                Q_p_cac = np.hstack([Q_p, np.zeros_like(Q_p)])
+                Q_q_cac = np.hstack([np.zeros_like(Q_p), Q_p])
+                Q_R_cac = chinfo_cac.make_valid(np.hstack([Q_R, -Q_R]))
+                leg_p = npc.LegCharge.from_qflat(chinfo_cac, Q_p_cac, qconj=+1)
+                leg_q = npc.LegCharge.from_qflat(chinfo_cac, Q_q_cac, qconj=-1)
+                leg_R = npc.LegCharge.from_qflat(chinfo_cac, Q_R_cac, qconj=-1)
+                s_cac = copy.copy(sites[i])
+                s_cac.change_charge(leg_p)  # note: if Q_p is sorted, so are Q_{p,q}_cac
+                sites_cac.append(s_cac)
+            B = npc.zeros([leg_L, leg_R, leg_p, leg_q],
                           dtype=dtype,
                           labels=['vL', 'vR', 'p', 'q'])
-            for p in range(p_leg.ind_len):
-                for vL in range(left_Q.shape[0]):
-                    Q_vR = left_Q[vL] + p_Q[p]
-                    vR = np.nonzero(np.all(Q_vR == right_Q, axis=1))[0]
-                    if len(vR) == 0:
-                        continue
-                    vR = vR.item()
-                    B[vL, vR, p, p] = 1.  # add an entry in the tensor
+            for j in range(leg_p.ind_len):
+                Q_p_j = Q_p[j]
+                for vL, Q_L_vL in enumerate(Q_L):
+                    Q_R_vR = tuple(chinfo.make_valid(Q_L_vL + Q_p_j))
+                    vR = Q_R_map.get(Q_R_vR, None)
+                    if vR is not None:
+                        B[vL, vR, j, j] = 1.  # add an entry in the tensor
+                    # else: dropped Q_R_vR since it can't reach charge_sector on the right any more
             Bs.append(B)
             Ss.append(np.ones(B.shape[1], np.float64))
-        res = cls(sites, Bs, Ss, 'finite', form)
+
+        if conserve_ancilla_charge:
+            sites = sites_cac
+        res = cls(sites, Bs, Ss, 'finite', form='B')
         res.canonical_form_finite()  # calculate S values and normalize
         return res
 
@@ -466,3 +506,85 @@ def _get_p_labels(self, ks, star=False):
             return [lbl + str(k) + '*' for k in range(ks) for lbl in self._p_label]
         else:
             return [lbl + str(k) for k in range(ks) for lbl in self._p_label]
+
+
+def convert_model_purification_canonical_conserve_ancilla_charge(model):
+    """Extend charges of model for :meth:`PurificationMPS.from_infiniteT_canonical`.
+
+    Parameters
+    ----------
+    model : :class:`tenpy.models.model.Model`
+        Model to be converted.
+
+    Returns
+    -------
+    model_with_extra_charges : :class:`tenpy.models.model.Model`
+        Shallow copy of the `model` with charges of sites, `H_MPO` and `H_bond` adjusted
+        to fit the doubled (with 0 extended) charges of the canonical ensemble of the
+        :class:`PurificationMPS`. The number of
+    """
+    # cac := conserve_ancilla_charge
+    model = model.copy()
+    chinfo = model.lat.unit_cell[0].leg.chinfo
+    chinfo_cac = npc.ChargeInfo(list(chinfo.mod) * 2,
+                                chinfo.names + [n + ' ancilla' for n in chinfo.names])
+
+    converted_sites_cache = {}
+    def _convert_site(site):
+        s = converted_sites_cache.get(site, None)
+        if s is not None:
+            return s
+        s_new = copy.copy(site)
+        leg_p = s_new.leg
+        Q_p = leg_p.charges
+        Q_p_cac = np.hstack([Q_p, np.zeros_like(Q_p)])  # still sorted
+        leg_p_cac = npc.LegCharge(chinfo_cac, leg_p.slices, Q_p_cac, leg_p.qconj)
+        s_new.change_charge(leg_p_cac)
+        converted_sites_cache[site] = s_new
+        return s_new
+
+    model.lat = model.lat.copy()
+    model.lat.unit_cell = [_convert_site(s) for s in model.lat.unit_cell]
+
+    if hasattr(model, 'H_MPO'):
+        model.H_MPO = H_MPO = model.H_MPO.copy()
+        H_MPO.sites = [_convert_site(s) for s in H_MPO.sites]
+        H_MPO.chinfo = chinfo_cac
+        new_W = []
+        for W in H_MPO._W:
+            W = W.copy()
+            W.itranspose(['wL', 'wR', 'p', 'p*'])
+            W.legs = W.legs[:]
+            for i in range(3):
+                leg = W.legs[i]
+                if i < 2:
+                    Q = np.hstack([leg.charges, -leg.charges])  # wL, wR
+                else:
+                    Q = np.hstack([leg.charges, np.zeros_like(leg.charges)])  # p
+                W.legs[i] = npc.LegCharge(chinfo_cac,
+                                          leg.slices,
+                                          chinfo_cac.make_valid(Q),
+                                          leg.qconj)
+            W.qtotal = np.hstack([W.qtotal, np.zeros_like(W.qtotal)])
+            W.legs[3] = W.legs[2].conj()
+            new_W.append(W)
+        H_MPO._W = new_W
+
+    if hasattr(model, 'H_bond'):
+        sites = model.lat.mps_sites()  # already updated!
+        model.H_bond = H_bond = model.H_bond[:]
+        L = len(sites)
+        assert len(sites) == len(H_bond)
+        for i, H in enumerate(H_bond):
+            if H is None:
+                continue
+            leg_p0 = sites[(i-1) % L].leg
+            leg_p1 = sites[i].leg
+            H = H.transpose(['p0', 'p1', 'p0*', 'p1*'])  # copy!
+            H.chinfo = chinfo_cac
+            H.legs = [leg_p0, leg_p1, leg_p0.conj(), leg_p1.conj()]
+            H.qtotal = np.hstack([H.qtotal, np.zeros_like(H.qtotal)])
+            H.test_sanity()
+            H_bond[i] = H
+    model.test_sanity()
+    return model
diff --git a/tenpy/networks/site.py b/tenpy/networks/site.py
index 152034aa2..9073d58ea 100644
--- a/tenpy/networks/site.py
+++ b/tenpy/networks/site.py
@@ -3,7 +3,7 @@
 The :class:`Site` is the prototype, read it's docstring.
 
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import itertools
diff --git a/tenpy/networks/terms.py b/tenpy/networks/terms.py
index 86386ec5b..0fd0b1614 100644
--- a/tenpy/networks/terms.py
+++ b/tenpy/networks/terms.py
@@ -5,7 +5,7 @@
 sites it acts on. Moreover, we associate a `strength` to each term, which corresponds to the
 prefactor when specifying e.g. a Hamiltonian.
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import warnings
diff --git a/tenpy/simulations/__init__.py b/tenpy/simulations/__init__.py
index 9a55fb73b..d905b9949 100644
--- a/tenpy/simulations/__init__.py
+++ b/tenpy/simulations/__init__.py
@@ -12,7 +12,7 @@
     ground_state_search
     time_evolution
 """
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 from . import measurement, simulation, ground_state_search, time_evolution
 
diff --git a/tenpy/simulations/ground_state_search.py b/tenpy/simulations/ground_state_search.py
index 8be0bbbc0..0d3504c8e 100644
--- a/tenpy/simulations/ground_state_search.py
+++ b/tenpy/simulations/ground_state_search.py
@@ -1,5 +1,5 @@
 """Simulations for ground state searches."""
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 from pathlib import Path
diff --git a/tenpy/simulations/measurement.py b/tenpy/simulations/measurement.py
index 8b5c97461..e16e58140 100644
--- a/tenpy/simulations/measurement.py
+++ b/tenpy/simulations/measurement.py
@@ -10,7 +10,7 @@
 
 Full description and details in :doc:`/intro/measurements`.
 """
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 from ..networks.mpo import MPOEnvironment
 from ..tools.misc import get_recursive
diff --git a/tenpy/simulations/simulation.py b/tenpy/simulations/simulation.py
index 346d6553f..6cd8e6249 100644
--- a/tenpy/simulations/simulation.py
+++ b/tenpy/simulations/simulation.py
@@ -7,7 +7,7 @@
 See :doc:`/intro/simulations` for an overview and
 :doc:`/examples` for a list of example parameter yaml files.
 """
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import os
 import sys
diff --git a/tenpy/simulations/time_evolution.py b/tenpy/simulations/time_evolution.py
index 0079f2118..0f7770a9d 100644
--- a/tenpy/simulations/time_evolution.py
+++ b/tenpy/simulations/time_evolution.py
@@ -1,5 +1,5 @@
 """Simulations for (real) time evolution."""
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tenpy/tools/__init__.py b/tenpy/tools/__init__.py
index 3f68d7ad4..737429d3d 100644
--- a/tenpy/tools/__init__.py
+++ b/tenpy/tools/__init__.py
@@ -21,7 +21,7 @@
     cache
     thread
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from . import (events, fit, hdf5_io, math, misc, params, process, string, optimization, cache,
                thread)
diff --git a/tenpy/tools/cache.py b/tenpy/tools/cache.py
index 4a1f8f00d..a035bfe9f 100644
--- a/tenpy/tools/cache.py
+++ b/tenpy/tools/cache.py
@@ -7,7 +7,7 @@
 this is easiest done through a ``with`` statement, see the example in :class:`DictCache`.
 """
 
-# Copyright 2021 TeNPy Developers, GNU GPLv3
+# Copyright 2021-2023 TeNPy Developers, GNU GPLv3
 
 import pickle
 import numpy as np
diff --git a/tenpy/tools/events.py b/tenpy/tools/events.py
index d2b54be46..76a523d15 100644
--- a/tenpy/tools/events.py
+++ b/tenpy/tools/events.py
@@ -4,7 +4,7 @@
 which can get called once a certain "event" happens.
 Examples are given in the class doc-string.
 """
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 from collections import namedtuple
 import warnings
diff --git a/tenpy/tools/fit.py b/tenpy/tools/fit.py
index 9266ba687..fd6e21c36 100644
--- a/tenpy/tools/fit.py
+++ b/tenpy/tools/fit.py
@@ -1,5 +1,5 @@
 """tools to fit to an algebraic decay."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import scipy.optimize as optimize
diff --git a/tenpy/tools/hdf5_io.py b/tenpy/tools/hdf5_io.py
index f5211d11a..53e841c1b 100644
--- a/tenpy/tools/hdf5_io.py
+++ b/tenpy/tools/hdf5_io.py
@@ -65,7 +65,7 @@
 
 .. autodata:: TYPES_FOR_HDF5_DATASETS
 """
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import pickle
 import gzip
diff --git a/tenpy/tools/math.py b/tenpy/tools/math.py
index 9c6478036..84d78da1c 100644
--- a/tenpy/tools/math.py
+++ b/tenpy/tools/math.py
@@ -2,7 +2,7 @@
 
 .. autodata:: LeviCivita3
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import warnings
diff --git a/tenpy/tools/misc.py b/tenpy/tools/misc.py
index e72bb181d..c6eef9753 100644
--- a/tenpy/tools/misc.py
+++ b/tenpy/tools/misc.py
@@ -1,5 +1,5 @@
 """Miscellaneous tools, somewhat random mix yet often helpful."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import logging
 import numpy as np
diff --git a/tenpy/tools/optimization.py b/tenpy/tools/optimization.py
index a02a859fd..91f9caf17 100644
--- a/tenpy/tools/optimization.py
+++ b/tenpy/tools/optimization.py
@@ -76,7 +76,7 @@
 .. autodata:: have_cython_functions
 .. autodata:: compiled_with_MKL
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from enum import IntEnum
 import warnings
diff --git a/tenpy/tools/params.py b/tenpy/tools/params.py
index 879238d81..329f5bbb1 100644
--- a/tenpy/tools/params.py
+++ b/tenpy/tools/params.py
@@ -2,7 +2,7 @@
 
 See the doc-string of :class:`Config` for details.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import warnings
 import numpy as np
diff --git a/tenpy/tools/process.py b/tenpy/tools/process.py
index dfc09c3f3..3c0f8e187 100644
--- a/tenpy/tools/process.py
+++ b/tenpy/tools/process.py
@@ -13,7 +13,7 @@
 while still being failsave if the shared OpenMP library is not found.  In the latter case,
 you might also try the equivalent :func:`mkl_get_nthreads` and :func:`mkl_set_nthreads`.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import warnings
 import ctypes
diff --git a/tenpy/tools/string.py b/tenpy/tools/string.py
index 180c58937..141d8b3e7 100644
--- a/tenpy/tools/string.py
+++ b/tenpy/tools/string.py
@@ -1,6 +1,6 @@
 """Tools for handling strings."""
 
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 __all__ = ['is_non_string_iterable', 'vert_join', 'to_mathematica_lists']
 
diff --git a/tenpy/tools/thread.py b/tenpy/tools/thread.py
index 44c142fac..a33a608b9 100644
--- a/tenpy/tools/thread.py
+++ b/tenpy/tools/thread.py
@@ -1,6 +1,6 @@
 """Tools for thread-based parallelization."""
 
-# Copyright 2021 TeNPy Developers, GNU GPLv3
+# Copyright 2021-2023 TeNPy Developers, GNU GPLv3
 
 import threading
 import queue
diff --git a/tenpy/version.py b/tenpy/version.py
index d8c40cb0e..188a875f6 100644
--- a/tenpy/version.py
+++ b/tenpy/version.py
@@ -10,7 +10,7 @@
 .. autodata :: full_version
 .. autodata :: version_summary
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import sys
 import subprocess
diff --git a/toycodes/a_mps.py b/toycodes/a_mps.py
index 11af397b9..d2c8428dd 100644
--- a/toycodes/a_mps.py
+++ b/toycodes/a_mps.py
@@ -1,5 +1,5 @@
 """Toy code implementing a matrix product state."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 from scipy.linalg import svd
diff --git a/toycodes/b_model.py b/toycodes/b_model.py
index fde9335fe..3ec2d5547 100644
--- a/toycodes/b_model.py
+++ b/toycodes/b_model.py
@@ -1,5 +1,5 @@
 """Toy code implementing the transverse-field ising model."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/toycodes/c_tebd.py b/toycodes/c_tebd.py
index bd29bc042..fa06750f3 100644
--- a/toycodes/c_tebd.py
+++ b/toycodes/c_tebd.py
@@ -1,5 +1,5 @@
 """Toy code implementing the time evolving block decimation (TEBD)."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 from scipy.linalg import expm
diff --git a/toycodes/d_dmrg.py b/toycodes/d_dmrg.py
index 059074949..d035b7c06 100644
--- a/toycodes/d_dmrg.py
+++ b/toycodes/d_dmrg.py
@@ -1,5 +1,5 @@
 """Toy code implementing the density-matrix renormalization group (DMRG)."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 from a_mps import split_truncate_theta
diff --git a/toycodes/e_tdvp.py b/toycodes/e_tdvp.py
index fbf80a4d4..534aae132 100644
--- a/toycodes/e_tdvp.py
+++ b/toycodes/e_tdvp.py
@@ -1,5 +1,5 @@
 """Toy code implementing the TDVP for *finite* MPS."""
-# Copyright 2021-2022 TeNPy Developers, GNU GPLv3
+# Copyright 2021-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import scipy.sparse.linalg
diff --git a/toycodes/lanczos.py b/toycodes/lanczos.py
index 1f0f00ab3..f3dcf9396 100644
--- a/toycodes/lanczos.py
+++ b/toycodes/lanczos.py
@@ -3,7 +3,7 @@
 These functions can be used in d_dmrg.py and e_tdvp to replace
 the eigsh and expm_mulitply calls, as commented out in the code.
 """
-# Copyright 2021-2022 TeNPy Developers, GNU GPLv3
+# Copyright 2021-2023 TeNPy Developers, GNU GPLv3
 
 from scipy.sparse.linalg import expm
 import numpy as np
diff --git a/toycodes/tfi_exact.py b/toycodes/tfi_exact.py
index bd995e82d..42d748429 100644
--- a/toycodes/tfi_exact.py
+++ b/toycodes/tfi_exact.py
@@ -4,7 +4,7 @@
 .. math ::
     H = - J \\sum_{i} \\sigma^x_i \\sigma^x_{i+1} - g \\sum_{i} \\sigma^z_i
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import scipy.sparse as sparse


Test patch:
diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index 2b87e6416..2c3ddf6e6 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -1,6 +1,10 @@
 name: pytest
 
-on: [push]
+on:
+  push: 
+    branches:
+      - main
+  pull_request:
 
 jobs:
   build:
diff --git a/doc/changelog/latest.rst b/doc/changelog/latest.rst
new file mode 100644
index 000000000..234e2c099
--- /dev/null
+++ b/doc/changelog/latest.rst
@@ -0,0 +1,27 @@
+[latest]
+========
+
+Release Notes
+-------------
+TODO: Summarize the most important changes
+
+Changelog
+---------
+
+Backwards incompatible changes
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+- nothing yet
+
+Added
+^^^^^
+- nothing yet
+
+Changed
+^^^^^^^
+- nothing yet
+
+Fixed
+^^^^^
+- :meth:`~tenpy.networks.purification_mps.PurificationMPS.from_infiniteT_canonical` should now work with arbitrary
+  charges of the original model, and has the option to double the number of charges to separately conserve the charges
+  on each the physical and ancialla legs.
diff --git a/tests/benchmark/benchmark.py b/tests/benchmark/benchmark.py
index ff50c1626..c6c28a57c 100755
--- a/tests/benchmark/benchmark.py
+++ b/tests/benchmark/benchmark.py
@@ -11,7 +11,7 @@
 Afterwards, you can plot the produced files::
     python benchmark.py -p tensordot_*_benchmark_*.txt
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from __future__ import print_function  # (backwards compatibility to python 2)
 
diff --git a/tests/benchmark/combine_npc.py b/tests/benchmark/combine_npc.py
index 264a98240..bedcdca28 100644
--- a/tests/benchmark/combine_npc.py
+++ b/tests/benchmark/combine_npc.py
@@ -1,5 +1,5 @@
 """To be used in the `-m` argument of benchmark.py."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import tenpy.linalg.np_conserved as npc
diff --git a/tests/benchmark/combine_numpy.py b/tests/benchmark/combine_numpy.py
index f8bffbe06..e83dde703 100644
--- a/tests/benchmark/combine_numpy.py
+++ b/tests/benchmark/combine_numpy.py
@@ -1,5 +1,5 @@
 """To be used in the `-m` argument of benchmark.py."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import tensordot_npc
diff --git a/tests/benchmark/dmrg_infinite.py b/tests/benchmark/dmrg_infinite.py
index 3fcfb3ae9..92d8aff4e 100644
--- a/tests/benchmark/dmrg_infinite.py
+++ b/tests/benchmark/dmrg_infinite.py
@@ -1,5 +1,5 @@
 """To be used in the `-m` argument of benchmark.py."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tests/benchmark/profiling.py b/tests/benchmark/profiling.py
index 78e89178a..72e061b42 100755
--- a/tests/benchmark/profiling.py
+++ b/tests/benchmark/profiling.py
@@ -8,7 +8,7 @@
 
     python profiling.py -p tensordot_*_profile_*.prof --sort time --limit 10
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import cProfile
diff --git a/tests/benchmark/singlerun.py b/tests/benchmark/singlerun.py
index 8c1a16bd5..649c3ad20 100755
--- a/tests/benchmark/singlerun.py
+++ b/tests/benchmark/singlerun.py
@@ -4,7 +4,7 @@
 -s 20 -S 50 -q 1 1
 """
 
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import cProfile
diff --git a/tests/benchmark/split_npc.py b/tests/benchmark/split_npc.py
index ad9b0b4de..114f0e246 100644
--- a/tests/benchmark/split_npc.py
+++ b/tests/benchmark/split_npc.py
@@ -1,5 +1,5 @@
 """To be used in the `-m` argument of benchmark.py."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import tenpy.linalg.np_conserved as npc
diff --git a/tests/benchmark/split_numpy.py b/tests/benchmark/split_numpy.py
index 5e56fcb42..d33f69175 100644
--- a/tests/benchmark/split_numpy.py
+++ b/tests/benchmark/split_numpy.py
@@ -1,5 +1,5 @@
 """To be used in the `-m` argument of benchmark.py."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import tensordot_npc
diff --git a/tests/benchmark/tebd_infinite.py b/tests/benchmark/tebd_infinite.py
index f183c37ce..8ba456557 100644
--- a/tests/benchmark/tebd_infinite.py
+++ b/tests/benchmark/tebd_infinite.py
@@ -1,5 +1,5 @@
 """To be used in the `-m` argument of benchmark.py."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 
diff --git a/tests/benchmark/tensordot_npc.py b/tests/benchmark/tensordot_npc.py
index b8b4db20f..cd9cc1d2b 100644
--- a/tests/benchmark/tensordot_npc.py
+++ b/tests/benchmark/tensordot_npc.py
@@ -1,5 +1,5 @@
 """To be used in the `-m` argument of benchmark.py."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import tenpy.linalg.np_conserved as npc
diff --git a/tests/benchmark/tensordot_numpy.py b/tests/benchmark/tensordot_numpy.py
index a8471f0f0..332d24a58 100644
--- a/tests/benchmark/tensordot_numpy.py
+++ b/tests/benchmark/tensordot_numpy.py
@@ -1,5 +1,5 @@
 """To be used in the `-m` argument of benchmark.py."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import tensordot_npc
diff --git a/tests/export_import_test/io_test.py b/tests/export_import_test/io_test.py
index 6b1e59bfa..ec9705649 100644
--- a/tests/export_import_test/io_test.py
+++ b/tests/export_import_test/io_test.py
@@ -7,7 +7,7 @@
 test_pickle.py``. This will generate the files with pre-defined data (see :func:`gen_example_data`)
 and the tenpy version in the filename.
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import os
diff --git a/tests/random_test.py b/tests/random_test.py
index 346190b21..ccef10817 100644
--- a/tests/random_test.py
+++ b/tests/random_test.py
@@ -1,5 +1,5 @@
 """Provide helper functions for test of random Arrays."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import itertools as it
diff --git a/tests/tdvp_numpy.py b/tests/tdvp_numpy.py
index 4414229ee..9de5442df 100644
--- a/tests/tdvp_numpy.py
+++ b/tests/tdvp_numpy.py
@@ -2,7 +2,7 @@
 
 .. todo ::     Make this a nice toycode! For the test, use some pre-defined values.
 """
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 from scipy.linalg import expm
 from scipy.sparse.linalg import expm_multiply
diff --git a/tests/test_charges.py b/tests/test_charges.py
index 44cf00250..b55950eab 100644
--- a/tests/test_charges.py
+++ b/tests/test_charges.py
@@ -1,5 +1,5 @@
 """A collection of tests for tenpy.linalg.charges."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy.linalg.charges as charges
 import numpy as np
diff --git a/tests/test_consistencies.py b/tests/test_consistencies.py
index 7303f2db2..3ddb571f0 100644
--- a/tests/test_consistencies.py
+++ b/tests/test_consistencies.py
@@ -1,10 +1,11 @@
 """Check for consistencies."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy
 import types
 import os
 import re
+from datetime import date
 
 
 def test_all(check_module=tenpy):
@@ -70,7 +71,8 @@ def test_copyright():
     #  you can use the following:
     # tenpy_files = get_python_files(os.path.dirname(os.path.dirname(tenpy.__file__)))
     #  (but this doesn't work for the pip-installed tenpy, so you can only do it temporary!)
-    regex = re.compile(r'#\s[Cc]opyright 20[0-9\-]+\s+(TeNPy|tenpy) [dD]evelopers, GNU GPLv3')
+    current_year = date.today().year
+    regex = re.compile(r'#\s[Cc]opyright (20[0-9]{2}-)?' + str(current_year) + r' (TeNPy|tenpy) [dD]evelopers, GNU GPLv3')
     for fn in tenpy_files:
         with open(fn, 'r') as f:
             for line in f:
@@ -79,4 +81,4 @@ def test_copyright():
                     if match is not None:
                         break
             else:  # no break
-                raise AssertionError("No/wrong copyright notice in {0!s}".format(fn))
+                raise AssertionError(f'No/wrong copyright notice in {fn}.')
diff --git a/tests/test_dmrg.py b/tests/test_dmrg.py
index 61fbb5c19..1b42f182c 100644
--- a/tests/test_dmrg.py
+++ b/tests/test_dmrg.py
@@ -1,5 +1,5 @@
 """A collection of tests to check the functionality of `tenpy.dmrg`"""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import itertools as it
 import tenpy.linalg.np_conserved as npc
diff --git a/tests/test_exact_diag.py b/tests/test_exact_diag.py
index ebb52539f..7c0d94768 100644
--- a/tests/test_exact_diag.py
+++ b/tests/test_exact_diag.py
@@ -1,5 +1,5 @@
 """A collection of tests to check the functionality of algorithms.exact_diagonalization."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy.linalg.np_conserved as npc
 import numpy as np
diff --git a/tests/test_lanczos.py b/tests/test_lanczos.py
index f52ff0b83..a64f6bbd3 100644
--- a/tests/test_lanczos.py
+++ b/tests/test_lanczos.py
@@ -1,5 +1,5 @@
 """A collection of tests for tenpy.linalg.lanczos."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy.linalg.np_conserved as npc
 import numpy as np
diff --git a/tests/test_lattice.py b/tests/test_lattice.py
index 42e217349..b5118e88a 100644
--- a/tests/test_lattice.py
+++ b/tests/test_lattice.py
@@ -1,5 +1,5 @@
 """A collection of tests for tenpy.models.lattice."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from tenpy.models import lattice
 import tenpy.linalg.np_conserved as npc
diff --git a/tests/test_model.py b/tests/test_model.py
index be8e5ea5f..b8733f7dd 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -1,5 +1,5 @@
 """A collection of tests for (classes in) :mod:`tenpy.models.model`."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import warnings
 import itertools
diff --git a/tests/test_model_aklt.py b/tests/test_model_aklt.py
index 7a92cd0b0..40b222c69 100644
--- a/tests/test_model_aklt.py
+++ b/tests/test_model_aklt.py
@@ -1,4 +1,4 @@
-# Copyright 2021 TeNPy Developers, GNU GPLv3
+# Copyright 2021-2023 TeNPy Developers, GNU GPLv3
 from tenpy.models import aklt
 from test_model import check_general_model
 from tenpy.algorithms.dmrg import TwoSiteDMRGEngine
diff --git a/tests/test_model_fermions_spinless.py b/tests/test_model_fermions_spinless.py
index 685587832..58ae358a7 100644
--- a/tests/test_model_fermions_spinless.py
+++ b/tests/test_model_fermions_spinless.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 import numpy as np
 from tenpy.models.fermions_spinless import FermionChain, FermionModel
 from test_model import check_general_model
diff --git a/tests/test_model_haldane.py b/tests/test_model_haldane.py
index a3875bf3b..b9b736af4 100644
--- a/tests/test_model_haldane.py
+++ b/tests/test_model_haldane.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 from tenpy.models.haldane import BosonicHaldaneModel, FermionicHaldaneModel
 from test_model import check_general_model
 import pytest
diff --git a/tests/test_model_hofstadter.py b/tests/test_model_hofstadter.py
index da02f3b31..b7b2dd5ec 100644
--- a/tests/test_model_hofstadter.py
+++ b/tests/test_model_hofstadter.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 from tenpy.models.hofstadter import HofstadterBosons, HofstadterFermions
 from test_model import check_general_model
 import pytest
diff --git a/tests/test_model_hubbard.py b/tests/test_model_hubbard.py
index 55dd17d1e..afbbd568f 100644
--- a/tests/test_model_hubbard.py
+++ b/tests/test_model_hubbard.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 from tenpy.models import hubbard
 from test_model import check_general_model
 
diff --git a/tests/test_model_mixed_xk.py b/tests/test_model_mixed_xk.py
index 0f2132c8f..13bec794a 100644
--- a/tests/test_model_mixed_xk.py
+++ b/tests/test_model_mixed_xk.py
@@ -1,3 +1,4 @@
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 from tenpy.models.mixed_xk import (MixedXKLattice, MixedXKModel, SpinlessMixedXKSquare,
                                    HubbardMixedXKSquare)
 from test_model import check_general_model
diff --git a/tests/test_model_spins.py b/tests/test_model_spins.py
index 6e1b382ee..b28b6d402 100644
--- a/tests/test_model_spins.py
+++ b/tests/test_model_spins.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 from tenpy.models import spins
 from test_model import check_general_model
 
diff --git a/tests/test_model_spins_nnn.py b/tests/test_model_spins_nnn.py
index 2919e4639..04da22ed1 100644
--- a/tests/test_model_spins_nnn.py
+++ b/tests/test_model_spins_nnn.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 from tenpy.models import spins_nnn
 from test_model import check_general_model
 from tenpy.models.model import NearestNeighborModel
diff --git a/tests/test_model_tf_ising.py b/tests/test_model_tf_ising.py
index 0c3e6def4..52a37f147 100644
--- a/tests/test_model_tf_ising.py
+++ b/tests/test_model_tf_ising.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from tenpy.models.tf_ising import TFIModel, TFIChain
 from test_model import check_general_model
diff --git a/tests/test_model_toric_code.py b/tests/test_model_toric_code.py
index 4f5845d6e..4b120bbfc 100644
--- a/tests/test_model_toric_code.py
+++ b/tests/test_model_toric_code.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 from tenpy.models.toric_code import ToricCode
 from test_model import check_general_model
diff --git a/tests/test_model_xxz_chain.py b/tests/test_model_xxz_chain.py
index 093ed902b..9367252ea 100644
--- a/tests/test_model_xxz_chain.py
+++ b/tests/test_model_xxz_chain.py
@@ -1,5 +1,5 @@
 """test of :class:`tenpy.models.XXZChain`."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import numpy.testing as npt
diff --git a/tests/test_mpo.py b/tests/test_mpo.py
index 79871b3f4..1dbb55c52 100644
--- a/tests/test_mpo.py
+++ b/tests/test_mpo.py
@@ -1,5 +1,5 @@
 """A collection of tests for (classes in) :module:`tenpy.networks.mpo`."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import numpy.testing as npt
diff --git a/tests/test_mps.py b/tests/test_mps.py
index 79deec0f5..427fb8536 100644
--- a/tests/test_mps.py
+++ b/tests/test_mps.py
@@ -1,5 +1,5 @@
 """A collection of tests for :module:`tenpy.networks.mps`."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import numpy.testing as npt
diff --git a/tests/test_network_contractor.py b/tests/test_network_contractor.py
index 8d169be5d..cd7625abd 100644
--- a/tests/test_network_contractor.py
+++ b/tests/test_network_contractor.py
@@ -1,11 +1,12 @@
 """a collection of tests to check the functionality of network_contractor.py."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
-from tenpy.algorithms.network_contractor import contract, outer_product
+from tenpy.algorithms.network_contractor import contract, ncon
 import numpy as np
 from tenpy.linalg import np_conserved as npc
 import pytest
 import warnings
+from random_test import random_Array
 
 # Contruct toy tensors
 # ====================
@@ -103,23 +104,13 @@ def test_contract_with_sequence():
     S = Sy
     S.iset_leg_labels(['U', 'L'])
 
-    # the supplied sequence should generate three suboptimal warnings
-    with warnings.catch_warnings(record=True) as cw:
-        # Cause all warnings to always be triggered.
-        warnings.simplefilter("always")
-
-        res = contract(tensor_list=[v, h2, S, h, w],
-                       tensor_names=['v', 'h2', 'S', 'h', 'w'],
-                       leg_contractions=[['v', 'L1', 'h2', 'p1*'], ['v', 'L2', 'h2', 'p2*'],
-                                         ['h2', 'p1', 'h', 'p1*'], ['h2', 'p2', 'S', 'U'],
-                                         ['S', 'L', 'h', 'p2*'], ['h', 'p1', 'w', 'U1'],
-                                         ['h', 'p2', 'w', 'U2']],
-                       sequence=[1, 3, 5, 6, 4, 2, 0])
-
-        assert len(cw) == 3
-        assert "Suboptimal contraction sequence" in str(cw[0].message)
-        assert "Suboptimal contraction sequence" in str(cw[1].message)
-        assert "Suboptimal contraction sequence" in str(cw[2].message)
+    res = contract(tensor_list=[v, h2, S, h, w],
+                tensor_names=['v', 'h2', 'S', 'h', 'w'],
+                leg_contractions=[['v', 'L1', 'h2', 'p1*'], ['v', 'L2', 'h2', 'p2*'],
+                                    ['h2', 'p1', 'h', 'p1*'], ['h2', 'p2', 'S', 'U'],
+                                    ['S', 'L', 'h', 'p2*'], ['h', 'p1', 'w', 'U1'],
+                                    ['h', 'p2', 'w', 'U2']],
+                sequence=[1, 3, 5, 6, 4, 2, 0])
 
     expected_result = -0.1735 - 0.5015j  # from MatLab
 
@@ -163,3 +154,24 @@ def test_outer_product():
     expected_result = expected_result.transpose([1, 3, 0, 2])
 
     assert np.linalg.norm(res.to_ndarray() - expected_result) < 1.e-10
+    
+    
+def test_fixes_issue131():
+    no_charge = npc.Array.from_ndarray_trivial(np.random.rand(8, 2, 4))
+    with_charge = random_Array((20, 15, 10), npc.ChargeInfo([1, 3, 1, 2]), sort=False)
+    
+    for A in [no_charge, with_charge]:
+
+        links1 = [[1, -1, -3], [1, -2, -4]]
+        # sum_x A[x, a, c] A[x, b, d] = result[a, b, c, d]
+        expect1 = npc.tensordot(A, A.conj(), (0, 0)).transpose([0, 2, 1, 3])
+        res1 = ncon([A, A.conj()], links1)
+        assert npc.norm(expect1 - res1) < 1e-10
+        
+        links2 = [[1, -1, -2], [1, -3, -4]]
+        # sum_x A[x, a, b] A[x, c, d] = result[a, b, c, d]
+        expect2 = npc.tensordot(A, A.conj(), (0, 0))
+        res2 = ncon([A, A.conj()], links2)
+        assert npc.norm(expect2 - res2) < 1e-10
+
+        assert res1.shape != res2.shape
diff --git a/tests/test_np_conserved.py b/tests/test_np_conserved.py
index 01cee50a9..e85faef7e 100644
--- a/tests/test_np_conserved.py
+++ b/tests/test_np_conserved.py
@@ -1,5 +1,5 @@
 """A collection of tests for tenpy.linalg.np_conserved."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import tenpy.linalg.np_conserved as npc
 import numpy as np
diff --git a/tests/test_params.py b/tests/test_params.py
index 24f3b07f9..2c2d3a070 100644
--- a/tests/test_params.py
+++ b/tests/test_params.py
@@ -1,5 +1,5 @@
 """A test for tenpy.tools.params."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import warnings
 from tenpy.tools.params import Config, asConfig
diff --git a/tests/test_purification.py b/tests/test_purification.py
index 3d02333bc..37d00ce9b 100644
--- a/tests/test_purification.py
+++ b/tests/test_purification.py
@@ -1,5 +1,5 @@
 """A collection of tests for :module:`tenpy.networks.purification_mps`."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import warnings
 import numpy as np
@@ -46,10 +46,15 @@ def test_purification_mps():
             npt.assert_allclose(C, 0.5 * 0.5 * np.eye(L), atol=1.e-13)
 
 
-def test_canoncial_purification(L=6, charge_sector=0, eps=1.e-14):
+
+@pytest.mark.parametrize('conserve_ancilla', [False, True])
+def test_canoncial_purification(conserve_ancilla, L=6, charge_sector=0, eps=1.e-14):
     site = spin_half
-    psi = purification_mps.PurificationMPS.from_infiniteT_canonical([site] * L, [charge_sector])
+    psi = purification_mps.PurificationMPS.from_infiniteT_canonical(
+        [site] * L, [charge_sector], conserve_ancilla_charge=conserve_ancilla)
     psi.test_sanity()
+    Szs = psi.expectation_value('Sz')
+    assert abs(sum(Szs) - charge_sector) < 1.e-13
     total_psi = psi.get_theta(0, L).take_slice(0, 'vL').take_slice(0, 'vR')
     total_psi.itranspose(['p' + str(i) for i in range(L)] + ['q' + str(i) for i in range(L)])
     # note: don't `combine_legs`: it will permute the p legs differently than q due to charges
@@ -80,6 +85,8 @@ def test_canoncial_purification(L=6, charge_sector=0, eps=1.e-14):
         'dt': 0.1,
         'N_steps': 2
     }
+    if conserve_ancilla:
+        M = purification_mps.convert_model_purification_canonical_conserve_ancilla_charge(M)
     eng = PurificationTEBD(psi, M, TEBD_params)
     eng.run_imaginary(0.2)
     eng.run()
diff --git a/tests/test_random_matrix.py b/tests/test_random_matrix.py
index d75ddca00..803ff96b4 100644
--- a/tests/test_random_matrix.py
+++ b/tests/test_random_matrix.py
@@ -1,5 +1,5 @@
 """A collection of tests for :module:`tenpy.linalg.random_matrix`."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import numpy.testing as npt
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 0812bff7b..7e08da97b 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -1,5 +1,5 @@
 """A collection of tests to check the functionality of modules in `tenpy.simulations`"""
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import copy
 import numpy as np
diff --git a/tests/test_site.py b/tests/test_site.py
index a716c2428..bdb31cc35 100644
--- a/tests/test_site.py
+++ b/tests/test_site.py
@@ -1,5 +1,5 @@
 """A collection of tests for :mod:`tenpy.models.site`."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import numpy.testing as npt
diff --git a/tests/test_sparse.py b/tests/test_sparse.py
index e9c3662d2..aaf1ac626 100644
--- a/tests/test_sparse.py
+++ b/tests/test_sparse.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 import tenpy.linalg.np_conserved as npc
 import numpy as np
 
diff --git a/tests/test_svd_robust.py b/tests/test_svd_robust.py
index e7156a77c..e2d7d1b97 100644
--- a/tests/test_svd_robust.py
+++ b/tests/test_svd_robust.py
@@ -1,5 +1,5 @@
 """A collection of tests for tenpy.linalg.svd_robust."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import numpy.testing as npt
diff --git a/tests/test_tdvp.py b/tests/test_tdvp.py
index 30292b764..f47182c6f 100644
--- a/tests/test_tdvp.py
+++ b/tests/test_tdvp.py
@@ -1,5 +1,5 @@
 #!/usr/bin/python2
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 import numpy as np
 import pytest
 from tenpy.models.spins import SpinChain
diff --git a/tests/test_tebd.py b/tests/test_tebd.py
index b891e8b14..d0d1d3917 100644
--- a/tests/test_tebd.py
+++ b/tests/test_tebd.py
@@ -1,5 +1,5 @@
 """A collection of tests to check the functionality of `tenpy.tebd`"""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy.testing as npt
 import tenpy.linalg.np_conserved as npc
diff --git a/tests/test_terms.py b/tests/test_terms.py
index 6637ea57f..5e3cd7e4d 100644
--- a/tests/test_terms.py
+++ b/tests/test_terms.py
@@ -1,5 +1,5 @@
 """A collection of tests for :module:`tenpy.networks.terms`."""
-# Copyright 2019-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2019-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import copy
diff --git a/tests/test_time_evolution.py b/tests/test_time_evolution.py
index aac520eae..44f84d991 100644
--- a/tests/test_time_evolution.py
+++ b/tests/test_time_evolution.py
@@ -2,7 +2,7 @@
 
 Setup is a time evolution with TFI Model starting from a spin polarized state.
 """
-# Copyright 2020-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2020-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import scipy.linalg as LA
diff --git a/tests/test_tools.py b/tests/test_tools.py
index 6d1d22b18..8a7a98afe 100644
--- a/tests/test_tools.py
+++ b/tests/test_tools.py
@@ -1,5 +1,5 @@
 """A collection of tests for tenpy.tools submodules."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import logging
 import numpy as np
diff --git a/tests/test_truncation.py b/tests/test_truncation.py
index 2ec17ed7e..1180dfdcb 100644
--- a/tests/test_truncation.py
+++ b/tests/test_truncation.py
@@ -1,5 +1,5 @@
 """A test for tenpy.algorithms.truncation."""
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import numpy as np
 import pytest
diff --git a/tests/z_test_examples.py b/tests/z_test_examples.py
index 234402f0a..664778efb 100644
--- a/tests/z_test_examples.py
+++ b/tests/z_test_examples.py
@@ -3,7 +3,7 @@
 The files are only imported, so please protect example code from running with standard ``if
 __name__ == "__main__": ...`` clauses, if you want to demonstrate an interactive code.
 """
-# Copyright 2018-2021 TeNPy Developers, GNU GPLv3
+# Copyright 2018-2023 TeNPy Developers, GNU GPLv3
 
 import sys
 import os