diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
deleted file mode 100644
index 18ed248..0000000
--- a/.gitlab-ci.yml
+++ /dev/null
@@ -1,99 +0,0 @@
-before_script:
-  - which git || (apt-get update -qq && apt-get install --no-install-recommends -qqqy git)
-  - git config --global url."https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.com/".insteadOf "git@gitlab.com:"
-  - git config --global url."https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.com/".insteadOf "https://gitlab.com/" --add
-  - |
-      julia -e '
-      using Pkg
-      Pkg.Registry.add(
-        RegistrySpec(url = "https://github.com/JuliaRegistries/General.git")
-      )
-      Pkg.Registry.add(
-        RegistrySpec(url = "https://gitlab.com/tensors4fields/tensors4fieldsregistry.git")
-      )'
-
-variables:
-  CI_JULIA_CACHE_DIR: ${CI_PROJECT_DIR}/julia_pkg
-  JULIA_DEPOT_PATH: ${CI_JULIA_CACHE_DIR}
-cache:
-  key:
-    files:
-      - Project.toml
-      - docs/Project.toml
-    prefix: ${CI_JOB_NAME}
-  paths:
-    - ${CI_JULIA_CACHE_DIR}
-
-.script:
-  script:
-    - |
-      julia --project=@. -e '
-        using Pkg
-        Pkg.build()
-        Pkg.test(coverage=true)'
-.coverage:
-  coverage: /Test coverage (\d+\.\d+%)/
-  after_script:
-    - |
-      julia -e '
-        using Pkg
-        Pkg.add("Coverage")
-        using Coverage
-        c, t = get_summary(process_folder())
-        using Printf
-        @printf "Test coverage %.2f%%\n" 100c / t'
-Julia test:
-  image: julia:${JULIA_VERSION}
-  extends:
-    - .script
-    - .coverage
-  parallel:
-    matrix:
-      - JULIA_VERSION:
-        - "1.9"
-        - "1" # latest stable version
-.doctest:
-  script:
-    - |
-      julia --project=docs -e '
-        using Pkg
-        Pkg.develop(PackageSpec(path=pwd()))
-        Pkg.instantiate()
-        using Documenter: doctest
-        using Quantics
-        doctest(Quantics)
-        include("docs/make.jl")'
-doctest:
-  image: julia:1.6
-  extends:
-    - .doctest
-pages:
-  image: julia:1.6
-  stage: deploy
-  extends:
-    - .doctest
-  after_script:
-    - mkdir -p public
-    - mv docs/build public/dev
-  artifacts:
-    paths:
-      - public
-  only:
-    - main
-
-CompatHelper:
-  image: julia:1 # Set to the Julia version you want to use
-  rules:
-    - if: $CI_PIPELINE_SOURCE == "schedule"
-  script:
-    - |
-      julia --color=yes -e "
-        import Pkg;
-        name = \"CompatHelper\";
-        uuid = \"aa819f21-2bde-4658-8897-bab36330d9b7\";
-        version = \"3\";
-        Pkg.add(; name, uuid, version)"
-    - |
-      julia --color=yes -e "
-        import CompatHelper;
-        CompatHelper.main(;use_existing_registries=true)"
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..83beb7a
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,18 @@
+[submodule "TensorCrossInterpolation.jl"]
+	path = TensorCrossInterpolation.jl
+	url = git@github.com:tensor4all/TensorCrossInterpolation.jl.git
+[submodule "QuanticsGrids.jl"]
+	path = QuanticsGrids.jl
+	url = git@github.com:tensor4all/QuanticsGrids.jl.git
+[submodule "QuanticsTCI.jl"]
+	path = QuanticsTCI.jl
+	url = git@github.com:tensor4all/QuanticsTCI.jl.git
+[submodule "FastMPOContractions.jl"]
+	path = FastMPOContractions.jl
+	url = git@github.com:tensor4all/FastMPOContractions.jl.git
+[submodule "T4AITensorCompat.jl"]
+	path = T4AITensorCompat.jl
+	url = git@github.com:tensor4all/T4AITensorCompat.jl.git
+[submodule "PartitionedMPSs.jl"]
+	path = PartitionedMPSs.jl
+	url = git@github.com:tensor4all/PartitionedMPSs.jl.git
diff --git a/AGENTS.md b/AGENTS.md
new file mode 100644
index 0000000..28c9b03
--- /dev/null
+++ b/AGENTS.md
@@ -0,0 +1,33 @@
+- Use the same language as in past conversations with the user (if it has been Japanese, use Japanese)
+
+- All source code and documentation must be in English
+
+- Each subdirectory is a submodule. If there is an AGENTS.md in each directory, read it when working on the corresponding library
+
+- Each submodule is an independent Julia package with its own Project.toml, src/, test/, and docs/ directories. Understand the package structure before making changes
+
+- When working on a submodule, navigate into its directory and work as if it were a standalone package. Be aware of dependencies between packages
+
+- To update submodules after pulling changes, use `git submodule update --init --recursive`
+
+- Run tests for a specific package by navigating to its directory and running `julia --project=. -e 'using Pkg; Pkg.test()'` or `julia --project=. test/runtests.jl`
+
+- For debugging specific tests, it's more efficient to run only the relevant test file. For packages using `include()` statements in `runtests.jl` (like TensorCrossInterpolation.jl), comment out unnecessary includes and run only the test file you're debugging
+
+- Some libraries use ReTestItems as their test framework (e.g., Quantics.jl, QuanticsGrids.jl, TreeTCI.jl, SimpleTensorTrains.jl). However, ReTestItems has compatibility issues with libraries that use Distributed for parallel computation, so those libraries use the standard Test.jl framework instead
+
+- If a package has a `.JuliaFormatter.toml` file, follow its formatting rules. Otherwise, follow standard Julia style guidelines
+
+- When making changes that affect multiple packages, consider the dependency graph and test affected packages accordingly
+
+- The `gh` (GitHub CLI) command is available locally and can be used for GitHub-related operations
+
+- **Never push directly to main branch**: All changes must be made through pull requests. Create a branch, commit changes, push the branch, and create a PR. Wait for CI workflows to pass before merging.
+
+- **Never use force push to main branch**: Force pushing (including `--force-with-lease`) to main is prohibited. If you need to rewrite history, do it on a feature branch and create a PR.
+
+- All libraries are under the [tensor4all GitHub organization](https://github.com/tensor4all)
+
+- Some libraries are registered in T4ARegistry. Use T4ARegistrator.jl to register them. T4ARegistrator.jl is a development tool that should be installed in the global environment, not added as a dependency in individual package Project.toml files. When manually registering packages in T4ARegistry, use HTTPS URLs (not SSH) in the `repo` field of Package.toml to ensure compatibility in environments without SSH access
+
+- Some libraries are already registered in the official Julia registry. To register a new version, comment `@JuliaRegistrator register` in the library's issue, and the bot will create a PR to the official registry
diff --git a/Project.toml b/Project.toml
index 8c649e2..7845a99 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,34 +1,29 @@
-name = "Quantics"
-uuid = "87f76fb3-a40a-40c9-a63c-29fcfe7b7547"
+name = "T4AQuantics"
+uuid = "99202d80-6772-4c79-995b-f2660cf2bd6d"
 authors = ["Hiroshi Shinaoka <h.shinaoka@gmail.com> and contributors"]
 version = "0.4.6"
 
 [deps]
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
-FastMPOContractions = "f6e391d2-8ffa-4d7a-98cd-7e70024481ca"
-ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuanticsTCI = "b11687fd-3a1c-4c41-97d0-998ab401d50e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseIR = "4fe2279e-80f0-4adb-8463-ee114ff56b7d"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
-[weakdeps]
-PartitionedMPSs = "17ce1de9-5131-45e3-8a48-9723b6e2dc23"
+T4AITensorCompat = "864b39ca-388c-4f43-a593-a1076cf4b253"
+T4APartitionedMPSs = "f7a3e09c-877e-4116-af24-45397c1413c2"
 
 [extensions]
-QuantiscPartitionedMPSsExt = "PartitionedMPSs"
+QuantiscPartitionedMPSsExt = "T4APartitionedMPSs"
 
 [compat]
 EllipsisNotation = "1"
-FastMPOContractions = "^0.2.2"
-ITensorMPS = "0.3.1"
 ITensors = "0.7, 0.8, 0.9"
-PartitionedMPSs = "0.5.2, 0.6"
 QuanticsTCI = "0.7"
 SparseIR = "^0.96, 0.97, 1"
 StaticArrays = "1"
+T4APartitionedMPSs = "0.7"
 julia = "1"
 
 [extras]
@@ -38,8 +33,9 @@ QuanticsTCI = "b11687fd-3a1c-4c41-97d0-998ab401d50e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 SparseIR = "4fe2279e-80f0-4adb-8463-ee114ff56b7d"
+T4APartitionedMPSs = "f7a3e09c-877e-4116-af24-45397c1413c2"
 TensorCrossInterpolation = "b261b2ec-6378-4871-b32e-9173bb050604"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "ReTestItems", "Aqua", "PartitionedMPSs", "TensorCrossInterpolation", "QuanticsTCI", "QuanticsGrids", "SparseIR"]
+test = ["Test", "Random", "ReTestItems", "Aqua", "T4APartitionedMPSs", "TensorCrossInterpolation", "QuanticsTCI", "QuanticsGrids", "SparseIR"]
diff --git a/README.md b/README.md
index b1b6137..96959c6 100644
--- a/README.md
+++ b/README.md
@@ -1,19 +1,19 @@
-# Quantics
+# T4AQuantics
 
