diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 77e2b0b..33c160d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -24,7 +24,7 @@ jobs: matrix: version: - '1.10' - - '1.7' + - '1.12' os: - ubuntu-latest arch: diff --git a/Project.toml b/Project.toml index 66eadd3..8ad94a0 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d" FiniteDifferenceMatrices = "a7a66f33-e7b8-47af-b618-f9b5bea05f3d" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Lux = "b2108857-7c20-44ae-9111-449ecde12c47" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2" @@ -17,6 +18,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e" +TensorTrainNumerics = "304e70ee-b558-420e-aed5-de4d05b6ca2b" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] @@ -29,5 +31,6 @@ Optimisers = "0.2, 0.3, 0.4" QuadGK = "2.11" SpecialFunctions = "2.3.1" Subscripts = "0.1.3" +TensorTrainNumerics = "1.1.4" Zygote = "0.6, 0.7" -julia = "1.7" +julia = "1.10" diff --git a/README.md b/README.md index 89bd1fb..ed8b1a1 100644 --- a/README.md +++ b/README.md @@ -18,11 +18,15 @@ flowchart TD A["Hamiltonian.jl"] C["Rayleigh-Ritz.jl"] F["FDM.jl"] + Q["QTT.jl"] N["VNN.jl"] G["VMC.jl"] H["DB.jl"] Z["TwoBody.jl"] A --> H + A --> C & F & Q & G + H --> C & F & Q & G + C & F & Q & G --> Z A --> C & F & N & G H --> C & F & N & G F --> N diff --git a/docs/make.jl b/docs/make.jl index 7479707..e3bffd7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -22,6 +22,7 @@ makedocs(; "Rayleigh-Ritz Method" => "Rayleigh-Ritz.md", "Free Complement Method" => "Free-Complement.md", "Finite Difference Method" => "FDM.md", + "Quantics Tensor Train" => "QTT.md", "Variational Neural Network" => "VNN.md", "Variational Monte Carlo" => "VMC.md", "API reference" => "API.md", diff --git a/docs/src/QTT.md b/docs/src/QTT.md new file mode 100644 index 0000000..153cc84 --- /dev/null +++ b/docs/src/QTT.md @@ -0,0 +1,160 @@ +```@meta +CurrentModule = TwoBody +``` + +# Quantics Tensor Train + +The quantics tensor train (QTT) solver rewrites a radial grid of ``N=2^q`` points as +``q`` binary sites and approximates the sampled wave function by a chain of tensor +cores, + +```math +f(r_i)=F_{b_1\ldots b_q} +\simeq G^{(1)}(b_1)\cdots G^{(q)}(b_q), +\qquad +i-1=\sum_{k=1}^{q}b_k2^{q-k}. +``` + +The Hamiltonian is stored as a matrix product operator (MPO). Two-site DMRG finds +its low-energy states, with projector penalties exposing successive excitations, + +```math +H^{(n)}_{\mathrm{QTT}}u_n=E_nu_n, +\qquad +H^{(n)}_{\mathrm{QTT}} +=H_{\mathrm{QTT}}+\mu\sum_{k=0}^{n-1}|u_k\rangle\langle u_k|. +``` + +## Theory + +Each vector core has dimensions +``G^{(k)}\in\mathbb{R}^{\chi_{k-1}\times2\times\chi_k}``, with +``\chi_0=\chi_q=1``. A matrix is quantized in both its row and column indices and +represented as + +```math +A_{b_1\ldots b_q,\,c_1\ldots c_q} +\simeq W^{(1)}(b_1,c_1)W^{(2)}(b_2,c_2)\cdots W^{(q)}(b_q,c_q). +``` + +The central second-difference operator can be written using the one-point shift +matrices ``S_-`` and ``S_+`` as + +```math +D^{(2)} = \frac{S_- - 2I + S_+}{\Delta r^2}. +``` + +This tridiagonal operator has an exact QTT/MPO representation with maximum bond +dimension three. Potential functions and initial wave functions are compressed by +tensor cross interpolation. Two-site DMRG sweeps optimize the tensor cores and adapt +their bond dimensions; intermediate MPOs are compressed between solves. + +If ``\chi`` and ``\rho`` bound the vector and MPO bond dimensions, their storage is +bounded by + +```math +\operatorname{storage}(F_{\mathrm{QTT}}) \leq 2q\chi^2, +\qquad +\operatorname{storage}(A_{\mathrm{QTT}}) \leq 4q\rho^2. +``` + +Thus, when the QTT ranks remain moderate, storage grows as ``O(\log N)`` instead of +``O(N)`` for a dense vector or ``O(N^2)`` for a dense matrix. + +## Usage + +For the three-dimensional spherical oscillator in atomic units, configure the QTT +grid with `quantics`; the number of interior grid points is `2^quantics`. Here a +coarse grid is used so the example runs quickly. + +```@example qtt +using TwoBody +using Random + +Random.seed!(1234) + +H = Hamiltonian( + Kinetic(hbar=1, m=1), + PowerLaw(coefficient=1 / 2, exponent=2), +) + +method = QuanticsTensorTrainMethod( + quantics=5, + r₀=0.0, + rₘₐₓ=8.0, + l=0, + tolerance=1e-8, + maxbonddim=16, + maxoperatorbonddim=32, + sweeps=4, +) + +result = solve( + H, + method; + initial=r -> exp(-r^2 / 2), + nₘₐₓ=2, + info=0, +) +result.E +``` + +The values approach the exact ``l=0`` oscillator energies ``E_0=3/2`` and +``E_1=7/2`` as the grid is refined. Inspect the eigenvector ranks without expanding +them: + +```@example qtt +ranks.(result.C) +``` + +The overlap matrix and residual norms provide convergence diagnostics: + +```@example qtt +result.overlaps +``` + +```@example qtt +result.residuals +``` + +Individual entries can also be contracted directly: + +```@example qtt +qttvalue(result.ψ[1], 1) +``` + +!!! note "Current scope" + `QuanticsTensorTrainMethod` currently supports `Kinetic`, `RestEnergy`, + `Constant`, `Linear`, `Coulomb`, `PowerLaw`, `Gaussian`, `Exponential`, `Yukawa`, + and `Custom` operators. Choose `deflation_shift` above the relevant spectral gaps, + monitor MPO ranks, overlaps, and residuals, and perform a grid-convergence study. + +## Acknowledgments + +The proof of concept for this QTT solver was developed at +[CompPhysHack 2026](https://qc-hybrid.github.io/CompPhysHack2026/) in collaboration +with Lucas Arenstein; its source is available in the +[CompPhysHack repository](https://github.com/ohno/CompPhysHack2026Ohno/blob/main/julia/qtt.jl). +The current solver builds on the tensor-train operations and two-site DMRG +eigensolver in +[TensorTrainNumerics.jl](https://github.com/MartinMikkelsen/TensorTrainNumerics.jl). +We thank Lucas Arenstein, the TensorTrainNumerics.jl developers, and the hackathon +organizers for their contributions and support. + +## Bibliography + +- L. Arenstein, M. Mikkelsen, and M. Kastoryano, + [*Fast and Flexible Quantum-Inspired Differential Equation Solvers with Data Integration*](https://doi.org/10.48550/arXiv.2505.17046), + arXiv:2505.17046 (2025). + +## API reference + +```@docs; canonical=false +QuanticsTensorTrainMethod +solve(hamiltonian::Hamiltonian, method::QuanticsTensorTrainMethod) +QTTVector +QTTMatrix +order +ranks +qttvalue +``` diff --git a/src/QTT.jl b/src/QTT.jl new file mode 100644 index 0000000..c15f1f2 --- /dev/null +++ b/src/QTT.jl @@ -0,0 +1,337 @@ +export QuanticsTensorTrainMethod, QTTVector, QTTMatrix, order, ranks, qttvalue + +import LinearAlgebra +import Logging +import TensorTrainNumerics + +"A vector stored as a quantics tensor train." +const QTTVector = TensorTrainNumerics.TTvector + +"A matrix stored as a quantics tensor train / matrix product operator." +const QTTMatrix = TensorTrainNumerics.TToperator + +"Return the number of binary sites in a QTT." +order(tt::Union{QTTVector,QTTMatrix}) = tt.N + +"Return the left boundary, internal, and right boundary ranks of a QTT." +ranks(tt::QTTVector) = copy(tt.ttv_rks) +ranks(tt::QTTMatrix) = copy(tt.tto_rks) + +struct QuanticsTensorTrainMethod + quantics::Int + r₀::Float64 + rₘₐₓ::Float64 + R::LinRange{Float64,Int} + Δr::Float64 + l::Int + tolerance::Float64 + maxbonddim::Int + maxoperatorbonddim::Int + maxiter::Int + sweeps::Int + deflation_shift::Float64 + + function QuanticsTensorTrainMethod(; + quantics=10, + r₀=0.0, + rₘₐₓ=50.0, + l=0, + tolerance=1e-10, + maxbonddim=64, + maxoperatorbonddim=2 * maxbonddim, + maxiter=100, + sweeps=4, + deflation_shift=10.0, + ) + quantics >= 2 || throw(ArgumentError("quantics must be at least 2")) + 0 <= r₀ < rₘₐₓ || throw(ArgumentError("the radial interval must satisfy 0 <= r₀ < rₘₐₓ")) + l >= 0 || throw(ArgumentError("l must be nonnegative")) + tolerance > 0 || throw(ArgumentError("tolerance must be positive")) + maxbonddim >= 1 || throw(ArgumentError("maxbonddim must be positive")) + maxoperatorbonddim >= 1 || throw(ArgumentError("maxoperatorbonddim must be positive")) + maxiter >= 1 || throw(ArgumentError("maxiter must be positive")) + sweeps >= 1 || throw(ArgumentError("sweeps must be positive")) + deflation_shift > 0 || throw(ArgumentError("deflation_shift must be positive")) + + npoints = 1 << quantics + Δr = Float64(rₘₐₓ - r₀) / (npoints + 1) + R = range(Float64(r₀) + Δr, Float64(rₘₐₓ) - Δr; length=npoints) + new( + quantics, + Float64(r₀), + Float64(rₘₐₓ), + R, + Δr, + l, + Float64(tolerance), + maxbonddim, + maxoperatorbonddim, + maxiter, + sweeps, + Float64(deflation_shift), + ) + end +end + +Base.string(method::QuanticsTensorTrainMethod) = + "QuanticsTensorTrainMethod(" * + join(["$(symbol)=$(getproperty(method, symbol))" for symbol in + (:quantics, :r₀, :rₘₐₓ, :l, :tolerance, :maxbonddim, + :maxoperatorbonddim, :maxiter, :sweeps, :deflation_shift)], ", ") * ")" +Base.show(io::IO, method::QuanticsTensorTrainMethod) = print(io, string(method)) + +function qttvalue(tt::QTTVector{T}, index::Integer) where {T} + n = 1 << order(tt) + 1 <= index <= n || throw(BoundsError(1:n, index)) + state = ones(T, 1) + index0 = index - 1 + for site in 1:order(tt) + physical = Int((index0 >> (order(tt) - site)) & 1) + 1 + state = vec(reshape(state, 1, :) * @view(tt.ttv_vec[site][physical, :, :])) + end + return only(state) +end + +function qttvalue(tt::QTTMatrix{T}, row::Integer, column::Integer) where {T} + n = 1 << order(tt) + 1 <= row <= n || throw(BoundsError(1:n, row)) + 1 <= column <= n || throw(BoundsError(1:n, column)) + state = ones(T, 1) + row0, column0 = row - 1, column - 1 + for site in 1:order(tt) + rowbit = Int((row0 >> (order(tt) - site)) & 1) + 1 + columnbit = Int((column0 >> (order(tt) - site)) & 1) + 1 + state = vec(reshape(state, 1, :) * @view(tt.tto_vec[site][rowbit, columnbit, :, :])) + end + return only(state) +end + +function _qttvector(f::Function, method::QuanticsTensorTrainMethod) + q = method.quantics + domain = [Float64[0, 1] for _ in 1:q] + + function evaluate(bits::AbstractMatrix) + indices = zeros(Int, size(bits, 1)) + for site in 1:q + indices .+= round.(Int, @view(bits[:, site])) .* (1 << (q - site)) + end + return [Float64(f(method.R[index + 1])) for index in indices] + end + + algorithm = TensorTrainNumerics.MaxVol( + tol=method.tolerance, + rmax=method.maxbonddim, + maxiter=method.maxiter, + verbose=false, + ) + return TensorTrainNumerics.tt_cross( + evaluate, + domain, + algorithm; + ranks=min(2, method.maxbonddim), + val_size=min(1000, 1 << q), + ) +end + +function _compress_vector(vector::QTTVector, method::QuanticsTensorTrainMethod) + compressed = copy(vector) + TensorTrainNumerics.tt_compress!( + compressed, + method.maxbonddim; + truncerr=method.tolerance, + ) + return compressed +end + +function _compress_operator(operator::QTTMatrix, method::QuanticsTensorTrainMethod) + fused = TensorTrainNumerics.tto_to_ttv(operator) + TensorTrainNumerics.tt_compress!( + fused, + method.maxoperatorbonddim; + truncerr=method.tolerance, + ) + return TensorTrainNumerics.ttv_to_tto(fused) +end + +function matrix(term::Kinetic, method::QuanticsTensorTrainMethod) + kinetic = (term.hbar^2 / (2 * term.m * method.Δr^2)) * + TensorTrainNumerics.Δ(method.quantics) + if method.l == 0 + return kinetic + end + centrifugal = _qttvector( + r -> term.hbar^2 * method.l * (method.l + 1) / (2 * term.m * r^2), + method, + ) + return kinetic + TensorTrainNumerics.ttv_to_diag_tto(centrifugal) +end + +matrix(term::RestEnergy, method::QuanticsTensorTrainMethod) = + TensorTrainNumerics.ttv_to_diag_tto(_qttvector(_ -> term.m * term.c^2, method)) + +const QTTPotentialTerm = Union{Constant,Linear,Coulomb,PowerLaw,Gaussian,Exponential,Yukawa,Custom} + +matrix(term::QTTPotentialTerm, method::QuanticsTensorTrainMethod) = + TensorTrainNumerics.ttv_to_diag_tto(_qttvector(r -> V(term, r), method)) + +matrix(term::PotentialTerm, ::QuanticsTensorTrainMethod) = + throw(ArgumentError("$(typeof(term)) is not supported by QuanticsTensorTrainMethod")) + +function matrix(hamiltonian::Hamiltonian, method::QuanticsTensorTrainMethod) + isempty(hamiltonian.terms) && throw(ArgumentError("the Hamiltonian must contain at least one term")) + matrices = map(hamiltonian.terms) do term + term isa Union{Kinetic,RestEnergy,QTTPotentialTerm} || + throw(ArgumentError("$(typeof(term)) is not supported by QuanticsTensorTrainMethod")) + matrix(term, method) + end + return _compress_operator(reduce(+, matrices), method) +end + +_ttnorm(vector::QTTVector) = TensorTrainNumerics.norm(vector) +_ttdot(left::QTTVector, right::QTTVector) = TensorTrainNumerics.dot(left, right) + +function _dmrg_groundstate(operator::QTTMatrix, guess::QTTVector, method::QuanticsTensorTrainMethod) + normalized_guess = guess / _ttnorm(guess) + history, state, rank_history = Logging.with_logger(Logging.NullLogger()) do + TensorTrainNumerics.dmrg_eigsolve( + operator, + normalized_guess; + N=2, + tol=method.tolerance, + sweep_schedule=[method.sweeps], + rmax_schedule=[method.maxbonddim], + it_solver=true, + linsolv_maxiter=method.maxiter, + linsolv_tol=max(sqrt(method.tolerance), 1e-8), + ) + end + state = state / _ttnorm(state) + return last(history), state, history, rank_history +end + +function _state_guess(initial::Function, level::Int, method::QuanticsTensorTrainMethod) + width = method.rₘₐₓ - method.r₀ + modulation(r) = level == 0 ? 1.0 : cos(level * π * (r - method.r₀) / width) + return _qttvector(r -> r * initial(r) * modulation(r), method) +end + +function _solve_qtt( + hamiltonian::Hamiltonian, + initial::Function, + method::QuanticsTensorTrainMethod; + info=4, + nₘₐₓ=4, +) + 1 <= nₘₐₓ <= (1 << method.quantics) || + throw(ArgumentError("nₘₐₓ must lie between 1 and the number of grid points")) + + H = matrix(hamiltonian, method) + deflated = H + energies = Float64[] + states = QTTVector[] + histories = Vector{Float64}[] + rank_histories = Vector{Int}[] + residuals = Float64[] + + for level in 0:(nₘₐₓ - 1) + guess = _state_guess(initial, level, method) + _, state, history, rank_history = _dmrg_groundstate(deflated, guess, method) + + Hstate = H * state + energy = real(_ttdot(state, Hstate) / _ttdot(state, state)) + residual = _ttnorm(Hstate - energy * state) + push!(energies, energy) + push!(states, state) + push!(histories, history) + push!(rank_histories, rank_history) + push!(residuals, residual) + + if level < nₘₐₓ - 1 + projector = TensorTrainNumerics.outer_product(state, state) + deflated = _compress_operator(deflated + method.deflation_shift * projector, method) + end + end + + overlaps = [_ttdot(states[i], states[j]) for i in eachindex(states), j in eachindex(states)] + scale = inv(sqrt(4π * method.Δr)) + u = [scale * state for state in states] + inverse_radius = _qttvector(r -> inv(r), method) + ψ = [_compress_vector(TensorTrainNumerics.hadamard(state, inverse_radius), method) for state in u] + + if info > 0 + println("\n# method\n") + println(method) + println("\n# eigenvalues\n") + for level in 1:min(nₘₐₓ, info) + println("E[$level] = $(energies[level])") + end + println() + end + + return info >= 0 ? ( + hamiltonian=hamiltonian, + method=method, + nₘₐₓ=nₘₐₓ, + H=H, + E=energies, + C=states, + ψ=ψ, + u=u, + overlaps=overlaps, + residuals=residuals, + histories=histories, + rank_histories=rank_histories, + ) : (E=energies,) +end + +function solve( + hamiltonian::Hamiltonian, + method::QuanticsTensorTrainMethod; + initial=r -> exp(-r), + info=4, + nₘₐₓ=4, +) + return _solve_qtt(hamiltonian, initial, method; info=info, nₘₐₓ=nₘₐₓ) +end + +function solve( + hamiltonian::Hamiltonian, + wavefunction::Function, + method::QuanticsTensorTrainMethod; + info=4, + nₘₐₓ=4, +) + return _solve_qtt(hamiltonian, wavefunction, method; info=info, nₘₐₓ=nₘₐₓ) +end + +@doc raw""" + QuanticsTensorTrainMethod(; quantics=10, r₀=0.0, rₘₐₓ=50.0, l=0, + tolerance=1e-10, maxbonddim=64, + maxoperatorbonddim=128, maxiter=100, sweeps=4, + deflation_shift=10.0) + +Configure a radial QTT grid with ``2^{\mathtt{quantics}}`` uniformly spaced interior +points. The excluded endpoints `r₀` and `rₘₐₓ` impose zero-value (Dirichlet) +boundaries. `maxbonddim` and `maxoperatorbonddim` bound state and MPO ranks; +`deflation_shift` is the penalty for each previously computed state. +""" QuanticsTensorTrainMethod + +@doc raw""" + solve(hamiltonian, method::QuanticsTensorTrainMethod; + initial=r -> exp(-r), info=4, nₘₐₓ=4) + solve(hamiltonian, initial::Function, method::QuanticsTensorTrainMethod; + info=4, nₘₐₓ=4) + +Compute the lowest `nₘₐₓ` QTT eigenstates. `initial` seeds the first DMRG solve; +`info` controls progress output. The result contains energies `E`, unit-norm tensor +trains `C`, radial functions `ψ` and `u`, and the diagnostics `overlaps`, `residuals`, +`histories`, and `rank_histories`. Each energy is evaluated with the original +Hamiltonian rather than its deflated counterpart. +""" solve(hamiltonian::Hamiltonian, method::QuanticsTensorTrainMethod; initial=r -> exp(-r), info=4, nₘₐₓ=4) + +@doc raw""" + qttvalue(tt::QTTVector, index) + qttvalue(tt::QTTMatrix, row, column) + +Contract a QTT at one vector or matrix index without materializing its ``2^q`` entries. +""" qttvalue diff --git a/src/TwoBody.jl b/src/TwoBody.jl index 8e2607e..354c55f 100644 --- a/src/TwoBody.jl +++ b/src/TwoBody.jl @@ -12,6 +12,7 @@ include("./Basis.jl") # Solvers include("./Rayleigh-Ritz.jl") include("./FDM.jl") +include("./QTT.jl") include("./VNN.jl") include("./VMC.jl") diff --git a/test/QTT.jl b/test/QTT.jl new file mode 100644 index 0000000..f62102e --- /dev/null +++ b/test/QTT.jl @@ -0,0 +1,57 @@ +@testset "QTT.jl" begin + @test_throws ArgumentError QuanticsTensorTrainMethod(quantics=1) + @test_throws ArgumentError QuanticsTensorTrainMethod(r₀=-0.1) + @test_throws ArgumentError QuanticsTensorTrainMethod(deflation_shift=0) + + method = QuanticsTensorTrainMethod( + quantics=5, + r₀=0.0, + rₘₐₓ=8.0, + tolerance=1e-8, + maxbonddim=16, + maxoperatorbonddim=32, + maxiter=100, + sweeps=4, + deflation_shift=10.0, + ) + @test length(method.R) == 2^method.quantics + @test method.R[1] ≈ method.r₀ + method.Δr + @test method.R[end] ≈ method.rₘₐₓ - method.Δr + + kinetic = TwoBody.matrix(Kinetic(hbar=1, m=1), method) + for row in 1:8, column in 1:8 + expected = row == column ? 1 / method.Δr^2 : + abs(row - column) == 1 ? -1 / (2 * method.Δr^2) : 0 + @test qttvalue(kinetic, row, column) ≈ expected atol=1e-10 + end + + H = Hamiltonian( + Kinetic(hbar=1, m=1), + PowerLaw(coefficient=1 / 2, exponent=2), + ) + Random.seed!(1234) + result = solve( + H, + method; + initial=r -> exp(-r^2 / 2), + nₘₐₓ=2, + info=0, + ) + @test result.E ≈ [1.5, 3.5] atol=6e-2 + @test abs(result.overlaps[1, 2]) < 1e-6 + @test result.residuals[1] < 1e-5 + @test result.residuals[2] < 1e-5 + @test length(result.ψ) == length(result.u) == length(result.C) == 2 + @test all(order(state) == method.quantics for state in result.C) + @test all(first(ranks(state)) == last(ranks(state)) == 1 for state in result.C) + @test qttvalue(result.ψ[1], 1) isa Real + + discrete_norm = 4π * method.Δr * + sum(abs2(qttvalue(result.u[1], index)) for index in eachindex(method.R)) + @test discrete_norm ≈ 1 atol=1e-8 + + potential = TwoBody.matrix(PowerLaw(coefficient=1 / 2, exponent=2), method) + @test qttvalue(potential, 5, 5) ≈ method.R[5]^2 / 2 atol=1e-7 + @test qttvalue(potential, 5, 6) ≈ 0 atol=1e-12 + @test_throws ArgumentError TwoBody.matrix(Delta(), method) +end diff --git a/test/runtests.jl b/test/runtests.jl index 2d84f92..c447c61 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ using Printf using Antique using SpecialFunctions using ForwardDiff +using Random @testset verbose = true "TwoBody.jl" begin include("Hamiltonian.jl") @@ -13,6 +14,7 @@ using ForwardDiff include("FC.jl") include("Rayleigh-Ritz.jl") include("FDM.jl") + include("QTT.jl") include("VNN.jl") include("VMC.jl") end