-[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://tensor4all.github.io/Quantics.jl/dev)
-[![CI](https://github.com/tensor4all/Quantics.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/tensor4all/Quantics.jl/actions/workflows/CI.yml)
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://tensor4all.github.io/T4AQuantics.jl/dev)
+[![CI](https://github.com/tensor4all/T4AQuantics.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/tensor4all/T4AQuantics.jl/actions/workflows/CI.yml)
 
 This library provides a high-level interface to manipulate quantics tensor train (QTT) format in Julia such as Fourier transform, convolution, and matrix-vector multiplication.
 This library is based on `ITensors.jl`.
 
 ## Installation
 
-The following will install `Quantics.jl`:
+The following will install `T4AQuantics.jl`:
 
 ```julia
-julia> using Pkg; Pkg.add("Quantics.jl")
+julia> using Pkg; Pkg.add("T4AQuantics.jl")
 ```
 
 ## Usage
 
-Please refer to the [documentation](https://tensor4all.github.io/Quantics.jl/) for usage.
\ No newline at end of file
+Please refer to the [documentation](https://tensor4all.github.io/T4AQuantics.jl/) for usage.
diff --git a/docs/make.jl b/docs/make.jl
index 1cb651b..fda66ab 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,7 +1,7 @@
-using Quantics
+using T4AQuantics
 using Documenter
 
-DocMeta.setdocmeta!(Quantics, :DocTestSetup, :(using Quantics); recursive=true)
+DocMeta.setdocmeta!(T4AQuantics, :DocTestSetup, :(using T4AQuantics); recursive=true)
 
 makedocs(;
     modules=[Quantics],
diff --git a/ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl b/ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl
index fd067e4..a31961a 100644
--- a/ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl
+++ b/ext/QuantiscPartitionedMPSsExt/QuantiscPartitionedMPSsExt.jl
@@ -4,26 +4,26 @@ using ITensors
 import ITensors
 using ITensors.SiteTypes: siteinds
 import ITensors.NDTensors: Tensor, BlockSparseTensor, blockview
-using ITensorMPS: ITensorMPS, MPS, MPO, AbstractMPS
-using ITensorMPS: findsite, linkinds, linkind, findsites
+using T4AITensorCompat: TensorTrain
+using T4AITensorCompat: findsite, linkinds, linkind, findsites
 
-import Quantics: Quantics, _find_site_allplevs, combinesites, extractdiagonal, _asdiagonal
-import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat,
+import T4AQuantics: T4AQuantics, _find_site_allplevs, combinesites, extractdiagonal, _asdiagonal
+import T4APartitionedMPSs: T4APartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat,
                         project
 
-function Quantics.makesitediagonal(subdmps::SubDomainMPS, site::Index)
+function T4AQuantics.makesitediagonal(subdmps::SubDomainMPS, site::Index)
     return _makesitediagonal(subdmps, site; baseplev=0)
 end
 
-function Quantics.makesitediagonal(subdmps::SubDomainMPS, sites::AbstractVector{Index})
+function T4AQuantics.makesitediagonal(subdmps::SubDomainMPS, sites::AbstractVector{Index})
     return _makesitediagonal(subdmps, sites; baseplev=0)
 end
 
-function Quantics.makesitediagonal(subdmps::SubDomainMPS, tag::String)
-    mps_diagonal = Quantics.makesitediagonal(MPS(subdmps), tag)
+function T4AQuantics.makesitediagonal(subdmps::SubDomainMPS, tag::String)
+    mps_diagonal = T4AQuantics.makesitediagonal(TensorTrain(subdmps), tag)
     subdmps_diagonal = SubDomainMPS(mps_diagonal)
 
-    target_sites = Quantics.findallsiteinds_by_tag(
+    target_sites = T4AQuantics.findallsiteinds_by_tag(
         unique(ITensors.noprime.(Iterators.flatten(siteinds(subdmps)))); tag=tag
     )
 
@@ -40,7 +40,7 @@ end
 function _makesitediagonal(
         subdmps::SubDomainMPS, sites::AbstractVector{Index{IndsT}}; baseplev=0
 ) where {IndsT}
-    M_ = deepcopy(MPO(collect(MPS(subdmps))))
+    M_ = deepcopy(TensorTrain(subdmps))
     for site in sites
         target_site::Int = only(findsites(M_, site))
         M_[target_site] = _asdiagonal(M_[target_site], site; baseplev=baseplev)
@@ -52,7 +52,7 @@ function _makesitediagonal(subdmps::SubDomainMPS, site::Index; baseplev=0)
     return _makesitediagonal(subdmps, [site]; baseplev=baseplev)
 end
 
-function Quantics.extractdiagonal(
+function T4AQuantics.extractdiagonal(
         subdmps::SubDomainMPS, sites::AbstractVector{Index{IndsT}}
 ) where {IndsT}
     tensors = collect(subdmps.data)
@@ -60,7 +60,7 @@ function Quantics.extractdiagonal(
         for site in intersect(sites, ITensors.inds(tensors[i]))
             sitewithallplevs = _find_site_allplevs(tensors[i], site)
             tensors[i] = if length(sitewithallplevs) > 1
-                tensors[i] = Quantics._extract_diagonal(tensors[i], sitewithallplevs...)
+                tensors[i] = T4AQuantics._extract_diagonal(tensors[i], sitewithallplevs...)
             else
                 tensors[i]
             end
@@ -73,32 +73,32 @@ function Quantics.extractdiagonal(
             delete!(projector.data, site')
         end
     end
-    return SubDomainMPS(MPS(tensors), projector)
+    return SubDomainMPS(TensorTrain(Vector{ITensor}(tensors)), projector)
 end
 
-function Quantics.extractdiagonal(subdmps::SubDomainMPS, site::Index{IndsT}) where {IndsT}
-    return Quantics.extractdiagonal(subdmps, [site])
+function T4AQuantics.extractdiagonal(subdmps::SubDomainMPS, site::Index{IndsT}) where {IndsT}
+    return T4AQuantics.extractdiagonal(subdmps, [site])
 end
 
-function Quantics.extractdiagonal(subdmps::SubDomainMPS, tag::String)::SubDomainMPS
-    targetsites = Quantics.findallsiteinds_by_tag(
-        unique(ITensors.noprime.(PartitionedMPSs._allsites(subdmps))); tag=tag
+function T4AQuantics.extractdiagonal(subdmps::SubDomainMPS, tag::String)::SubDomainMPS
+    targetsites = T4AQuantics.findallsiteinds_by_tag(
+        unique(ITensors.noprime.(T4APartitionedMPSs._allsites(subdmps))); tag=tag
     )
-    return Quantics.extractdiagonal(subdmps, targetsites)
+    return T4AQuantics.extractdiagonal(subdmps, targetsites)
 end
 
-function Quantics.rearrange_siteinds(subdmps::SubDomainMPS, sites)
-    return PartitionedMPSs.rearrange_siteinds(subdmps, sites)
+function T4AQuantics.rearrange_siteinds(subdmps::SubDomainMPS, sites)
+    return T4APartitionedMPSs.rearrange_siteinds(subdmps, sites)
 end
 
-function Quantics.rearrange_siteinds(partmps::PartitionedMPS, sites)
-    return PartitionedMPSs.rearrange_siteinds(partmps, sites)
+function T4AQuantics.rearrange_siteinds(partmps::PartitionedMPS, sites)
+    return T4APartitionedMPSs.rearrange_siteinds(partmps, sites)
 end
 
 """
 Make the PartitionedMPS diagonal for a given site index `s` by introducing a dummy index `s'`.
 """
-function Quantics.makesitediagonal(obj::PartitionedMPS, site)
+function T4AQuantics.makesitediagonal(obj::PartitionedMPS, site)
     return PartitionedMPS([_makesitediagonal(prjmps, site; baseplev=baseplev)
                            for prjmps in values(obj)])
 end
@@ -112,14 +112,14 @@ end
 Extract diagonal of the PartitionedMPS for `s`, `s'`, ... for a given site index `s`,
 where `s` must have a prime level of 0.
 """
-function Quantics.extractdiagonal(obj::PartitionedMPS, site)
+function T4AQuantics.extractdiagonal(obj::PartitionedMPS, site)
     return PartitionedMPS([extractdiagonal(prjmps, site) for prjmps in values(obj)])
 end
 
 """
 By default, elementwise multiplication will be performed.
 """
-function Quantics.automul(
+function T4AQuantics.automul(
         M1::PartitionedMPS,
         M2::PartitionedMPS;
         tag_row::String="",
@@ -146,16 +146,16 @@ function Quantics.automul(
     M1 = _makesitediagonal(M1, sites1_ewmul; baseplev=1)
     M2 = _makesitediagonal(M2, sites2_ewmul; baseplev=0)
 
-    sites_M1_diag = [collect(x) for x in siteinds(M1)]
-    sites_M2_diag = [collect(x) for x in siteinds(M2)]
+    sites_M1_diag = Vector{Vector{Index{Int}}}([collect(x) for x in siteinds(M1)])
+    sites_M2_diag = Vector{Vector{Index{Int}}}([collect(x) for x in siteinds(M2)])
 
-    M1 = Quantics.rearrange_siteinds(
-        M1, combinesites(sites_M1_diag, sites_row, sites_shared))
+    M1 = T4AQuantics.rearrange_siteinds(
+        M1, combinesites(sites_M1_diag, Vector{Index{Int}}(sites_row), Vector{Index{Int}}(sites_shared)))
 
-    M2 = Quantics.rearrange_siteinds(
-        M2, combinesites(sites_M2_diag, sites_shared, sites_col))
+    M2 = T4AQuantics.rearrange_siteinds(
+        M2, combinesites(sites_M2_diag, Vector{Index{Int}}(sites_shared), Vector{Index{Int}}(sites_col)))
 
-    M = PartitionedMPSs.contract(M1, M2; alg=alg, kwargs...)
+    M = T4APartitionedMPSs.contract(M1, M2; alg=alg, kwargs...)
 
     M = extractdiagonal(M, sites1_ewmul)
 
@@ -174,12 +174,12 @@ function Quantics.automul(
             end
         end
     end
-    return PartitionedMPSs.truncate(
-        Quantics.rearrange_siteinds(M, ressites); cutoff=cutoff, maxdim=maxdim)
+    return T4APartitionedMPSs.truncate(
+        T4AQuantics.rearrange_siteinds(M, ressites); cutoff=cutoff, maxdim=maxdim)
 end
 
 function _findallsiteinds_by_tag(M::PartitionedMPS; tag=tag)
-    return Quantics.findallsiteinds_by_tag(only.(siteinds(M)); tag=tag)
+    return T4AQuantics.findallsiteinds_by_tag(only.(siteinds(M)); tag=tag)
 end
 
 end
diff --git a/src/Quantics.jl b/src/T4AQuantics.jl
similarity index 62%
rename from src/Quantics.jl
rename to src/T4AQuantics.jl
index bb31360..a5d3001 100644
--- a/src/Quantics.jl
+++ b/src/T4AQuantics.jl
@@ -1,11 +1,12 @@
-module Quantics
+module T4AQuantics
 
 using ITensors
 import ITensors
-using ITensors.SiteTypes: siteinds
+using ITensors.SiteTypes: siteinds as ITensorsSiteTypes_siteinds
 import ITensors.NDTensors: Tensor, BlockSparseTensor, blockview
-using ITensorMPS: ITensorMPS, MPS, MPO, AbstractMPS
-using ITensorMPS: findsite, linkinds, linkind, findsites, truncate!
+import T4AITensorCompat
+using T4AITensorCompat: TensorTrain
+using T4AITensorCompat: findsite, linkinds, linkind, findsites, truncate!, siteinds, truncate
 
 import SparseIR: Fermionic, Bosonic, Statistics
 import LinearAlgebra: I, inv
@@ -13,7 +14,7 @@ using StaticArrays
 import SparseArrays: sparse
 
 import QuanticsTCI
-import FastMPOContractions
+using T4AITensorCompat: contract, Algorithm, apply, product
 
 using EllipsisNotation
 
diff --git a/src/affine.jl b/src/affine.jl
index c47f728..9f22663 100644
--- a/src/affine.jl
+++ b/src/affine.jl
@@ -49,7 +49,7 @@ function affine_transform_mpo(
         A::AbstractMatrix{<:Union{Integer,Rational}},
         b::AbstractVector{<:Union{Integer,Rational}},
         boundary::AbstractBoundaryConditions=PeriodicBoundaryConditions()
-)::MPO
+)::TensorTrain
     R = size(y, 1)
     M, N = size(A)
     size(x) == (R, N) ||
@@ -66,7 +66,7 @@ function affine_transform_mpo(
     link = [Index(size(tensors[r], 2); tags="link $r") for r in 1:(R - 1)]
 
     # Fill the MPO, taking care to not include auxiliary links at the edges
-    mpo = MPO(R)
+    mpo = TensorTrain(R)
     spin_dims = ntuple(_ -> 2, M + N)
     if R == 1
         mpo[1] = ITensor(reshape(tensors[1], spin_dims...),
@@ -321,7 +321,7 @@ end
 
 function affine_mpo_to_matrix(
         outsite::AbstractMatrix{<:Index}, insite::AbstractMatrix{<:Index},
-        mpo::MPO)
+        mpo::TensorTrain)
     prev_warn_order = ITensors.disable_warn_order()
     try
         mpo_contr = reduce(*, mpo)
diff --git a/src/binaryop.jl b/src/binaryop.jl
index 2cfa359..8a79eca 100644
--- a/src/binaryop.jl
+++ b/src/binaryop.jl
@@ -121,13 +121,15 @@ where a, b, c, d = +/- 1, 0, and s1, s1 are arbitrary integers.
 
 `bc` is a vector of boundary conditions for each arguments of `g` (not of `f`).
 """
-function affinetransform(M::MPS,
+function affinetransform(M::TensorTrain,
         tags::AbstractVector{String},
         coeffs_dic::AbstractVector{Dict{String,Int}},
         shift::AbstractVector{Int},
         bc::AbstractVector{Int};
         kwargs...)
-    transformer = affinetransformmpo(siteinds(M), tags, coeffs_dic, shift, bc)
+    sites_flat = collect(Iterators.flatten(siteinds(M)))
+    sites_vec = Vector{Index{Int}}(sites_flat)
+    transformer = affinetransformmpo(sites_vec, tags, coeffs_dic, shift, bc)
     return _apply(transformer, M; kwargs...)
 end
 
@@ -135,13 +137,13 @@ function affinetransformmpo(sites::AbstractVector{Index{T}},
         tags::AbstractVector{String},
         coeffs_dic::AbstractVector{Dict{String,Int}},
         shift::AbstractVector{Int},
-        bc::AbstractVector{Int})::MPO where {T}
+        bc::AbstractVector{Int})::TensorTrain where {T}
     # f(x, y) = g(a * x + b * y + s1, c * x + d * y + s2)
     #         = h(a * x + b * y,      c * x + d * y),
     # where h(x, y) = g(x + s1, y + s2).
     # The transformation is executed in this order: g -> h -> f.
 
-    mpos = MPO[]
+    mpos = TensorTrain[]
 
     # Number of variables involved in transformation
     ntransvars = length(tags)
@@ -185,12 +187,14 @@ end
 Affine transform of a MPS with no shift
 Significant bits are assumed to be aligned from left to right for all tags.
 """
-function affinetransform(M::MPS,
+function affinetransform(M::TensorTrain,
         tags::AbstractVector{String},
         coeffs_dic::AbstractVector{Dict{String,Int}},
         bc::AbstractVector{Int};
         kwargs...)
-    transformer = affinetransformmpo(siteinds(M), tags, coeffs_dic, bc)
+    sites_flat = collect(Iterators.flatten(siteinds(M)))
+    sites_vec = Vector{Index{Int}}(sites_flat)
+    transformer = affinetransformmpo(sites_vec, tags, coeffs_dic, bc)
     return _apply(transformer, M; kwargs...)
 end
 
@@ -201,8 +205,8 @@ Significant bits are assumed to be aligned from left to right for all tags.
 function affinetransformmpo(sites::AbstractVector{Index{T}},
         tags::AbstractVector{String},
         coeffs_dic::AbstractVector{Dict{String,Int}},
-        bc::AbstractVector{Int})::MPO where {T}
-    mpos = MPO[]
+        bc::AbstractVector{Int})::TensorTrain where {T}
+    mpos = TensorTrain[]
 
     # f(x, y) = g(a * x + b * y + s1, c * x + d * y + s2)
     #         = h(a * x + b * y,      c * x + d * y),
@@ -275,7 +279,7 @@ function affinetransformmpo(sites::AbstractVector{Index{T}},
     if !rev_carrydirec
         transformer_ = affinetransformmpo(
             reverse(sites), reverse(tags), reverse(coeffs_dic), reverse(bc))
-        return MPO([transformer_[n] for n in reverse(1:length(transformer_))])
+        return TensorTrain([transformer_[n] for n in reverse(1:length(transformer_))])
     end
 
     # First check transformations with -1 and -1; e.g., (a, b) = (-1, -1)
@@ -426,7 +430,7 @@ function _binaryop_mpo_backend(sites::Vector{Index{T}},
 
     _removeedges!(tensors, sites)
 
-    M = truncate(MPO(tensors); cutoff=1e-25)
+    M = truncate(TensorTrain(tensors); cutoff=1e-25)
     cleanup_linkinds!(M)
 
     return M
@@ -443,7 +447,7 @@ function _shift_mpo(sites::Vector{Index{T}}, shift::Int; bc::Int=1) where {T<:Nu
     R = length(sites)
     0 <= shift <= 2^R - 1 || error("Invalid shift")
 
-    ys = Quantics.tobin(shift, R)
+    ys = T4AQuantics.tobin(shift, R)
 
     links = Index{T}[]
     tensors = ITensor[]
@@ -452,7 +456,7 @@ function _shift_mpo(sites::Vector{Index{T}}, shift::Int; bc::Int=1) where {T<:Nu
         cin_on = n != R
         cout_on = n != 1
         sitey = Index(2, "Qubit, y")
-        t, link_in, link_out = Quantics._binaryop_tensor(1, 1, sites[n]', sitey, sites[n],
+        t, link_in, link_out = T4AQuantics._binaryop_tensor(1, 1, sites[n]', sitey, sites[n],
             cin_on, cout_on, bc)
         t *= onehot(sitey => ys[n] + 1)
         if n < R
@@ -471,5 +475,5 @@ function _shift_mpo(sites::Vector{Index{T}}, shift::Int; bc::Int=1) where {T<:Nu
         push!(tensors, t)
     end
 
-    MPO(tensors)
+    TensorTrain(tensors)
 end
diff --git a/src/fouriertransform.jl b/src/fouriertransform.jl
index 40540f1..7e9aac3 100644
--- a/src/fouriertransform.jl
+++ b/src/fouriertransform.jl
@@ -29,7 +29,7 @@ function _qft(sites; cutoff::Float64=1e-25, sign::Int=1)
 
     sites_MPO = collect.(zip(prime.(sites), sites))
     fouriertt = QuanticsTCI.quanticsfouriermpo(R; sign=Float64(sign), normalize=true)
-    M = MPO(fouriertt; sites=sites_MPO)
+    M = TensorTrain(fouriertt; sites=sites_MPO)
 
     return truncate(M; cutoff)
 end
@@ -37,7 +37,7 @@ end
 abstract type AbstractFT end
 
 struct FTCore
-    forward::MPO
+    forward::TensorTrain
 
     function FTCore(sites; kwargs...)
         new(_qft(sites; kwargs...))
@@ -62,12 +62,12 @@ function forwardmpo(ftcore::FTCore, sites)
 end
 
 function backwardmpo(ftcore::FTCore, sites)
-    M = conj(MPO(reverse([x for x in ftcore.forward])))
+    M = conj(TensorTrain(reverse([x for x in ftcore.forward])))
     _replace_mpo_siteinds!(M, _extractsites(M), sites)
     return M
 end
 
-function _apply_qft(M::MPO, gsrc::MPS, target_sites, sitepos, sitesdst; kwargs...)
+function _apply_qft(M::TensorTrain, gsrc::TensorTrain, target_sites, sitepos, sitesdst; kwargs...)
     _replace_mpo_siteinds!(M, _extractsites(M), target_sites)
     M = matchsiteinds(M, siteinds(gsrc))
     gdst = _apply(M, gsrc; kwargs...)
@@ -103,7 +103,7 @@ where ``s = \pm 1``, ``x_0`` and ``y_0`` are constants, ``N=2^R``.
 Instead of specifying `sitessrc`, one can specify the source sites by setting `tag`.
 If `tag` = `x`, all sites with tags `x=1`, `x=2`, ... are used as `sitessrc`.
 """
-function fouriertransform(M::MPS;
+function fouriertransform(M::TensorTrain;
         sign::Int=1,
         tag::String="",
         sitessrc=nothing,
@@ -111,7 +111,8 @@ function fouriertransform(M::MPS;
         originsrc::Real=0.0,
         origindst::Real=0.0,
         cutoff_MPO=1e-25, kwargs...)
-    sites = siteinds(M)
+    sites_flat = collect(Iterators.flatten(siteinds(M)))
+    sites_vec = Vector{Index{Int}}(sites_flat)
     sitepos, target_sites = _find_target_sites(M; sitessrc=sitessrc, tag=tag)
 
     if sitesdst === nothing
@@ -127,14 +128,14 @@ function fouriertransform(M::MPS;
 
     # Prepare MPO for QFT
     MQ_ = _qft(target_sites; sign=sign, cutoff=cutoff_MPO)
-    MQ = matchsiteinds(MQ_, sites)
+    MQ = matchsiteinds(MQ_, sites_vec)
 
     # Phase shift from origindst
     M_result = phase_rotation(M, sign * 2π * BigFloat(origindst) / (BigFloat(2)^length(sitepos));
         targetsites=target_sites, kwargs...)
 
     # Apply QFT
-    M_result = _apply(MQ, M_result; kwargs...)
+    M_result = _apply(MQ, M_result; alg="naive", kwargs...)
 
     N = length(target_sites)
     for n in eachindex(target_sites)
diff --git a/src/imaginarytime.jl b/src/imaginarytime.jl
index 5a2257c..114b6de 100644
--- a/src/imaginarytime.jl
+++ b/src/imaginarytime.jl
@@ -16,8 +16,8 @@ _stat_shift(::Bosonic) = 0
 _stat_sign(::Fermionic) = -1
 _stat_sign(::Bosonic) = 1
 
-function to_wn(stat::Statistics, gtau::MPS, beta::Float64; sitessrc=nothing, tag="",
-        sitesdst=nothing, kwargs...)::MPS
+function to_wn(stat::Statistics, gtau::TensorTrain, beta::Float64; sitessrc=nothing, tag="",
+        sitesdst=nothing, kwargs...)::TensorTrain
     sitepos, _ = _find_target_sites(gtau; sitessrc=sitessrc, tag=tag)
     nqbit_t = length(sitepos)
     p_back = precision(BigFloat)
@@ -30,8 +30,8 @@ function to_wn(stat::Statistics, gtau::MPS, beta::Float64; sitessrc=nothing, tag
     return giv
 end
 
-function to_tau(stat::Statistics, giv::MPS, beta::Float64; sitessrc=nothing, tag="",
-        sitesdst=nothing, kwargs...)::MPS
+function to_tau(stat::Statistics, giv::TensorTrain, beta::Float64; sitessrc=nothing, tag="",
+        sitesdst=nothing, kwargs...)::TensorTrain
     sitepos, _ = _find_target_sites(giv; sitessrc=sitessrc, tag=tag)
     p_back = precision(BigFloat)
     setprecision(BigFloat, 256)
@@ -52,7 +52,7 @@ function decompose_gtau(gtau_smpl::Vector{ComplexF64}, sites; kwargs...)
     gtau_smpl = reshape(gtau_smpl, repeat([2], nbit)...)
     gtau_smpl = permutedims(gtau_smpl, reverse(collect(1:nbit)))
 
-    return MPS(gtau_smpl, sites; kwargs...)
+    return TensorTrain(gtau_smpl, sites; kwargs...)
 end
 
 """
@@ -66,7 +66,7 @@ function decompose_giv(giv_smpl::Vector{ComplexF64}, sites; kwargs...)
     nbit = length(sites)
     length(giv_smpl) == 2^nbit || error("Length mismatch")
     tensor = ITensor(giv_smpl, reverse(sites))
-    return MPS(tensor, reverse(sites); kwargs...)
+    return TensorTrain(tensor, reverse(sites); kwargs...)
 end
 
 """
@@ -83,7 +83,7 @@ function poletomps(stat::Statistics, sites, β, ω)
     tensors[1] *= -1 / (1 - _stat_sign(stat) * exp(-β * ω))
     tensors[1] *= onehot(links[1] => 1)
     tensors[end] *= onehot(links[end] => 1)
-    return MPS(tensors)
+    return TensorTrain(tensors)
 end
 
 """
@@ -93,14 +93,16 @@ function poletomps(::Fermionic, sites, β, ω)
     if β * ω < -100.0
         R = length(sites)
         sites_ = [Index(2, "Qubit,τ=$n") for n in 1:R]
-        tmp = reverseaxis(Quantics.expqtt(sites_, β * ω); tag="τ", bc=1)
+        tmp = reverseaxis(T4AQuantics.expqtt(sites_, β * ω); tag="τ", bc=1)
         res = shiftaxis(tmp, +1; tag="τ", bc=1) * (-exp(β * ω / 2^R))
         for n in 1:R
             replaceind!(res[n], sites_[n] => sites[n])
         end
         return res
     end
-    return expqtt(sites, -β * ω) / (-1 - exp(-β * ω))
+    result = expqtt(sites, -β * ω)
+    divisor = -1 - exp(-β * ω)
+    return result * (1.0 / divisor)
 end
 
 poletomps(sites, β, ω) = poletomps(Fermionic(), sites, β, ω)
diff --git a/src/mps.jl b/src/mps.jl
index d0fdc9b..5f84ba9 100644
--- a/src/mps.jl
+++ b/src/mps.jl
@@ -2,7 +2,7 @@
 Create a MPS filled with one
 """
 function onemps(::Type{T}, sites) where {T<:Number}
-    M = MPS(T, sites; linkdims=1)
+    M = TensorTrain(T, sites; linkdims=1)
     l = linkinds(M)
     for n in eachindex(M)
         if n == 1
@@ -32,5 +32,5 @@ function expqtt(sites, a::Float64)
     end
     tensors[1] *= onehot(links[1] => 1)
     tensors[end] *= onehot(links[end] => 1)
-    return MPS(tensors)
+    return TensorTrain(tensors)
 end
diff --git a/src/mul.jl b/src/mul.jl
index a6c75b5..7c761cf 100644
--- a/src/mul.jl
+++ b/src/mul.jl
@@ -21,7 +21,7 @@ function MatrixMultiplier(site_row::Index{T},
     return MatrixMultiplier([site_row], [site_shared], [site_col])
 end
 
-function preprocess(mul::MatrixMultiplier{T}, M1::MPO, M2::MPO) where {T}
+function preprocess(mul::MatrixMultiplier{T}, M1::TensorTrain, M2::TensorTrain) where {T}
     for (site_row, site_shared, site_col) in zip(mul.sites_row, mul.sites_shared,
         mul.sites_col)
         M1, M2 = combinesites(M1, site_row, site_shared),
@@ -30,21 +30,24 @@ function preprocess(mul::MatrixMultiplier{T}, M1::MPO, M2::MPO) where {T}
     return M1, M2
 end
 
-function postprocess(mul::MatrixMultiplier{T}, M::MPO)::MPO where {T}
-    tensors = ITensors.data(M)
+function postprocess(mul::MatrixMultiplier{T}, M::TensorTrain)::TensorTrain where {T}
+    tensors = collect(M)
     for (site_row, site_col) in zip(mul.sites_row, mul.sites_col)
         p = findfirst(hasind(site_row), tensors)
         hasind(tensors[p], site_col) ||
             error("$site_row and $site_col are not on the same site")
 
+        # Recompute links after each modification
+        links = length(tensors) > 1 ? [commoninds(tensors[n], tensors[n + 1])[1] for n in 1:(length(tensors) - 1)] : Index[]
+
         indsl = [site_row]
-        if p > 1
-            push!(indsl, linkind(M, p - 1))
+        if p > 1 && (p - 1) <= length(links)
+            push!(indsl, links[p - 1])
         end
 
         indsr = [site_col]
-        if p < length(M)
-            push!(indsr, linkind(M, p))
+        if p < length(tensors) && p <= length(links)
+            push!(indsr, links[p])
         end
 
         Ml, Mr = split_tensor(tensors[p], [indsl, indsr])
@@ -54,7 +57,7 @@ function postprocess(mul::MatrixMultiplier{T}, M::MPO)::MPO where {T}
         insert!(tensors, p + 1, Mr)
     end
 
-    return MPO(tensors)
+    return TensorTrain(Vector{ITensor}(tensors))
 end
 
 #===
@@ -103,9 +106,9 @@ function _todense(t, site::Index{T}) where {T<:Number}
     return ITensor(newdata, links..., site)
 end
 
-function preprocess(mul::ElementwiseMultiplier{T}, M1::MPO, M2::MPO) where {T}
-    tensors1 = ITensors.data(M1)
-    tensors2 = ITensors.data(M2)
+function preprocess(mul::ElementwiseMultiplier{T}, M1::TensorTrain, M2::TensorTrain) where {T}
+    tensors1 = collect(M1)
+    tensors2 = collect(M2)
     for s in mul.sites
         p = findfirst(hasind(s), tensors1)
         hasinds(tensors2[p], s) || error("ITensor of M2 at $p does not have $s")
@@ -115,52 +118,57 @@ function preprocess(mul::ElementwiseMultiplier{T}, M1::MPO, M2::MPO) where {T}
         replaceind!(tensors1[p], s => s')
         tensors2[p] = _asdiagonal(tensors2[p], s)
     end
-    return MPO(tensors1), MPO(tensors2)
+    return TensorTrain(Vector{ITensor}(tensors1)), TensorTrain(Vector{ITensor}(tensors2))
 end
 
-function postprocess(mul::ElementwiseMultiplier{T}, M::MPO)::MPO where {T}
-    tensors = ITensors.data(M)
+function postprocess(mul::ElementwiseMultiplier{T}, M::TensorTrain)::TensorTrain where {T}
+    tensors = collect(M)
     for s in mul.sites
         p = findfirst(hasind(s), tensors)
         tensors[p] = _todense(tensors[p], s)
     end
-    return MPO(tensors)
+    return TensorTrain(Vector{ITensor}(tensors))
 end
 
 """
 By default, elementwise multiplication will be performed.
 """
-function automul(M1::MPS, M2::MPS; tag_row::String="", tag_shared::String="",
+function automul(M1::TensorTrain, M2::TensorTrain; tag_row::String="", tag_shared::String="",
         tag_col::String="", alg="naive", cutoff=1e-30, kwargs...)
     if in(:maxbonddim, keys(kwargs))
         error("Illegal keyward parameter: maxbonddim. Use maxdim instead!")
     end
 
-    sites_row = findallsiteinds_by_tag(siteinds(M1); tag=tag_row)
-    sites_shared = findallsiteinds_by_tag(siteinds(M1); tag=tag_shared)
-    sites_col = findallsiteinds_by_tag(siteinds(M2); tag=tag_col)
+    sites1_flat = collect(Iterators.flatten(siteinds(M1)))
+    sites2_flat = collect(Iterators.flatten(siteinds(M2)))
+    sites1_vec = Vector{Index{Int}}(sites1_flat)
+    sites2_vec = Vector{Index{Int}}(sites2_flat)
+    
+    sites_row = findallsiteinds_by_tag(sites1_vec; tag=tag_row)
+    sites_shared = findallsiteinds_by_tag(sites1_vec; tag=tag_shared)
+    sites_col = findallsiteinds_by_tag(sites2_vec; tag=tag_col)
     sites_matmul = Set(Iterators.flatten([sites_row, sites_shared, sites_col]))
 
-    if sites_shared != findallsiteinds_by_tag(siteinds(M2); tag=tag_shared)
+    if sites_shared != findallsiteinds_by_tag(sites2_vec; tag=tag_shared)
         error("Invalid shared sites for MatrixMultiplier")
     end
 
     matmul = MatrixMultiplier(sites_row, sites_shared, sites_col)
-    ewmul = ElementwiseMultiplier([s for s in siteinds(M1) if s ∉ sites_matmul])
+    ewmul = ElementwiseMultiplier([s for s in sites1_vec if s ∉ sites_matmul])
 
-    M1_ = Quantics.asMPO(M1)
-    M2_ = Quantics.asMPO(M2)
-    M1_, M2_ = preprocess(matmul, M1_, M2_)
+    M1_, M2_ = preprocess(matmul, M1, M2)
     M1_, M2_ = preprocess(ewmul, M1_, M2_)
 
-    M = FastMPOContractions.contract_mpo_mpo(M1_, M2_; alg=alg, cutoff=cutoff, kwargs...)
+    # Convert alg to Algorithm type
+    alg_ = alg isa Algorithm ? alg : Algorithm(alg)
+    M = contract(M1_, M2_; alg=alg_, cutoff=cutoff, kwargs...)
 
-    M = Quantics.postprocess(matmul, M)
-    M = Quantics.postprocess(ewmul, M)
+    M = T4AQuantics.postprocess(matmul, M)
+    M = T4AQuantics.postprocess(ewmul, M)
 
     if in(:maxdim, keys(kwargs))
         truncate!(M; maxdim=kwargs[:maxdim], cutoff=cutoff)
     end
 
-    return asMPS(M)
+    return M
 end
diff --git a/src/transformer.jl b/src/transformer.jl
index 9d03d82..5150534 100644
--- a/src/transformer.jl
+++ b/src/transformer.jl
@@ -19,7 +19,7 @@ f(x) = g(-x)
 where f(x) = M * g(x) for x = 0, 1, ..., 2^R-1.
 """
 function flipop_to_negativedomain(sites::Vector{Index{T}}; rev_carrydirec=false,
-        bc::Int=1)::MPO where {T}
+        bc::Int=1)::TensorTrain where {T}
     return flipop(sites; rev_carrydirec=rev_carrydirec, bc=bc) * bc
 end
 
@@ -30,10 +30,10 @@ where f(x) = M * g(x) for x = 0, 1, ..., 2^R-1.
 
 `sites`: the sites of the output MPS
 """
-function flipop(sites::Vector{Index{T}}; rev_carrydirec=false, bc::Int=1)::MPO where {T}
+function flipop(sites::Vector{Index{T}}; rev_carrydirec=false, bc::Int=1)::TensorTrain where {T}
     if rev_carrydirec
         M = flipop(reverse(sites); rev_carrydirec=false, bc=bc)
-        return MPO([M[n] for n in reverse(1:length(M))])
+        return TensorTrain([M[n] for n in reverse(1:length(M))])
     end
 
     N = length(sites)
@@ -41,7 +41,7 @@ function flipop(sites::Vector{Index{T}}; rev_carrydirec=false, bc::Int=1)::MPO w
     N > 1 || error("MPO with one tensor is not supported")
 
     t = _single_tensor_flip()
-    M = MPO(N)
+    M = TensorTrain(N)
     links = [Index(2, "Link,l=$l") for l in 1:(N + 1)]
     for n in 1:N
         M[n] = ITensor(t, (links[n], links[n + 1], sites[n]', sites[n]))
@@ -63,12 +63,14 @@ where x = 0, 1, ..., N-1 and N = 2^R.
 
 Note that x = 0, 1, 2, ..., N-1 are mapped to x = 0, N-1, N-2, ..., 1 mod N.
 """
-function reverseaxis(M::MPS; tag="x", bc::Int=1, kwargs...)
+function reverseaxis(M::TensorTrain; tag="x", bc::Int=1, kwargs...)
     bc ∈ [1, -1] || error("bc must be either 1 or -1")
-    return _apply(reverseaxismpo(siteinds(M); tag=tag, bc=bc), M; kwargs...)
+    sites_flat = collect(Iterators.flatten(siteinds(M)))
+    sites_vec = Vector{Index{Int}}(sites_flat)
+    return _apply(reverseaxismpo(sites_vec; tag=tag, bc=bc), M; alg="naive", kwargs...)
 end
 
-function reverseaxismpo(sites::AbstractVector{Index{T}}; tag="x", bc::Int=1)::MPO where {T}
+function reverseaxismpo(sites::AbstractVector{Index{T}}; tag="x", bc::Int=1)::TensorTrain where {T}
     bc ∈ [1, -1] || error("bc must be either 1 or -1")
     targetsites = findallsiteinds_by_tag(sites; tag=tag)
     pos = findallsites_by_tag(sites; tag=tag)
@@ -82,16 +84,18 @@ end
 """
 f(x) = g(x + shift) for x = 0, 1, ..., 2^R-1 and 0 <= shift < 2^R.
 """
-function shiftaxis(M::MPS, shift::Int; tag="x", bc::Int=1, kwargs...)
+function shiftaxis(M::TensorTrain, shift::Int; tag="x", bc::Int=1, kwargs...)
     bc ∈ [1, -1] || error("bc must be either 1 or -1")
-    return _apply(shiftaxismpo(siteinds(M), shift; tag=tag, bc=bc), M; kwargs...)
+    sites_flat = collect(Iterators.flatten(siteinds(M)))
+    sites_vec = Vector{Index{Int}}(sites_flat)
+    return _apply(shiftaxismpo(sites_vec, shift; tag=tag, bc=bc), M; alg="naive", kwargs...)
 end
 
 """
 f(x) = g(x + shift) for x = 0, 1, ..., 2^R-1 and 0 <= shift < 2^R.
 """
 function shiftaxismpo(
-        sites::AbstractVector{Index{T}}, shift::Int; tag="x", bc::Int=1)::MPO where {T}
+        sites::AbstractVector{Index{T}}, shift::Int; tag="x", bc::Int=1)::TensorTrain where {T}
     bc ∈ [1, -1] || error("bc must be either 1 or -1")
     targetsites = findallsiteinds_by_tag(sites; tag=tag) # From left to right: x=1, 2, ...
     pos = findallsites_by_tag(sites; tag=tag)
@@ -105,7 +109,7 @@ function shiftaxismpo(
         transformer = _shift_mpo(targetsites, shift_mod; bc=bc)
     else
         transformer = _shift_mpo(targetsites, shift_mod; bc=bc)
-        transformer = MPO([transformer[n] for n in reverse(1:length(transformer))])
+        transformer = TensorTrain([transformer[n] for n in reverse(1:length(transformer))])
     end
     transformer = matchsiteinds(transformer, sites)
     transformer *= bc^nbc
@@ -116,9 +120,11 @@ end
 """
 Multiply by exp(i θ x), where x = (x_1, ..., x_R)_2.
 """
-function phase_rotation(M::MPS, θ::Real; targetsites=nothing, tag="", kwargs...)::MPS
-    transformer = phase_rotation_mpo(siteinds(M), θ; targetsites=targetsites, tag=tag)
-    _apply(transformer, M; kwargs...)
+function phase_rotation(M::TensorTrain, θ::Real; targetsites=nothing, tag="", kwargs...)::TensorTrain
+    sites_flat = collect(Iterators.flatten(siteinds(M)))
+    sites_vec = Vector{Index{Int}}(sites_flat)
+    transformer = phase_rotation_mpo(sites_vec, θ; targetsites=targetsites, tag=tag)
+    _apply(transformer, M; alg="naive", kwargs...)
 end
 
 """
@@ -127,13 +133,13 @@ Create an MPO for multiplication by `exp(i θ x)`, where `x = (x_1, ..., x_R)_2`
 `sites`: site indices for `x_1`, `x_2`, ..., `x_R`.
 """
 function phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Real;
-        targetsites=nothing, tag="")::MPO where {T}
+        targetsites=nothing, tag="")::TensorTrain where {T}
     _, target_sites = _find_target_sites(sites; sitessrc=targetsites, tag=tag)
     transformer = _phase_rotation_mpo(target_sites, θ)
     return matchsiteinds(transformer, sites)
 end
 
-function _phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Real)::MPO where {T}
+function _phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Real)::TensorTrain where {T}
     R = length(sites)
     tensors = [ITensor(true) for _ in 1:R]
     θ_mod = mod(θ, 2π)
@@ -154,7 +160,7 @@ function _phase_rotation_mpo(sites::AbstractVector{Index{T}}, θ::Real)::MPO whe
         Array(tensors[end], sites[end]', sites[end]), links[end], sites[end], sites[end]')
 
     setprecision(BigFloat, p_back)
-    return MPO(tensors)
+    return TensorTrain(tensors)
 end
 
 function _upper_lower_triangle(upper_or_lower::Symbol)::Array{Float64,4}
@@ -183,13 +189,13 @@ end
 Create QTT for a upper/lower triangle matrix filled with one except the diagonal line
 """
 function upper_lower_triangle_matrix(sites::Vector{Index{T}}, value::S;
-        upper_or_lower::Symbol=:upper)::MPO where {T,S}
+        upper_or_lower::Symbol=:upper)::TensorTrain where {T,S}
     upper_or_lower ∈ [:upper, :lower] || error("Invalid upper_or_lower $(upper_or_lower)")
     N = length(sites)
 
     t = _upper_lower_triangle(upper_or_lower)
 
-    M = MPO(N)
+    M = TensorTrain(N)
     links = [Index(2, "Link,l=$l") for l in 1:(N + 1)]
     for n in 1:N
         M[n] = ITensor(t, (links[n], links[n + 1], sites[n]', sites[n]))
@@ -216,12 +222,12 @@ end
 Add new site indices to an MPS
 """
 #==
-function asdiagonal(M::MPS, newsites; which_new="right", targetsites=nothing, tag="")
+function asdiagonal(M::TensorTrain, newsites; which_new="right", targetsites=nothing, tag="")
     which_new ∈ ["left", "right"] || error("Invalid which_new: left or right")
-    sitepos, target_sites = Quantics._find_target_sites(M; sitessrc=targetsites, tag=tag)
+    sitepos, target_sites = T4AQuantics._find_target_sites(M; sitessrc=targetsites, tag=tag)
     length(sitepos) == length(newsites) ||
         error("Length mismatch: $(newsites) vs $(target_sites)")
-    M_ = Quantics._addedges(M)
+    M_ = T4AQuantics._addedges(M)
     links = linkinds(M_)
 
     tensors = ITensor[]
@@ -247,8 +253,8 @@ function asdiagonal(M::MPS, newsites; which_new="right", targetsites=nothing, ta
     tensors[1] *= onehot(links[1] => 1)
     tensors[end] *= onehot(links[end] => 1)
 
-    M_result = MPS(tensors)
-    Quantics.cleanup_linkinds!(M_result)
+    M_result = TensorTrain(tensors)
+    T4AQuantics.cleanup_linkinds!(M_result)
     return M_result
 end
 ==#
diff --git a/src/util.jl b/src/util.jl
index 3da70d1..95a1695 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -1,4 +1,4 @@
-function _extractsite(x::Union{MPS,MPO}, n::Int)
+function _extractsite(x::TensorTrain, n::Int)
     if n == 1
         return noprime(copy(uniqueind(x[n], x[n + 1])))
     elseif n == length(x)
@@ -8,9 +8,9 @@ function _extractsite(x::Union{MPS,MPO}, n::Int)
     end
 end
 
-_extractsites(x::Union{MPS,MPO}) = [_extractsite(x, n) for n in eachindex(x)]
+_extractsites(x::TensorTrain) = [_extractsite(x, n) for n in eachindex(x)]
 
-function _replace_mpo_siteinds!(M::MPO, sites_src, sites_dst)
+function _replace_mpo_siteinds!(M::TensorTrain, sites_src, sites_dst)
     sites_src = noprime(sites_src)
     sites_dst = noprime(sites_dst)
     for j in eachindex(M)
@@ -24,10 +24,10 @@ end
 Reverse the order of the MPS/MPO tensors
 The order of the siteinds are reversed in the returned object.
 """
-function _reverse(M::MPO)
+function _reverse(M::TensorTrain)
     sites = _extractsites(M)
     N = length(M)
-    M_ = MPO([M[n] for n in reverse(1:length(M))])
+    M_ = TensorTrain([M[n] for n in reverse(1:length(M))])
     for n in 1:N
         replaceind!(M_[n], sites[N - n + 1], sites[n])
         replaceind!(M_[n], sites[N - n + 1]', sites[n]')
@@ -41,7 +41,7 @@ Create a MPO with ITensor objects of ElType ComplexF64 filled with zero
 function _zero_mpo(sites; linkdims=ones(Int, length(sites) - 1))
     length(linkdims) == length(sites) - 1 ||
         error("Length mismatch $(length(linkdims)) != $(length(sites)) - 1")
-    M = MPO(sites)
+    M = TensorTrain(sites)
     N = length(M)
     links = [Index(1, "n=0,Link")]
     for n in 1:(N - 1)
@@ -92,8 +92,8 @@ new_sites: Vector of vectors of new siteinds
 
 When splitting MPS tensors, the column major is assumed.
 """
-function unfuse_siteinds(M::MPS, targetsites::Vector{Index{T}},
-        newsites::AbstractVector{Vector{Index{T}}})::MPS where {T}
+function unfuse_siteinds(M::TensorTrain, targetsites::Vector{Index{T}},
+        newsites::AbstractVector{Vector{Index{T}}})::TensorTrain where {T}
     length(targetsites) == length(newsites) || error("Length mismatch")
     links = linkinds(M)
     L = length(M)
@@ -129,26 +129,29 @@ function unfuse_siteinds(M::MPS, targetsites::Vector{Index{T}},
         end
     end
 
-    M_ = MPS(tensors_)
+    M_ = TensorTrain(tensors_)
     cleanup_linkinds!(M_)
     return M_
 end
 
-function _removeedges!(x::MPS, sites)
-    length(inds(x[1])) == 3 || error("Dim of the first tensor must be 3")
-    length(inds(x[end])) == 3 || error("Dim of the last tensor must be 3")
-    elt = eltype(x[1])
-    x[1] *= onehot(elt, uniqueind(x[1], x[2], sites) => 1)
-    x[end] *= onehot(elt, uniqueind(x[end], x[end - 1], sites) => 1)
-    return nothing
-end
-
-function _removeedges!(x::MPO, sites)
-    length(inds(x[1])) == 4 || error("Dim of the first tensor must be 4")
-    length(inds(x[end])) == 4 || error("Dim of the last tensor must be 4")
-    elt = eltype(x[1])
-    x[1] *= onehot(elt, uniqueind(x[1], x[2], sites, prime.(sites)) => 1)
-    x[end] *= onehot(elt, uniqueind(x[end], x[end - 1], sites, prime.(sites)) => 1)
+function _removeedges!(x::TensorTrain, sites)
+    # Check if it's MPS-like (3 indices) or MPO-like (4 indices)
+    n_inds_first = length(inds(x[1]))
+    n_inds_last = length(inds(x[end]))
+    
+    if n_inds_first == 3 && n_inds_last == 3
+        # MPS-like: 3 indices (link, site, link)
+        elt = eltype(x[1])
+        x[1] *= onehot(elt, uniqueind(x[1], x[2], sites) => 1)
+        x[end] *= onehot(elt, uniqueind(x[end], x[end - 1], sites) => 1)
+    elseif n_inds_first == 4 && n_inds_last == 4
+        # MPO-like: 4 indices (link, site, site', link)
+        elt = eltype(x[1])
+        x[1] *= onehot(elt, uniqueind(x[1], x[2], sites, prime.(sites)) => 1)
+        x[end] *= onehot(elt, uniqueind(x[end], x[end - 1], sites, prime.(sites)) => 1)
+    else
+        error("Unsupported tensor structure: first tensor has $n_inds_first indices, last tensor has $n_inds_last indices")
+    end
     return nothing
 end
 
@@ -159,23 +162,26 @@ function _removeedges!(tensors::Vector{ITensor}, sites)
         uniqueind(tensors[end], tensors[end - 1], sites, prime.(sites)) => 1)
 end
 
-function _addedges!(x::MPS)
-    length(inds(x[1])) == 2 || error("Dim of the first tensor must be 2")
-    length(inds(x[end])) == 2 || error("Dim of the last tensor must be 2")
-    linkl = Index(1, "Link,l=0")
-    linkr = Index(1, "Link,l=$(length(x))")
-    x[1] = ITensor(ITensors.data(x[1]), [linkl, inds(x[1])...])
-    x[end] = ITensor(ITensors.data(x[end]), [inds(x[end])..., linkr])
-    return nothing
-end
-
-function _addedges!(x::MPO)
-    length(inds(x[1])) == 3 || error("Dim of the first tensor must be 3")
-    length(inds(x[end])) == 3 || error("Dim of the last tensor must be 3")
-    linkl = Index(1, "Link,l=0")
-    linkr = Index(1, "Link,l=$(length(x))")
-    x[1] = ITensor(ITensors.data(x[1]), [linkl, inds(x[1])...])
-    x[end] = ITensor(ITensors.data(x[end]), [inds(x[end])..., linkr])
+function _addedges!(x::TensorTrain)
+    # Check if it's MPS-like (2 indices) or MPO-like (3 indices)
+    n_inds_first = length(inds(x[1]))
+    n_inds_last = length(inds(x[end]))
+    
+    if n_inds_first == 2 && n_inds_last == 2
+        # MPS-like: 2 indices (site, site)
+        linkl = Index(1, "Link,l=0")
+        linkr = Index(1, "Link,l=$(length(x))")
+        x[1] = ITensor(ITensors.data(x[1]), [linkl, inds(x[1])...])
+        x[end] = ITensor(ITensors.data(x[end]), [inds(x[end])..., linkr])
+    elseif n_inds_first == 3 && n_inds_last == 3
+        # MPO-like: 3 indices (site, site', site)
+        linkl = Index(1, "Link,l=0")
+        linkr = Index(1, "Link,l=$(length(x))")
+        x[1] = ITensor(ITensors.data(x[1]), [linkl, inds(x[1])...])
+        x[end] = ITensor(ITensors.data(x[end]), [inds(x[end])..., linkr])
+    else
+        error("Unsupported tensor structure: first tensor has $n_inds_first indices, last tensor has $n_inds_last indices")
+    end
     return nothing
 end
 
@@ -245,7 +251,7 @@ isascendingordescending(x) = isascendingorder(x) || isdecendingorder(x)
 function kronecker_deltas(sitesin; sitesout=prime.(noprime.(sitesin)))
     N = length(sitesout)
     links = [Index(1, "Link,l=$l") for l in 0:N]
-    M = MPO([delta(links[n], links[n + 1], sitesout[n], sitesin[n]) for n in 1:N])
+    M = TensorTrain([delta(links[n], links[n + 1], sitesout[n], sitesin[n]) for n in 1:N])
     M[1] *= onehot(links[1] => 1)
     M[end] *= onehot(links[end] => 1)
     return M
@@ -260,17 +266,21 @@ The resultant MPS do not depends on the missing site indices.
 MPO:
 For missing site indices, identity operators are inserted.
 """
-function matchsiteinds(M::Union{MPS,MPO}, sites)
+function matchsiteinds(M::TensorTrain, sites)
     N = length(sites)
     sites = noprime.(sites)
     positions = Int[findfirst(sites, s) for s in siteinds(M)]
     if length(M) > 1 && issorted(positions; lt=Base.isgreater)
-        return matchsiteinds(MPO([M[n] for n in reverse(1:length(M))]), sites)
+        return matchsiteinds(TensorTrain([M[n] for n in reverse(1:length(M))]), sites)
     end
 
-    Quantics.isascendingorder(positions) ||
+    T4AQuantics.isascendingorder(positions) ||
         error("siteinds are not in ascending order!")
 
+    # Detect if M is MPO-like (2 indices per site) or MPS-like (1 index per site)
+    sites_per_tensor = siteinds(M)
+    is_mpo = length(M) > 0 && all(length(s) == 2 for s in sites_per_tensor)
+
     # Add edges
     M_ = deepcopy(M)
     linkl = Index(1, "Link,l=0")
@@ -300,10 +310,10 @@ function matchsiteinds(M::Union{MPS,MPO}, sites)
     end
 
     links = [Index(linkdims_new[l], "Link,l=$(l-1)") for l in eachindex(linkdims_new)]
-    if M isa MPO
+    if is_mpo
         tensors = [delta(links[n], links[n + 1]) * delta(sites[n], sites[n]')
                    for n in eachindex(sites)]
-    elseif M isa MPS
+    else
         tensors = [delta(links[n], links[n + 1]) * ITensor(1, sites[n])
                    for n in eachindex(sites)]
     end
@@ -314,9 +324,9 @@ function matchsiteinds(M::Union{MPS,MPO}, sites)
         tensor = copy(M_[n])
         replaceind!(tensor, links_old[n], links[p])
         replaceind!(tensor, links_old[n + 1], links[p + 1])
-        if M isa MPO
+        if is_mpo
             tensors[p] = permute(tensor, [links[p], links[p + 1], sites[p], sites[p]'])
-        elseif M isa MPS
+        else
             tensors[p] = permute(tensor, [links[p], links[p + 1], sites[p]])
         end
     end
@@ -324,44 +334,25 @@ function matchsiteinds(M::Union{MPS,MPO}, sites)
     tensors[1] *= onehot(links[1] => 1)
     tensors[end] *= onehot(links[end] => 1)
 
-    return typeof(M)(tensors)
-end
-
-asMPO(M::MPO) = M
-
-function asMPO(tensors::Vector{ITensor})
-    N = length(tensors)
-    M = MPO(N)
-    for n in 1:N
-        M[n] = tensors[n]
-    end
-    return M
-end
-
-function asMPO(M::MPS)
-    return asMPO(ITensors.data(M))
-end
-
-function asMPS(M::MPO)
-    return MPS([t for t in M])
+    return TensorTrain(tensors)
 end
 
 """
 Contract two adjacent tensors in MPO
 """
-function combinesites(M::MPO, site1::Index, site2::Index)
+function combinesites(M::TensorTrain, site1::Index, site2::Index)
     p1 = findsite(M, site1)
     p2 = findsite(M, site2)
     p1 === nothing && error("Not found $site1")
     p2 === nothing && error("Not found $site2")
     abs(p1 - p2) == 1 ||
         error("$site1 and $site2 are found at indices $p1 and $p2. They must be on two adjacent sites.")
-    tensors = ITensors.data(M)
+    tensors = collect(M)
     idx = min(p1, p2)
     tensor = tensors[idx] * tensors[idx + 1]
     deleteat!(tensors, idx:(idx + 1))
     insert!(tensors, idx, tensor)
-    return MPO(tensors)
+    return TensorTrain(Vector{ITensor}(tensors))
 end
 
 function combinesites(
@@ -406,11 +397,13 @@ function directprod(::Type{T}, sites, indices) where {T}
     end
     tensors[1] *= onehot(links[1] => 1)
     tensors[end] *= onehot(links[end] => 1)
-    return MPS(tensors)
+    return TensorTrain(tensors)
 end
 
-function _find_target_sites(M::MPS; sitessrc=nothing, tag="")
-    _find_target_sites(siteinds(M); sitessrc, tag)
+function _find_target_sites(M::TensorTrain; sitessrc=nothing, tag="")
+    sites_flat = collect(Iterators.flatten(siteinds(M)))
+    sites_vec = Vector{Index{Int}}(sites_flat)
+    _find_target_sites(sites_vec; sitessrc, tag)
 end
 
 function _find_target_sites(
@@ -433,7 +426,7 @@ function _find_target_sites(
     return sitepos, target_sites
 end
 
-function replace_siteinds_part!(M::MPS, sitesold, sitesnew)
+function replace_siteinds_part!(M::TensorTrain, sitesold, sitesnew)
     length(sitesold) == length(sitesnew) ||
         error("Length mismatch between sitesold and sitesnew")
 
@@ -452,7 +445,7 @@ end
 Connect two MPS's
 ITensor objects are deepcopied.
 """
-function _directprod(M1::MPS, Mx::MPS...)::MPS
+function _directprod(M1::TensorTrain, Mx::TensorTrain...)::TensorTrain
     M2 = Mx[1]
     l = Index(1, "Link")
     tensors1 = [deepcopy(x) for x in M1]
@@ -460,7 +453,7 @@ function _directprod(M1::MPS, Mx::MPS...)::MPS
     tensors1[end] = ITensor(ITensors.data(last(tensors1)), [inds(last(tensors1))..., l])
     tensors2[1] = ITensor(ITensors.data(first(tensors2)), [l, inds(first(tensors2))...])
 
-    M12 = MPS([tensors1..., tensors2...])
+    M12 = TensorTrain([tensors1..., tensors2...])
     if length(Mx) == 1
         return M12
     else
@@ -468,8 +461,8 @@ function _directprod(M1::MPS, Mx::MPS...)::MPS
     end
 end
 
-function rearrange_siteinds(M::AbstractMPS, sites::Vector{Vector{Index{T}}})::MPS where {T}
-    sitesold = siteinds(MPO(collect(M)))
+function rearrange_siteinds(M::TensorTrain, sites::Vector{Vector{Index{T}}})::TensorTrain where {T}
+    sitesold = siteinds(M)
 
     Set(Iterators.flatten(sites)) == Set(Iterators.flatten(sitesold)) ||
         error("siteinds do not match $(sites) != $(sitesold)")
@@ -500,34 +493,36 @@ function rearrange_siteinds(M::AbstractMPS, sites::Vector{Vector{Index{T}}})::MP
         tensors[i], t, _ = qr(t, linds)
     end
     tensors[end] *= t
-    cleanup_linkinds!(MPS(tensors))
+    return cleanup_linkinds!(TensorTrain(tensors))
 end
 
 """
 Makes an MPS/MPO diagonal for a specified a site index `s`.
 On return, the data will be deep copied and the target core tensor will be diagonalized with an additional site index `s'`.
 """
-function makesitediagonal(M::AbstractMPS, site::Index{T})::MPS where {T}
-    M_ = deepcopy(MPO(collect(M)))
+function makesitediagonal(M::TensorTrain, site::Index{T})::TensorTrain where {T}
+    M_ = deepcopy(M)
 
     target_site::Int = only(findsites(M_, site))
     M_[target_site] = _asdiagonal(M_[target_site], site)
 
-    return MPS(collect(M_))
+    return M_
 end
 
-function makesitediagonal(M::AbstractMPS, tag::String)::MPS
-    M_ = deepcopy(MPO(collect(M)))
+function makesitediagonal(M::TensorTrain, tag::String)::TensorTrain
+    M_ = deepcopy(M)
     sites = siteinds(M_)
 
-    target_positions = findallsites_by_tag(siteinds(M_); tag=tag)
+    # Convert to Vector{Vector{Index{Int}}} for type compatibility
+    sites_converted = Vector{Vector{Index{Int}}}([Vector{Index{Int}}(collect(s)) for s in sites])
+    target_positions = findallsites_by_tag(sites_converted; tag=tag)
 
     for t in eachindex(target_positions)
         i, j = target_positions[t]
         M_[i] = _asdiagonal(M_[i], sites[i][j])
     end
 
-    return MPS(collect(M_))
+    return M_
 end
 
 # FIXME: may be type unstable
@@ -541,18 +536,20 @@ end
 """
 Extract diagonal components
 """
-function extractdiagonal(M::AbstractMPS, tag::String)::MPS
-    M_ = deepcopy(MPO(collect(M)))
+function extractdiagonal(M::TensorTrain, tag::String)::TensorTrain
+    M_ = deepcopy(M)
     sites = siteinds(M_)
 
-    target_positions = findallsites_by_tag(siteinds(M_); tag=tag)
+    # Convert to Vector{Vector{Index{Int}}} for type compatibility
+    sites_converted = Vector{Vector{Index{Int}}}([Vector{Index{Int}}(collect(s)) for s in sites])
+    target_positions = findallsites_by_tag(sites_converted; tag=tag)
 
     for t in eachindex(target_positions)
         i, j = target_positions[t]
         M_[i] = _extract_diagonal(M_[i], sites[i][j], sites[i][j]')
     end
 
-    return MPS(collect(M_))
+    return M_
 end
 
 function _extract_diagonal(t, site::Index{T}, site2::Index{T}) where {T<:Number}
@@ -566,25 +563,15 @@ function _extract_diagonal(t, site::Index{T}, site2::Index{T}) where {T<:Number}
     return ITensor(newdata, restinds..., site)
 end
 
-function _apply(A::MPO, Ψ::MPO; alg::String="fit", cutoff::Real=1e-25, kwargs...)::MPO
-    if :algorithm ∈ keys(kwargs)
-        error("keyword argument :algorithm is not allowed")
-    end
-    if alg == "densitymatrix" && cutoff <= 1e-10
-        @warn "cutoff is too small for densitymatrix algorithm. Use fit algorithm instead."
-    end
-    AΨ = replaceprime(
-        FastMPOContractions.contract_mpo_mpo(A', asMPO(Ψ); alg, cutoff, kwargs...), 2 => 1)
-    MPO(collect(AΨ))
-end
-
-function _apply(A::MPO, Ψ::MPS; alg::String="fit", cutoff::Real=1e-25, kwargs...)::MPS
+function _apply(A::TensorTrain, Ψ::TensorTrain; alg::String="fit", cutoff::Real=1e-25, kwargs...)::TensorTrain
     if :algorithm ∈ keys(kwargs)
         error("keyword argument :algorithm is not allowed")
     end
     if alg == "densitymatrix" && cutoff <= 1e-10
         @warn "cutoff is too small for densitymatrix algorithm. Use fit algorithm instead."
     end
-    AΨ = noprime.(FastMPOContractions.contract_mpo_mpo(A, asMPO(Ψ); alg, cutoff, kwargs...))
-    MPS(collect(AΨ))
+    # Convert alg to Algorithm type
+    alg_ = alg isa Algorithm ? alg : Algorithm(alg)
+    # Use T4AITensorCompat.product for MPO * MPO
+    return product(A, Ψ; alg=alg_, cutoff=cutoff, kwargs...)
 end
diff --git a/test/_util.jl b/test/_util.jl
index ae09349..ea86ca7 100644
--- a/test/_util.jl
+++ b/test/_util.jl
@@ -1,5 +1,5 @@
 using ITensors
-import ITensorMPS: random_mps
+import T4AITensorCompat: random_mps, TensorTrain
 using Random
 
 function _random_mpo(
@@ -13,7 +13,7 @@ function _random_mpo(
         push!(tensors, prod(Ψ[pos:(pos + length(sites[i]) - 1)]))
         pos += length(sites[i])
     end
-    return MPO(tensors)
+    return TensorTrain(tensors)
 end
 
 function _random_mpo(
diff --git a/test/affine_tests.jl b/test/affine_tests.jl
index 7e49df4..a6201f8 100644
--- a/test/affine_tests.jl
+++ b/test/affine_tests.jl
@@ -1,7 +1,7 @@
 @testitem "affine" begin
     using Test
     using ITensors
-    using Quantics
+    using T4AQuantics
     using LinearAlgebra
 
     # Test results of affine_transform_matrix()
@@ -15,9 +15,9 @@
         for x_vals in Iterators.product(xranges...), y_vals in Iterators.product(yranges...)
             x = collect(x_vals)  # Convert input tuple to a vector
             Axb = A * x .+ b  # Apply the affine transformation
-            if boundary == Quantics.PeriodicBoundaryConditions()
+            if boundary == T4AQuantics.PeriodicBoundaryConditions()
                 Axb = mod.(Axb, 2^R)
-            elseif boundary == Quantics.OpenBoundaryConditions()
+            elseif boundary == T4AQuantics.OpenBoundaryConditions()
                 Axb = map(x -> 0 <= x <= 2^R - 1 ? x : nothing, Axb)
             end
             ref = Axb == collect(y_vals)  # Compare to the reference output
@@ -40,18 +40,18 @@
     vars = ("x", "y", "z")
     insite = [Index(2; tags="$v$l") for l in 1:10, v in vars]
     outsite = [Index(2; tags="$v$l")' for l in 1:10, v in vars]
-    boundaries = Quantics.OpenBoundaryConditions(), Quantics.PeriodicBoundaryConditions()
+    boundaries = T4AQuantics.OpenBoundaryConditions(), T4AQuantics.PeriodicBoundaryConditions()
 
     @testset "full" begin
         A = [1 0; 1 1]
         b = [0; 0]
 
-        T = Quantics.affine_transform_matrix(4, A, b)
+        T = T4AQuantics.affine_transform_matrix(4, A, b)
         @test T' * T == I
         @test T[219, 59] == 1
 
         b = [4; 1]
-        T = Quantics.affine_transform_matrix(4, A, b)
+        T = T4AQuantics.affine_transform_matrix(4, A, b)
         @test T * T' == I
     end
 
@@ -61,7 +61,7 @@
         for (A, b) in testtests[(M, N)]
             A_ = reshape(A, M, N)
             b_ = reshape(b, M)
-            T = Quantics.affine_transform_matrix(R, A_, b_, boundary)
+            T = T4AQuantics.affine_transform_matrix(R, A_, b_, boundary)
             test_affine_transform_matrix_multi_variables(R, A_, b_, T, boundary)
         end
     end
@@ -70,10 +70,10 @@
         b = [0; 0]
         R = 3
 
-        T = Quantics.affine_transform_matrix(R, A, b)
-        mpo = Quantics.affine_transform_mpo(
+        T = T4AQuantics.affine_transform_matrix(R, A, b)
+        mpo = T4AQuantics.affine_transform_mpo(
             outsite[1:R, 1:2], insite[1:R, 1:2], A, b)
-        Trec = Quantics.affine_mpo_to_matrix(
+        Trec = T4AQuantics.affine_mpo_to_matrix(
             outsite[1:R, 1:2], insite[1:R, 1:2], mpo)
         @test T == Trec
     end
@@ -83,11 +83,11 @@
         b = [11; 23; -15]
         R = 4
 
-        T = Quantics.affine_transform_matrix(R, A, b)
+        T = T4AQuantics.affine_transform_matrix(R, A, b)
         M, N = size(A)
-        mpo = Quantics.affine_transform_mpo(
+        mpo = T4AQuantics.affine_transform_mpo(
             outsite[1:R, 1:M], insite[1:R, 1:N], A, b)
-        Trec = Quantics.affine_mpo_to_matrix(
+        Trec = T4AQuantics.affine_mpo_to_matrix(
             outsite[1:R, 1:M], insite[1:R, 1:N], mpo)
         @test T == Trec
     end
@@ -97,11 +97,11 @@
         b = [11; -3]
         R = 4
 
-        T = Quantics.affine_transform_matrix(R, A, b)
+        T = T4AQuantics.affine_transform_matrix(R, A, b)
         M, N = size(A)
-        mpo = Quantics.affine_transform_mpo(
+        mpo = T4AQuantics.affine_transform_mpo(
             outsite[1:R, 1:M], insite[1:R, 1:N], A, b)
-        Trec = Quantics.affine_mpo_to_matrix(
+        Trec = T4AQuantics.affine_mpo_to_matrix(
             outsite[1:R, 1:M], insite[1:R, 1:N], mpo)
         @test T == Trec
     end
@@ -112,11 +112,11 @@
 
         for R in [1, 3, 6]
             for bc in boundaries
-                T = Quantics.affine_transform_matrix(R, A, b, bc)
+                T = T4AQuantics.affine_transform_matrix(R, A, b, bc)
                 M, N = size(A)
-                mpo = Quantics.affine_transform_mpo(
+                mpo = T4AQuantics.affine_transform_mpo(
                     outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc)
-                Trec = Quantics.affine_mpo_to_matrix(
+                Trec = T4AQuantics.affine_mpo_to_matrix(
                     outsite[1:R, 1:M], insite[1:R, 1:N], mpo)
                 @test T == Trec
             end
@@ -129,12 +129,12 @@
         for b in [[3], [5], [-3], [-5]]
             for R in [3, 5]
                 #for bc in boundaries
-                for bc in [Quantics.PeriodicBoundaryConditions()]
-                    T = Quantics.affine_transform_matrix(R, A, b, bc)
+                for bc in [T4AQuantics.PeriodicBoundaryConditions()]
+                    T = T4AQuantics.affine_transform_matrix(R, A, b, bc)
                     M, N = size(A)
-                    mpo = Quantics.affine_transform_mpo(
+                    mpo = T4AQuantics.affine_transform_mpo(
                         outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc)
-                    Trec = Quantics.affine_mpo_to_matrix(
+                    Trec = T4AQuantics.affine_mpo_to_matrix(
                         outsite[1:R, 1:M], insite[1:R, 1:N], mpo)
                     @test T == Trec
                 end
@@ -148,11 +148,11 @@
 
         for R in [3, 4]
             for bc in boundaries
-                T = Quantics.affine_transform_matrix(R, A, b, bc)
+                T = T4AQuantics.affine_transform_matrix(R, A, b, bc)
                 M, N = size(A)
-                mpo = Quantics.affine_transform_mpo(
+                mpo = T4AQuantics.affine_transform_mpo(
                     outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc)
-                Trec = Quantics.affine_mpo_to_matrix(
+                Trec = T4AQuantics.affine_mpo_to_matrix(
                     outsite[1:R, 1:M], insite[1:R, 1:N], mpo)
                 @test T == Trec
             end
diff --git a/test/aqua_tests.jl b/test/aqua_tests.jl
index 370317d..0b98558 100644
--- a/test/aqua_tests.jl
+++ b/test/aqua_tests.jl
@@ -1,9 +1,9 @@
 @testitem "code quality test" begin
     using Test
     using Aqua
-    import Quantics
+    import T4AQuantics
 
     @testset "Aqua" begin
-        Aqua.test_stale_deps(Quantics)
+        Aqua.test_stale_deps(T4AQuantics)
     end
 end
diff --git a/test/binaryop_tests.jl b/test/binaryop_tests.jl
index 9be359d..47c2e31 100644
--- a/test/binaryop_tests.jl
+++ b/test/binaryop_tests.jl
@@ -1,9 +1,9 @@
 @testitem "binaryop_tests.jl/_binaryop" begin
     using Test
     using ITensors
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps, apply
     ITensors.disable_warn_order()
-    using Quantics
+    using T4AQuantics
     import Random
 
     @testset "_binaryop" for rev_carrydirec in [true], nbit in 2:3
@@ -28,10 +28,10 @@
         rsites = reverse(sites)
 
         for a in -1:1, b in -1:1, c in -1:1, d in -1:1, bc_x in [1, -1], bc_y in [1, -1]
-            g = randomMPS(sites)
-            M = Quantics._binaryop_mpo(sites, [(a, b), (c, d)], [(1, 2), (1, 2)];
+            g = random_mps(sites)
+            M = T4AQuantics._binaryop_mpo(sites, [(a, b), (c, d)], [(1, 2), (1, 2)];
                 rev_carrydirec=rev_carrydirec, bc=[bc_x, bc_y])
-            f = apply(M, g)
+            f = apply(M, g; alg="naive")
 
             # f[x_R, ..., x_1, y_R, ..., y_1] and f[x, y]
             f_arr = Array(reduce(*, f), vcat(reverse(sitesx), reverse(sitesy)))
@@ -65,9 +65,9 @@ end
 @testitem "binaryop_tests.jl/affinetransform" begin
     using Test
     using ITensors
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps
     ITensors.disable_warn_order()
-    using Quantics
+    using T4AQuantics
     import Random
 
     @testset "affinetransform" for rev_carrydirec in [true, false], nbit in 2:3
@@ -91,10 +91,10 @@ end
         shift = rand((-2 * 2^nbit):(2 * 2^nbit), 2)
 
         for a in -1:1, b in -1:1, c in -1:1, d in -1:1, bc_x in [1, -1], bc_y in [1, -1]
-            g = randomMPS(sites)
-            f = Quantics.affinetransform(g, ["x", "y"],
+            g = random_mps(sites)
+            f = T4AQuantics.affinetransform(g, ["x", "y"],
                 [Dict("x" => a, "y" => b), Dict("x" => c, "y" => d)],
-                shift, [bc_x, bc_y]; cutoff=1e-25)
+                shift, [bc_x, bc_y]; cutoff=1e-25, alg="naive")
 
             # f[x_R, ..., x_1, y_R, ..., y_1] and f[x, y]
             f_arr = Array(reduce(*, f), vcat(reverse(sitesx), reverse(sitesy)))
@@ -128,9 +128,9 @@ end
 @testitem "binaryop_tests.jl/affinetransform_three_vars" begin
     using Test
     using ITensors
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps
     ITensors.disable_warn_order()
-    using Quantics
+    using T4AQuantics
     import Random
 
     affinetransform_testsets = []
@@ -212,10 +212,10 @@ end
 
         shift = rand((-2 * 2^nbit):(2 * 2^nbit), 3)
 
-        g = randomMPS(sites)
-        f = Quantics.affinetransform(g, ["x", "y", "z"],
+        g = random_mps(sites)
+        f = T4AQuantics.affinetransform(g, ["x", "y", "z"],
             coeffs_dic,
-            shift, [bc_x, bc_y, bc_z]; cutoff=1e-25)
+            shift, [bc_x, bc_y, bc_z]; cutoff=1e-25, alg="naive")
 
         # f[x_R, ..., x_1, y_R, ..., y_1, z_R, ..., z_1] and f[x, y, z]
         f_arr = Array(reduce(*, f),
@@ -251,18 +251,18 @@ end
 @testitem "binaryop_tests.jl/shiftop" begin
     using Test
     using ITensors
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps, apply, product
     ITensors.disable_warn_order()
-    using Quantics
+    using T4AQuantics
     import Random
 
     @testset "shiftop" for R in [3], bc in [1, -1]
         sites = [Index(2, "Qubit, x=$n") for n in 1:R]
-        g = randomMPS(sites)
+        g = random_mps(sites)
 
         for shift in [0, 1, 2, 2^R - 1]
-            M = Quantics._shift_mpo(sites, shift; bc=bc)
-            f = apply(M, g)
+            M = T4AQuantics._shift_mpo(sites, shift; bc=bc)
+            f = apply(M, g; alg="naive")
 
             f_vec = vec(Array(reduce(*, f), reverse(sites)))
             g_vec = vec(Array(reduce(*, g), reverse(sites)))
diff --git a/test/fouriertransform_tests.jl b/test/fouriertransform_tests.jl
index daab5a4..5ceb4b8 100644
--- a/test/fouriertransform_tests.jl
+++ b/test/fouriertransform_tests.jl
@@ -1,10 +1,10 @@
 
 @testitem "fouriertransform_tests.jl/qft_mpo" begin
     using Test
-    using Quantics
+    using T4AQuantics
     using ITensors
-    using ITensorMPS: randomMPS, MPO
-    using ITensors.SiteTypes: siteinds
+    import T4AITensorCompat: TensorTrain
+    import ITensorMPS
     # A brute-force implementation of _qft (only for tests)
     function _qft_ref(sites; cutoff::Float64=1e-14, sign::Int=1)
         abs(sign) == 1 || error("sign must either 1 or -1")
@@ -23,15 +23,15 @@
         tmat = reshape(tmat, ntuple(x -> 2, 2 * nbit))
 
         trans_t = ITensor(tmat, reverse(sites)..., prime(sites)...)
-        M = MPO(trans_t, sites; cutoff=cutoff)
+        M = TensorTrain(ITensorMPS.MPO(trans_t, sites; cutoff=cutoff))
         return M
     end
 
     @testset "qft_mpo" for sign in [1, -1], nbit in [2, 3]
         N = 2^nbit
 
-        sites = siteinds("Qubit", nbit)
-        M = Quantics._qft(sites; sign=sign)
+        sites = [Index(2, "Qubit,n=$n") for n in 1:nbit]
+        M = T4AQuantics._qft(sites; sign=sign)
         M_ref = _qft_ref(sites; sign=sign)
 
         @test M ≈ M_ref
@@ -40,9 +40,9 @@ end
 
 @testitem "fouriertransform_tests.jl/fouriertransform" begin
     using Test
-    using Quantics
+    using T4AQuantics
     using ITensors
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps
 
     function _ft_1d_ref(X, sign, originx, origink)
         N = length(X)
@@ -66,11 +66,11 @@ end
         sitesk = [Index(2, "Qubit,k=$k") for k in 1:nbit]
 
         # X(x)
-        X = randomMPS(sitesx)
+        X = random_mps(sitesx)
         X_vec = Array(reduce(*, X), reverse(sitesx))
 
         # Y(k)
-        Y = Quantics.fouriertransform(X; sign=sign, tag="x", sitesdst=sitesk,
+        Y = T4AQuantics.fouriertransform(X; sign=sign, tag="x", sitesdst=sitesk,
             originsrc=originx, origindst=originy)
 
         Y_vec_ref = _ft_1d_ref(X_vec, sign, originx, originy)
@@ -104,13 +104,13 @@ end
 
         # F(x, y)
         # F(x_1, y_1, ..., x_R, y_R)
-        F = randomMPS(sitesin)
+        F = random_mps(sitesin)
         F_mat = reshape(Array(reduce(*, F), vcat(reverse(sitesx), reverse(sitesy))), N, N)
 
         # G(kx, ky)
         # G(kx_R, ky_R, ..., kx_1, ky_1)
-        G_ = Quantics.fouriertransform(F; sign=sign, tag="x", sitesdst=siteskx)
-        G = Quantics.fouriertransform(G_; sign=sign, tag="y", sitesdst=sitesky)
+        G_ = T4AQuantics.fouriertransform(F; sign=sign, tag="x", sitesdst=siteskx)
+        G = T4AQuantics.fouriertransform(G_; sign=sign, tag="y", sitesdst=sitesky)
 
         G_mat_ref = _ft_2d_ref(F_mat, sign)
         G_mat = reshape(Array(reduce(*, G), vcat(reverse(siteskx), reverse(sitesky))), N, N)
diff --git a/test/imaginarytime_tests.jl b/test/imaginarytime_tests.jl
index c1608e1..ad16528 100644
--- a/test/imaginarytime_tests.jl
+++ b/test/imaginarytime_tests.jl
@@ -1,16 +1,21 @@
 @testitem "imaginarytime_tests.jl/imaginarytime" begin
     using Test
-    using Quantics
-    import Quantics
-    using ITensors: siteinds, Index
-    using ITensors.SiteTypes: op
+    using T4AQuantics
+    import T4AQuantics
+    using ITensors: Index
     import ITensors
     import SparseIR: Fermionic, Bosonic, FermionicFreq, valueim
 
-    import ITensorMPS: MPS, onehot
-    import QuanticsGrids as QG
+    import T4AITensorCompat: TensorTrain
+    # Conditionally import QuanticsGrids and QuanticsTCI
+    try
+        import QuanticsGrids as QG
+        import QuanticsTCI: quanticscrossinterpolate
+        has_qg = true
+    catch
+        has_qg = false
+    end
     import TensorCrossInterpolation as TCI
-    import QuanticsTCI: quanticscrossinterpolate
 
     function _test_data_imaginarytime(nbit, β)
         ω = 0.5
@@ -38,8 +43,8 @@
 
         gtau_smpl, giv_smpl = _test_data_imaginarytime(nbit, β)
 
-        sites = siteinds("Qubit", nbit)
-        gtau_mps = Quantics.decompose_gtau(gtau_smpl, sites; cutoff=1e-20)
+        sites = [Index(2, "Qubit,n=$n") for n in 1:nbit]
+        gtau_mps = T4AQuantics.decompose_gtau(gtau_smpl, sites; cutoff=1e-20)
 
         gtau_smpl_reconst = vec(Array(reduce(*, gtau_mps), reverse(sites)...))
 
@@ -56,8 +61,8 @@
 
         sitesτ = [Index(2, "Qubit,τ=$n") for n in 1:nbit]
         sitesiω = [Index(2, "Qubit,iω=$n") for n in 1:nbit]
-        gtau_mps = Quantics.decompose_gtau(gtau_smpl, sitesτ; cutoff=1e-20)
-        giv_mps = Quantics.to_wn(Fermionic(), gtau_mps, β; cutoff=1e-20, tag="τ",
+        gtau_mps = T4AQuantics.decompose_gtau(gtau_smpl, sitesτ; cutoff=1e-20)
+        giv_mps = T4AQuantics.to_wn(Fermionic(), gtau_mps, β; cutoff=1e-20, tag="τ",
             sitesdst=sitesiω)
 
         giv = vec(Array(reduce(*, giv_mps), reverse(sitesiω)...))
@@ -75,9 +80,9 @@
 
         sitesτ = [Index(2, "Qubit,τ=$n") for n in 1:nbit]
         sitesiω = [Index(2, "Qubit,iω=$n") for n in 1:nbit]
-        giv_mps = Quantics.decompose_giv(giv_smpl, sitesiω; cutoff=1e-20)
+        giv_mps = T4AQuantics.decompose_giv(giv_smpl, sitesiω; cutoff=1e-20)
 
-        gtau_mps = Quantics.to_tau(Fermionic(), giv_mps, β; cutoff=1e-20, tag="iω",
+        gtau_mps = T4AQuantics.to_tau(Fermionic(), giv_mps, β; cutoff=1e-20, tag="iω",
             sitesdst=sitesτ)
 
         gtau = vec(Array(reduce(*, gtau_mps), reverse(sitesτ)...))
@@ -89,6 +94,15 @@
 
 
     @testset "ImaginaryTimeFT.to_tau with large R" begin
+        # Skip this test if QuanticsGrids or QuanticsTCI is not available
+        try
+            import QuanticsGrids as QG
+            import QuanticsTCI: quanticscrossinterpolate
+        catch
+            @test_skip "QuanticsGrids or QuanticsTCI not available"
+            return
+        end
+        
         function fermionic_wn(n, β)
             return (2 * n + 1) * π / β
         end
@@ -97,8 +111,8 @@
             return 1.0 / (im * fermionic_wn(n, β) - ϵ)
         end
 
-        _evaluate(Ψ::MPS, sites, index::Vector{Int}) = only(reduce(
-            *, Ψ[n] * onehot(sites[n] => index[n]) for n in 1:length(Ψ)))
+        _evaluate(Ψ::TensorTrain, sites, index::Vector{Int}) = only(reduce(
+            *, Ψ[n] * ITensors.onehot(sites[n] => index[n]) for n in 1:length(Ψ)))
 
         β = 100.0
         R = 50
@@ -122,9 +136,12 @@
             ComplexF64, inv_iwn_tci, ngrid; tolerance=tol, maxbonddim=maxdim_TCI)
 
         inv_iwn_tt = TCI.TensorTrain(qtci2.tci)
-        iwmps = MPS(inv_iwn_tt; sites=sitesiω)
+        # Convert TensorCrossInterpolation.TensorTrain to ITensorMPS.MPS first
+        import ITensorMPS
+        iwmps_mps = ITensorMPS.MPS(inv_iwn_tt; sites=sitesiω)
+        iwmps = TensorTrain(iwmps_mps)
 
-        fourier_inv_iw = Quantics.to_tau(
+        fourier_inv_iw = T4AQuantics.to_tau(
             Fermionic(), iwmps, β; tag="iω", sitesdst=sitesτ, cutoff_MPO=cutoff_mpo,
             cutoff=cutoff_contract, maxdim=maxdim_contract, alg="naive")
 
@@ -136,17 +153,17 @@ end
 
 @testitem "imaginarytime_tests.jl/poletomps" begin
     using Test
-    using Quantics
-    import ITensors: siteinds, Index
+    using T4AQuantics
+    import ITensors: Index
     import ITensors
     import SparseIR: Fermionic, Bosonic, FermionicFreq, valueim
 
     @testset "poletomps" begin
         nqubit = 10
-        sites = siteinds("Qubit", nqubit)
+        sites = [Index(2, "Qubit,n=$n") for n in 1:nqubit]
         β = 10.0
         ω = 1.2
-        gtau = Quantics.poletomps(sites, β, ω)
+        gtau = T4AQuantics.poletomps(sites, β, ω)
         gtauvec = vec(Array(reduce(*, gtau), reverse(sites)))
         gtauf(τ) = -exp(-τ * ω) / (1 + exp(-β * ω))
         gtauref = gtauf.(LinRange(0, β, 2^nqubit + 1)[1:(end - 1)])
@@ -155,10 +172,10 @@ end
 
     @testset "poletomps_negative_pole" begin
         nqubit = 16
-        sites = siteinds("Qubit", nqubit)
+        sites = [Index(2, "Qubit,n=$n") for n in 1:nqubit]
         β = 1000.0
         ω = -10.0
-        gtau = Quantics.poletomps(Fermionic(), sites, β, ω)
+        gtau = T4AQuantics.poletomps(Fermionic(), sites, β, ω)
         gtauvec = vec(Array(reduce(*, gtau), reverse(sites)))
         gtauf(τ) = -exp((β - τ) * ω)
         gtauref = gtauf.(LinRange(0, β, 2^nqubit + 1)[1:(end - 1)])
diff --git a/test/mps_tests.jl b/test/mps_tests.jl
index 964acf9..ecfa044 100644
--- a/test/mps_tests.jl
+++ b/test/mps_tests.jl
@@ -1,25 +1,25 @@
 @testitem "mps_tests.jl/onemps" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
     using ITensors.SiteTypes: siteinds
     @testset "onemps" begin
         nbit = 3
         sites = siteinds("Qubit", nbit)
-        M = Quantics.onemps(Float64, sites)
+        M = T4AQuantics.onemps(Float64, sites)
         @test vec(Array(reduce(*, M), sites)) ≈ ones(2^nbit)
     end
 end
 
 @testitem "mps_tests.jl/expqtt" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
     using ITensors.SiteTypes: siteinds
     @testset "expqtt" begin
         R = 10
         sites = siteinds("Qubit", 10)
-        f = Quantics.expqtt(sites, -1.0)
+        f = T4AQuantics.expqtt(sites, -1.0)
         f_values = vec(Array(reduce(*, f), reverse(sites)))
         xs = collect(LinRange(0, 1, 2^R + 1)[1:(end - 1)])
         f_values_ref = (x -> exp(-x)).(xs)
diff --git a/test/mul_tests.jl b/test/mul_tests.jl
index ad4a83f..a96b028 100644
--- a/test/mul_tests.jl
+++ b/test/mul_tests.jl
@@ -1,8 +1,8 @@
 @testitem "test_mul.jl/preprocess_matmul" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
-    using ITensorMPS: random_mps, MPO
+    import T4AITensorCompat: random_mps, TensorTrain
 
     @testset "_preprocess_matmul" begin
         N = 2
@@ -11,12 +11,12 @@
         sitesz = [Index(2, "z=$n") for n in 1:N]
         sites1 = collect(Iterators.flatten(zip(sitesx, sitesy)))
         sites2 = collect(Iterators.flatten(zip(sitesy, sitesz)))
-        M1 = Quantics.asMPO(random_mps(sites1))
-        M2 = Quantics.asMPO(random_mps(sites2))
+        M1 = random_mps(sites1)
+        M2 = random_mps(sites2)
 
-        mul = Quantics.MatrixMultiplier(sitesx, sitesy, sitesz)
+        mul = T4AQuantics.MatrixMultiplier(sitesx, sitesy, sitesz)
 
-        M1, M2 = Quantics.preprocess(mul, M1, M2)
+        M1, M2 = T4AQuantics.preprocess(mul, M1, M2)
 
         flag = true
         for n in 1:N
@@ -31,15 +31,17 @@
         sitesx = [Index(2, "x=$n") for n in 1:N]
         sitesy = [Index(2, "y=$n") for n in 1:N]
         sitesz = [Index(2, "z=$n") for n in 1:N]
-        mul = Quantics.MatrixMultiplier(sitesx, sitesy, sitesz)
+        mul = T4AQuantics.MatrixMultiplier(sitesx, sitesy, sitesz)
 
         links = [Index(1, "Link,l=$l") for l in 0:N]
-        M = MPO(N)
+        sites_pairs = [[sitesx[n], sitesz[n]] for n in 1:N]
+        M = TensorTrain(ComplexF64, sites_pairs, 1)
         for n in 1:N
-            M[n] = randomITensor(links[n], links[n + 1], sitesx[n], sitesz[n])
+            inds_list = [links[n], links[n + 1], sitesx[n], sitesz[n]]
+            M[n] = randomITensor(inds_list...)
         end
 
-        M = Quantics.postprocess(mul, M)
+        M = T4AQuantics.postprocess(mul, M)
 
         flag = true
         for n in 1:N
@@ -52,30 +54,30 @@ end
 
 @testitem "mul_tests.jl/matmul" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
-    using ITensorMPS: random_mps
+    import T4AITensorCompat: random_mps, contract, Algorithm
 
     @testset "matmul" for T in [Float64, ComplexF64]
         N = 3
         sitesx = [Index(2, "x=$n") for n in 1:N]
         sitesy = [Index(2, "y=$n") for n in 1:N]
         sitesz = [Index(2, "z=$n") for n in 1:N]
-        mul = Quantics.MatrixMultiplier(sitesx, sitesy, sitesz)
+        mul = T4AQuantics.MatrixMultiplier(sitesx, sitesy, sitesz)
 
         sites1 = collect(Iterators.flatten(zip(sitesx, sitesy)))
         sites2 = collect(Iterators.flatten(zip(sitesy, sitesz)))
-        M1 = Quantics.asMPO(random_mps(T, sites1))
-        M2 = Quantics.asMPO(random_mps(T, sites2))
+        M1 = random_mps(T, sites1)
+        M2 = random_mps(T, sites2)
 
         # preprocess
-        M1, M2 = Quantics.preprocess(mul, M1, M2)
+        M1, M2 = T4AQuantics.preprocess(mul, M1, M2)
 
         # MPO-MPO contraction
-        M = Quantics.asMPO(contract(M1, M2; alg="naive"))
+        M = contract(M1, M2; alg=Algorithm("naive"))
 
         # postprocess
-        M = Quantics.postprocess(mul, M)
+        M = T4AQuantics.postprocess(mul, M)
 
         M_mat_reconst = reshape(Array(reduce(*, M), [reverse(sitesx)..., reverse(sitesz)]),
             2^N, 2^N)
@@ -93,21 +95,21 @@ end
     @testset "elementwisemul" for T in [Float64, ComplexF64]
         N = 5
         sites = [Index(2, "n=$n") for n in 1:N]
-        mul = Quantics.ElementwiseMultiplier(sites)
+        mul = T4AQuantics.ElementwiseMultiplier(sites)
 
         M1_ = random_mps(T, sites)
         M2_ = random_mps(T, sites)
-        M1 = Quantics.asMPO(M1_)
-        M2 = Quantics.asMPO(M2_)
+        M1 = M1_
+        M2 = M2_
 
         # preprocess
-        M1, M2 = Quantics.preprocess(mul, M1, M2)
+        M1, M2 = T4AQuantics.preprocess(mul, M1, M2)
 
         # MPO-MPO contraction
-        M = Quantics.asMPO(contract(M1, M2; alg="naive"))
+        M = contract(M1, M2; alg=Algorithm("naive"))
 
         # postprocess
-        M = Quantics.postprocess(mul, M)
+        M = T4AQuantics.postprocess(mul, M)
 
         # Comparison with reference data
         M_reconst = Array(reduce(*, M), sites)
@@ -120,22 +122,23 @@ end
 
 @testitem "mul_tests.jl/batchedmatmul" begin
     using Test
-    import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, isprojectedat,
-                            project, Projector
-    import Quantics
+    import PartitionedMPSs: PartitionedMPSs, isprojectedat
+    import T4APartitionedMPSs: project, Projector, SubDomainMPS, PartitionedMPS
+    import T4AQuantics
     using ITensors
-    using ITensors.SiteTypes: siteinds
-    using ITensorMPS: random_mps, MPS
+    import T4AITensorCompat: random_mps, TensorTrain, siteinds
 
     """
     Reconstruct 3D matrix
     """
     function _tomat3(a)
-        sites = siteinds(a)
-        N = length(sites)
+        sites_vec = siteinds(a)
+        sites_flat = collect(Iterators.flatten(sites_vec))
+        N = length(sites_flat)
         Nreduced = N ÷ 3
-        sites_ = [sites[1:3:N]..., sites[2:3:N]..., sites[3:3:N]...]
-        return reshape(Array(reduce(*, a), sites_), 2^Nreduced, 2^Nreduced, 2^Nreduced)
+        sites_ = [sites_flat[1:3:N]..., sites_flat[2:3:N]..., sites_flat[3:3:N]...]
+        result = reduce(*, a)
+        return reshape(Array(result, sites_), 2^Nreduced, 2^Nreduced, 2^Nreduced)
     end
 
     @testset "batchedmatmul" for T in [Float64, ComplexF64]
@@ -163,7 +166,7 @@ end
             ab_arr[:, :, k] .= a_arr[:, :, k] * b_arr[:, :, k]
         end
 
-        ab = Quantics.automul(a, b; tag_row="x", tag_shared="y", tag_col="z", alg="fit")
+        ab = T4AQuantics.automul(a, b; tag_row="x", tag_shared="y", tag_col="z", alg="fit")
         ab_arr_reconst = _tomat3(ab)
         @test ab_arr ≈ ab_arr_reconst
     end
@@ -210,17 +213,17 @@ end
                 project(b, Projector(Dict(sy[1] => 2, sz[1] => 2)))
             ])
 
-            @test a ≈ MPS(a_)
-            @test b ≈ MPS(b_)
+            @test a ≈ TensorTrain(a_)
+            @test b ≈ TensorTrain(b_)
 
-            ab = Quantics.automul(
+            ab = T4AQuantics.automul(
                 a_, b_; tag_row="x", tag_shared="y", tag_col="z", alg="fit", cutoff, maxdim
             )
-            ab_ref = Quantics.automul(
+            ab_ref = T4AQuantics.automul(
                 a, b; tag_row="x", tag_shared="y", tag_col="z", alg="fit", cutoff, maxdim
             )
 
-            @test MPS(ab)≈ab_ref rtol=10 * sqrt(cutoff)
+            @test TensorTrain(ab)≈ab_ref rtol=10 * sqrt(cutoff)
         end
     end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 1947886..ca6007a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
 using ReTestItems: runtests, @testitem
-using Quantics: Quantics
+using T4AQuantics: T4AQuantics
 
-runtests(Quantics)
+runtests(T4AQuantics)
diff --git a/test/tag_tests.jl b/test/tag_tests.jl
index f69a4bd..c55e578 100644
--- a/test/tag_tests.jl
+++ b/test/tag_tests.jl
@@ -1,38 +1,38 @@
 @testitem "tag_tests.jl/tag" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
 
     @testset "findallsites_by_tag" for tag in ["x", "y"]
         nbit = 4
         sites = [Index(2, "Qubit,$(tag)=$x") for x in 1:nbit]
-        @test Quantics.findallsites_by_tag(sites; tag=tag) == [1, 2, 3, 4]
-        @test isempty(Quantics.findallsites_by_tag(sites; tag="notfound"))
+        @test T4AQuantics.findallsites_by_tag(sites; tag=tag) == [1, 2, 3, 4]
+        @test isempty(T4AQuantics.findallsites_by_tag(sites; tag="notfound"))
 
         invalid_tag = "$(tag)="
-        @test_throws "Invalid tag: $(tag)=" Quantics.findallsites_by_tag(sites,
+        @test_throws "Invalid tag: $(tag)=" T4AQuantics.findallsites_by_tag(sites,
             tag=invalid_tag)
 
         invalid_sites = [Index(2, "Qubit,$(tag)=1"), Index(2, "Qubit,$(tag)=1")]
-        @test_throws "with $(tag)=1!" Quantics.findallsites_by_tag(invalid_sites, tag=tag)
-        @test_throws "Invalid tag: $(tag)=" Quantics.findallsites_by_tag(invalid_sites,
+        @test_throws "with $(tag)=1!" T4AQuantics.findallsites_by_tag(invalid_sites, tag=tag)
+        @test_throws "Invalid tag: $(tag)=" T4AQuantics.findallsites_by_tag(invalid_sites,
             tag="$(tag)=")
     end
 
     @testset "findallsiteinds_by_tag" for tag in ["x", "y"]
         nbit = 4
         sites = [Index(2, "Qubit,$(tag)=$x") for x in 1:nbit]
-        @test Quantics.findallsiteinds_by_tag(sites; tag=tag) == sites
-        @test isempty(Quantics.findallsiteinds_by_tag(sites; tag="notfound"))
+        @test T4AQuantics.findallsiteinds_by_tag(sites; tag=tag) == sites
+        @test isempty(T4AQuantics.findallsiteinds_by_tag(sites; tag="notfound"))
 
         invalid_tag = "$(tag)="
-        @test_throws "Invalid tag: $(tag)=" Quantics.findallsiteinds_by_tag(sites,
+        @test_throws "Invalid tag: $(tag)=" T4AQuantics.findallsiteinds_by_tag(sites,
             tag=invalid_tag)
 
         invalid_sites = [Index(2, "Qubit,$(tag)=1"), Index(2, "Qubit,$(tag)=1")]
-        @test_throws "with $(tag)=1!" Quantics.findallsiteinds_by_tag(invalid_sites,
+        @test_throws "with $(tag)=1!" T4AQuantics.findallsiteinds_by_tag(invalid_sites,
             tag=tag)
-        @test_throws "Invalid tag: $(tag)=" Quantics.findallsiteinds_by_tag(invalid_sites,
+        @test_throws "Invalid tag: $(tag)=" T4AQuantics.findallsiteinds_by_tag(invalid_sites,
             tag="$(tag)=")
     end
 end
diff --git a/test/transformer_tests.jl b/test/transformer_tests.jl
index cd525c2..1469282 100644
--- a/test/transformer_tests.jl
+++ b/test/transformer_tests.jl
@@ -1,16 +1,15 @@
 
 @testitem "transformer_tests.jl/functions" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
-    using ITensors.SiteTypes: siteinds
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps, TensorTrain, siteinds, apply, product
     using LinearAlgebra
 
     @testset "upper_lower_triangle" for upper_or_lower in [:upper, :lower]
         R = 3
-        sites = siteinds("Qubit", R)
-        trimat = Quantics.upper_lower_triangle_matrix(sites, 1.0;
+        sites = [Index(2, "Qubit, n=$n") for n in 1:R]
+        trimat = T4AQuantics.upper_lower_triangle_matrix(sites, 1.0;
             upper_or_lower=upper_or_lower)
         trimatdata = Array(reduce(*, trimat), [reverse(sites')..., reverse(sites)...])
         trimatdata = reshape(trimatdata, 2^R, 2^R)
@@ -23,13 +22,13 @@
 
     @testset "cusum" begin
         R = 3
-        sites = siteinds("Qubit", R)
-        UT = Quantics.upper_lower_triangle_matrix(sites, 1.0; upper_or_lower=:lower)
-        f = Quantics.expqtt(sites, -1.0)
+        sites = [Index(2, "Qubit, n=$n") for n in 1:R]
+        UT = T4AQuantics.upper_lower_triangle_matrix(sites, 1.0; upper_or_lower=:lower)
+        f = T4AQuantics.expqtt(sites, -1.0)
         f_values = vec(Array(reduce(*, f), reverse(sites)))
         xs = collect(LinRange(0, 1, 2^R + 1)[1:(end - 1)])
 
-        g = apply(UT, f)
+        g = apply(UT, f; nsweeps=4)
         g_values = vec(Array(reduce(*, g), reverse(sites)))
 
         g_values_ref = cumsum(f_values) .- f_values # Second term remove the own values
@@ -38,11 +37,11 @@
     end
 
     @testset "flipop" for nbit in 2:3, rev_carrydirec in [true, false], bc in [1, -1]
-        sites = siteinds("Qubit", nbit)
+        sites = [Index(2, "Qubit, n=$n") for n in 1:nbit]
 
-        g = randomMPS(rev_carrydirec ? sites : reverse(sites))
+        g = random_mps(rev_carrydirec ? sites : reverse(sites))
 
-        op = Quantics.flipop(siteinds(g); rev_carrydirec=rev_carrydirec, bc=bc)
+        op = T4AQuantics.flipop(Vector{Index{Int}}(collect(Iterators.flatten(siteinds(g)))); rev_carrydirec=rev_carrydirec, bc=bc)
         f = apply(op, g; alg="naive")
         g_reconst = vec(Array(reduce(*, g), reverse(sites)))
         f_reconst = vec(Array(reduce(*, f), reverse(sites)))
@@ -59,18 +58,18 @@ end
 
 @testitem "transformer_tests.jl/reverseaxis" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps, apply
     using LinearAlgebra
 
     @testset "reverseaxis" for bc in [1], nbit in 2:2, rev_carrydirec in [true, false]
         sitesx = [Index(2, "x=$x") for x in 1:nbit]
         sites = rev_carrydirec ? sitesx : reverse(sitesx)
 
-        g = randomMPS(sites)
+        g = random_mps(sites)
 
-        f = Quantics.reverseaxis(g; tag="x", alg="naive", bc=bc)
+        f = T4AQuantics.reverseaxis(g; tag="x", alg="naive", bc=bc)
         g_reconst = vec(Array(reduce(*, g), reverse(sitesx)))
         f_reconst = vec(Array(reduce(*, f), reverse(sitesx)))
 
@@ -95,7 +94,7 @@ end
             sites = collect(Iterators.flatten(zip(reverse(sitesx), reverse(sitesy))))
         end
 
-        g = randomMPS(sites)
+        g = random_mps(sites)
 
         function _reconst(M)
             arr = Array(reduce(*, M), [reverse(sitesx)..., reverse(sitesy)...])
@@ -104,10 +103,10 @@ end
 
         g_reconst = _reconst(g)
 
-        fx = Quantics.reverseaxis(g; tag="x", alg="naive")
+        fx = T4AQuantics.reverseaxis(g; tag="x", alg="naive")
         fx_reconst = _reconst(fx)
 
-        fy = Quantics.reverseaxis(g; tag="y", alg="naive")
+        fy = T4AQuantics.reverseaxis(g; tag="y", alg="naive")
         fy_reconst = _reconst(fy)
 
         fx_ref = similar(fx_reconst)
@@ -136,7 +135,7 @@ end
                 reverse(sitesz))))
         end
 
-        g = randomMPS(sites)
+        g = random_mps(sites)
 
         function _reconst(M)
             arr = Array(reduce(*, M),
@@ -146,13 +145,13 @@ end
 
         g_reconst = _reconst(g)
 
-        fx = Quantics.reverseaxis(g; tag="x", alg="naive")
+        fx = T4AQuantics.reverseaxis(g; tag="x", alg="naive")
         fx_reconst = _reconst(fx)
 
-        fy = Quantics.reverseaxis(g; tag="y", alg="naive")
+        fy = T4AQuantics.reverseaxis(g; tag="y", alg="naive")
         fy_reconst = _reconst(fy)
 
-        fz = Quantics.reverseaxis(g; tag="z", alg="naive")
+        fz = T4AQuantics.reverseaxis(g; tag="z", alg="naive")
         fz_reconst = _reconst(fz)
 
         fx_ref = similar(fx_reconst)
@@ -178,7 +177,7 @@ end
         M = randomMPS(sites)
 
         for which_new in ["left", "right"]
-            Mnew = Quantics.asdiagonal(M, sites′; tag="n", which_new=which_new)
+            Mnew = T4AQuantics.asdiagonal(M, sites′; tag="n", which_new=which_new)
 
             M_reconst = reshape(Array(reduce(*, M), reverse(sites)), 2^R)
             Mnew_reconst = reshape(Array(reduce(*, Mnew),
@@ -193,9 +192,9 @@ end
 
 @testitem "transformer_tests.jl/phase_rotation" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps, apply
     using LinearAlgebra
 
     @testset "phase_rotation" begin
@@ -205,30 +204,30 @@ end
         sites = [Index(2, "Qubit,x=$x") for x in 1:nqbit]
         _reconst(x) = vec(Array(reduce(*, x), reverse(sites)))
 
-        f = randomMPS(sites)
+        f = random_mps(sites)
         f_vec = _reconst(f)
 
         ref = exp.(im * θ * xvec) .* f_vec
 
-        @test ref ≈ _reconst(Quantics.phase_rotation(f, θ; tag="x"))
-        @test ref ≈ _reconst(Quantics.phase_rotation(f, θ; targetsites=sites))
+        @test ref ≈ _reconst(T4AQuantics.phase_rotation(f, θ; tag="x"))
+        @test ref ≈ _reconst(T4AQuantics.phase_rotation(f, θ; targetsites=sites))
     end
 end
 
 @testitem "transformer_tests.jl/shiftaxis" begin
     using Test
-    import Quantics
+    import T4AQuantics
     using ITensors
-    using ITensorMPS: randomMPS
+    import T4AITensorCompat: random_mps, apply
     using LinearAlgebra
 
     @testset "shiftaxis" for R in [3], bc in [1, -1], rev_carrydirec in [true, false]
         sitesx = [Index(2, "Qubit, x=$n") for n in 1:R]
         sites = rev_carrydirec ? sitesx : reverse(sitesx)
-        g = randomMPS(sites)
+        g = random_mps(sites)
 
         for shift in [0, 1, 2, 2^R - 1]
-            f = Quantics.shiftaxis(g, shift; bc=bc, tag="x")
+            f = T4AQuantics.shiftaxis(g, shift; bc=bc, tag="x")
 
             f_vec = vec(Array(reduce(*, f), reverse(sitesx)))
             g_vec = vec(Array(reduce(*, g), reverse(sitesx)))
@@ -252,10 +251,10 @@ end
         else
             sites = collect(Iterators.flatten(zip(reverse(sitesx), reverse(sitesy))))
         end
-        g = randomMPS(sites)
+        g = random_mps(sites)
 
         for shift in [-4^R + 1, -1, 0, 1, 2^R - 1, 2^R, 2^R + 1, 4^R + 1]
-            f = Quantics.shiftaxis(g, shift; tag="x", bc=bc)
+            f = T4AQuantics.shiftaxis(g, shift; tag="x", bc=bc)
 
             f_mat = reshape(Array(reduce(*, f), vcat(reverse(sitesx), reverse(sitesy))),
                 2^R, 2^R)
diff --git a/test/util_tests.jl b/test/util_tests.jl
index 389d33e..f26cab6 100644
--- a/test/util_tests.jl
+++ b/test/util_tests.jl
@@ -1,21 +1,35 @@
 @testitem "util.jl" begin
     using Test
-    import PartitionedMPSs: PartitionedMPSs, SubDomainMPS, PartitionedMPS, project,
+    import T4APartitionedMPSs: T4APartitionedMPSs, SubDomainMPS, PartitionedMPS, project,
                             isprojectedat
-    import Quantics
+    import T4AQuantics
     using ITensors
-    using ITensors.SiteTypes: siteinds
-    using ITensorMPS: randomMPS, randomMPO, random_mps, MPO, MPS
+    import T4AITensorCompat: random_mps, random_mpo, TensorTrain, siteinds, linkind, linkinds
 
     include("_util.jl")
 
     @testset "_replace_mpo_siteinds!" begin
         nbit = 3
-        sites = siteinds("Qubit", nbit)
-        M = MPO(ComplexF64, sites, ["Y" for n in 1:nbit])
+        sites = [Index(2, "Qubit, n=$n") for n in 1:nbit]
+        # Create MPO manually with Y operators
+        sites_pairs = [[s, s'] for s in sites]
+        M = TensorTrain(ComplexF64, sites_pairs, 1)
+        # Set Y operator on each site
+        Y = [0 -im; im 0]
+        links = linkinds(M)
+        for n in 1:nbit
+            inds_list = [sites[n]', sites[n]]
+            if n > 1
+                push!(inds_list, links[n-1])
+            end
+            if n < nbit
+                push!(inds_list, links[n])
+            end
+            M[n] = ITensor(Y, inds_list...)
+        end
 
         sites2 = [Index(2, "n=$n") for n in 1:nbit]
-        Quantics._replace_mpo_siteinds!(M, sites, sites2)
+        T4AQuantics._replace_mpo_siteinds!(M, sites, sites2)
 
         @test all([!hasind(M[n], sites[n]) for n in 1:nbit])
         @test all([!hasind(M[n], sites[n]') for n in 1:nbit])
@@ -31,7 +45,7 @@
         csites = [Index(4, "csite=$s") for s in 1:2]
         M = randomMPS(sites; linkdims=2)
 
-        Mc = Quantics.combinesiteinds(M, csites; targetsites=sites[2:5])
+        Mc = T4AQuantics.combinesiteinds(M, csites; targetsites=sites[2:5])
 
         @test length(Mc) == 4
         @test all(dim.(siteinds(Mc)) .== [2, 4, 4, 2])
@@ -42,10 +56,10 @@
         csites = [Index(4, "csite=$s") for s in 1:(nbit ÷ 2)]
         D = 3
         mps = randomMPS(csites; linkdims=D)
-        mps_split = Quantics.splitsiteind(mps, sites)
+        mps_split = T4AQuantics.splitsiteind(mps, sites)
         @test vec(Array(reduce(*, mps_split), sites)) ≈ vec(Array(reduce(*, mps), csites))
 
-        mps_reconst = Quantics.combinesiteinds(mps_split, csites)
+        mps_reconst = T4AQuantics.combinesiteinds(mps_split, csites)
         @test vec(Array(reduce(*, mps_reconst), csites)) ≈
               vec(Array(reduce(*, mps), csites))
     end
@@ -55,14 +69,14 @@
         sites = [Index(2^R, "csite=$s") for s in 1:nsites]
 
         bonddim = 3
-        mps = randomMPS(sites; linkdims=bonddim)
+        mps = random_mps(sites; linkdims=bonddim)
 
         newsites = [[Index(2, "n=$n,m=$m") for m in 1:R] for n in 1:nsites]
 
-        mps_split = Quantics.unfuse_siteinds(mps, sites, newsites)
+        mps_split = T4AQuantics.unfuse_siteinds(mps, sites, newsites)
 
         newsites_flatten = collect(Iterators.flatten(newsites))
-        @test newsites_flatten == siteinds(mps_split)
+        @test newsites_flatten == collect(Iterators.flatten(siteinds(mps_split)))
         @test vec(Array(reduce(*, mps_split), newsites_flatten)) ≈
               vec(Array(reduce(*, mps), sites))
     end
@@ -71,7 +85,7 @@
         nsite = 6
         sites = [Index(2, "Qubit, site=$n") for n in 1:nsite]
         tensor = randomITensor(sites)
-        tensors = Quantics.split_tensor(tensor, [sites[1:2], sites[3:4], sites[5:6]])
+        tensors = T4AQuantics.split_tensor(tensor, [sites[1:2], sites[3:4], sites[5:6]])
         @test tensor ≈ reduce(*, tensors)
     end
 
@@ -79,7 +93,7 @@
         nsite = 8
         sites = [Index(2, "Qubit, site=$n") for n in 1:nsite]
         tensor = randomITensor(sites)
-        tensors = Quantics.split_tensor(tensor, [sites[1:3], sites[4:5], sites[6:8]])
+        tensors = T4AQuantics.split_tensor(tensor, [sites[1:3], sites[4:5], sites[6:8]])
         @test length(inds(tensors[1])) == 4
         @test length(inds(tensors[2])) == 4
         @test length(inds(tensors[3])) == 4
@@ -92,9 +106,9 @@
 
         sites = [Index(physdim, "n=$n") for n in 1:(2N)]
         sites_sub = sites[1:2:end]
-        M = randomMPS(sites_sub) + randomMPS(sites_sub)
+        M = random_mps(sites_sub) + random_mps(sites_sub)
 
-        M_ext = Quantics.matchsiteinds(M, sites)
+        M_ext = T4AQuantics.matchsiteinds(M, sites)
 
         tensor = Array(reduce(*, M), sites_sub)
         tensor_reconst = zeros(Float64, fill(physdim, 2N)...)
@@ -111,9 +125,9 @@
         sites = [Index(physdim, "n=$n") for n in 1:(2N)]
         sites_A = sites[1:2:end]
         sites_B = sites[2:2:end]
-        M = randomMPO(sites_A) + randomMPO(sites_A)
+        M = random_mpo(sites_A) + random_mpo(sites_A)
 
-        M_ext = Quantics.matchsiteinds(M, sites)
+        M_ext = T4AQuantics.matchsiteinds(M, sites)
 
         tensor_ref = reduce(*, M) * reduce(*, [delta(s, s') for s in sites_B])
         tensor_reconst = reduce(*, M_ext)
@@ -129,9 +143,9 @@
         sites_B = sites[2:3:end]
         sites_C = sites[3:3:end]
         sites_BC = vcat(sites_B, sites_C)
-        M = randomMPO(sites_A) + randomMPO(sites_A)
+        M = random_mpo(sites_A) + random_mpo(sites_A)
 
-        M_ext = Quantics.matchsiteinds(M, sites)
+        M_ext = T4AQuantics.matchsiteinds(M, sites)
 
         tensor_ref = reduce(*, M) * reduce(*, [delta(s, s') for s in sites_BC])
         tensor_reconst = reduce(*, M_ext)
@@ -141,13 +155,14 @@
     @testset "combinsite" begin
         nrepeat = 3
         N = 3 * nrepeat
-        sites = siteinds("Qubit", N)
-        M = MPO(randomMPS(sites))
+        sites = [Index(2, "Qubit, n=$n") for n in 1:N]
+        # Create MPO directly
+        M = random_mpo(sites, 1)
         sites1 = sites[1:3:end]
         sites2 = sites[2:3:end]
         sites3 = sites[3:3:end]
         for n in 1:nrepeat
-            M = Quantics.combinesites(M, sites1[n], sites2[n])
+            M = T4AQuantics.combinesites(M, sites1[n], sites2[n])
         end
         flag = true
         for n in 1:nrepeat
@@ -158,11 +173,11 @@
     end
 
     @testset "_directprod" begin
-        sites1 = siteinds("Qubit", 2)
-        sites2 = siteinds("Qubit", 2)
-        M1 = randomMPS(sites1)
-        M2 = randomMPS(sites2)
-        M12 = Quantics._directprod(M1, M2)
+        sites1 = [Index(2, "Qubit, n=$n") for n in 1:2]
+        sites2 = [Index(2, "Qubit, n=$n") for n in 1:2]
+        M1 = random_mps(sites1)
+        M2 = random_mps(sites2)
+        M12 = T4AQuantics._directprod(M1, M2)
 
         M1_reconst = Array(reduce(*, M1), sites1)
         M2_reconst = Array(reduce(*, M2), sites2)
@@ -185,17 +200,17 @@
 
         sitesxy_fused = [[x, y] for (x, y) in zip(sitesx, sitesy)]
 
-        Ψ_fused = Quantics.rearrange_siteinds(Ψ, sitesxy_fused)
+        Ψ_fused = T4AQuantics.rearrange_siteinds(Ψ, sitesxy_fused)
 
         @test prod(Ψ) ≈ prod(Ψ_fused) # We reconstruct a full tensor, do not use it for large L
 
-        sitesxy_fused_ = siteinds(MPO(collect(Ψ_fused)))
+        sitesxy_fused_ = siteinds(Ψ_fused)
 
         for (x, y) in zip(sitesxy_fused, sitesxy_fused_)
             @test Set(x) == Set(y)
         end
 
-        Ψ_reconst = Quantics.rearrange_siteinds(Ψ_fused, [[x] for x in sitesxy])
+        Ψ_reconst = T4AQuantics.rearrange_siteinds(Ψ_fused, [[x] for x in sitesxy])
 
         @test Ψ ≈ Ψ_reconst
     end
@@ -216,17 +231,17 @@
             push!(sitesxyz_fused, [sitesz[i]])
         end
 
-        Ψ_fused = Quantics.rearrange_siteinds(Ψ, sitesxyz_fused)
+        Ψ_fused = T4AQuantics.rearrange_siteinds(Ψ, sitesxyz_fused)
 
         @test prod(Ψ) ≈ prod(Ψ_fused)
 
-        sitesxyz_fused_ = siteinds(MPO(collect(Ψ_fused)))
+        sitesxyz_fused_ = siteinds(Ψ_fused)
 
         for (x, y) in zip(sitesxyz_fused, sitesxyz_fused_)
             @test Set(x) == Set(y)
         end
 
-        Ψ_reconst = Quantics.rearrange_siteinds(Ψ_fused, [[x] for x in sitesxyz])
+        Ψ_reconst = T4AQuantics.rearrange_siteinds(Ψ_fused, [[x] for x in sitesxyz])
 
         @test Ψ ≈ Ψ_reconst
     end
@@ -238,7 +253,7 @@
 
         Ψ = random_mps(sitesx)
 
-        M = Quantics.makesitediagonal(Ψ, "x")
+        M = T4AQuantics.makesitediagonal(Ψ, "x")
 
         Ψ_recost = Array(prod(Ψ), sitesx...)
         M_recost = Array(prod(M), prime.(sitesx)..., sitesx...)
@@ -260,7 +275,7 @@
             sitesz = [Index(2, "z=$n") for n in 1:N]
             sites = collect(collect.(zip(sitesx, sitesy, sitesz)))
 
-            Ψ = MPS(collect(_random_mpo(sites)))
+            Ψ = _random_mpo(sites)
 
             prjΨ = SubDomainMPS(Ψ)
             prjΨ1 = project(prjΨ, Dict(sitesx[1] => 1))
@@ -271,10 +286,10 @@
                 push!(sites_rearranged, sitesxy[i])
                 push!(sites_rearranged, [sitesz[i]])
             end
-            prjΨ1_rearranged = Quantics.rearrange_siteinds(prjΨ1, sites_rearranged)
+            prjΨ1_rearranged = T4AQuantics.rearrange_siteinds(prjΨ1, sites_rearranged)
 
-            @test reduce(*, MPS(prjΨ1)) ≈ reduce(*, MPS(prjΨ1_rearranged))
-            @test PartitionedMPSs.siteinds(prjΨ1_rearranged) == sites_rearranged
+            @test reduce(*, TensorTrain(prjΨ1)) ≈ reduce(*, TensorTrain(prjΨ1_rearranged))
+            @test T4APartitionedMPSs.siteinds(prjΨ1_rearranged) == sites_rearranged
         end
 
         @testset "makesitediagonal and extractdiagonal" begin
@@ -287,18 +302,18 @@
             sitesz_vec = [[z] for z in sitesz]
             sites = [x for pair in zip(sitesxy_vec, sitesz_vec) for x in pair]
 
-            Ψ = MPS(collect(_random_mpo(sites)))
+            Ψ = _random_mpo(sites)
 
             prjΨ = SubDomainMPS(Ψ)
             prjΨ1 = project(prjΨ, Dict(sitesx[1] => 1))
 
-            prjΨ1_diagonalz = Quantics.makesitediagonal(prjΨ1, "y")
-            sites_diagonalz = Iterators.flatten(siteinds(prjΨ1_diagonalz))
+            prjΨ1_diagonalz = T4AQuantics.makesitediagonal(prjΨ1, "y")
+            sites_diagonalz = Iterators.flatten(T4APartitionedMPSs.siteinds(prjΨ1_diagonalz))
 
             psi_diag = prod(prjΨ1_diagonalz.data)
             psi = prod(prjΨ1.data)
 
-            @test Quantics.extractdiagonal(prjΨ1_diagonalz, "y") ≈ prjΨ1
+            @test T4AQuantics.extractdiagonal(prjΨ1_diagonalz, "y") ≈ prjΨ1
 
             for indval in eachindval(sites_diagonalz...)
                 ind = first.(indval)