From 0eeb2fe2c757a5734f4f88a2aca5f50540afd462 Mon Sep 17 00:00:00 2001 From: qmortier Date: Thu, 8 Feb 2024 17:53:01 +0100 Subject: [PATCH] ctmrg (#22) * ctmrg * Create UniformPeps.m * remove plot tsvd + ctmrg tests * Add CTMRG to docs --------- Co-authored-by: Lander Burgelman <39218680+leburgel@users.noreply.github.com> Co-authored-by: leburgel --- docs/src/index.rst | 2 +- docs/src/lib/algorithms.rst | 9 + docs/src/lib/environments.rst | 3 +- docs/src/lib/mps.rst | 3 + src/algorithms/Ctmrg.m | 369 ++++++++++++++++++++++++++++ src/environments/CtmrgEnvironment.m | 64 +++++ src/mps/PepsSandwich.m | 4 +- src/mps/PepsTensor.m | 5 + src/mps/UniformPeps.m | 163 ++++++++++++ src/tensors/Tensor.m | 18 +- src/utility/plot_schmidtspectrum.m | 59 +++++ test/TestCtmrg.m | 76 ++++++ 12 files changed, 768 insertions(+), 7 deletions(-) create mode 100644 src/algorithms/Ctmrg.m create mode 100644 src/environments/CtmrgEnvironment.m create mode 100644 src/mps/UniformPeps.m create mode 100644 src/utility/plot_schmidtspectrum.m create mode 100644 test/TestCtmrg.m diff --git a/docs/src/index.rst b/docs/src/index.rst index 9a63f69..63e72a7 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -35,8 +35,8 @@ Additionally, for tensors which are invariant under general global symmetries, v lib/tensors lib/sparse lib/mps - lib/environments lib/algorithms + lib/environments lib/models lib/utility lib/caches diff --git a/docs/src/lib/algorithms.rst b/docs/src/lib/algorithms.rst index 446c01d..204c96c 100644 --- a/docs/src/lib/algorithms.rst +++ b/docs/src/lib/algorithms.rst @@ -44,6 +44,15 @@ Infinite MPS algorithms :members: changebonds +Infinite PEPS algorithms +------------------------ + +.. autoclass:: src.algorithms.Ctmrg + :no-members: + :members: fixedpoint + + + Eigsolvers ---------- diff --git a/docs/src/lib/environments.rst b/docs/src/lib/environments.rst index cde31f8..696ad17 100644 --- a/docs/src/lib/environments.rst +++ b/docs/src/lib/environments.rst @@ -9,4 +9,5 @@ Environments This section contains the API documentation for the :mod:`.environments` module. .. autoclass:: src.environments.FiniteEnvironment - +.. autoclass:: src.environments.CtmrgEnvironment + :no-members: diff --git a/docs/src/lib/mps.rst b/docs/src/lib/mps.rst index 53e7d7e..569c9d6 100644 --- a/docs/src/lib/mps.rst +++ b/docs/src/lib/mps.rst @@ -28,3 +28,6 @@ Operators .. autoclass:: src.mps.PepsTensor .. autoclass:: src.mps.PepsSandwich :no-members: +.. autoclass:: src.mps.UniformPeps + :no-members: + :members: UniformPeps diff --git a/src/algorithms/Ctmrg.m b/src/algorithms/Ctmrg.m new file mode 100644 index 0000000..2c15eeb --- /dev/null +++ b/src/algorithms/Ctmrg.m @@ -0,0 +1,369 @@ +classdef Ctmrg + % `Corner Transfer Matrix Renormalisation Group algorithm `_ for PEPS. + % + % Properties + % ---------- + % tol : :class:`double` + % tolerance for convergence criterion, defaults to :code:`1e-10`. + % + % miniter : :class:`int` + % minimum number of iteration, defaults to :code:`5`. + % + % maxiter : :class:`int` + % maximum number of iteration, defaults to :code:`100`. + % + % projectortype : :class:`char` + % projector scheme used in the algorithm, currently only supports a single default + % value. + % + % trunc : :class:`struct` + % specification of truncation options, see :meth:`.Tensor.tsvd` for details. + % + % verbosity : :class:`.Verbosity` + % verbosity level of the algorithm, defaults to :code:`Verbosity.iter`. + % + % doplot : :class:`logical` + % plot progress, defaults to :code:`false`. + + %% Options + properties + + tol = 1e-10 + miniter = 5 + maxiter = 100 + + projectortype = 'cornersboris' + trunc + + verbosity = Verbosity.iter + doPlot = false + + doSave = false + saveIterations = 1 + saveMethod = 'full' + name = 'CTMRG' + + end + + methods + function alg = Ctmrg(kwargs) + arguments + kwargs.?Ctmrg + end + + fields = fieldnames(kwargs); + if ~isempty(fields) + for field = fields.' + alg.(field{1}) = kwargs.(field{1}); + end + end + + end + + function [envs, new_norm, err] = fixedpoint(alg, peps_top, peps_bot, envs) + % Find the fixed point CTMRG environment of an infinite PEPS overlap. + % + % Usage + % ----- + % :code:`[envs, new_norm, err] = fixedpoint(alg, peps_top, peps_bot, envs)` + % + % Arguments + % --------- + % alg : :class:`.Ctmrg` + % CTMRG algorithm. + % + % peps_top : :class:`.UniformPeps` + % top-layer uniform PEPS, usually interpreted as the 'ket' in the overlap. + % + % peps_bot : :class:`.UniformPeps` + % bottom-layer uniform PEPS, usually interpreted as the 'bra' in the overlap, + % defaults to the top layer state. + % + % envs : :class:`.CtmrgEnv` + % initial guess for the fixed point CTMRG environment. + % + % Returns + % ------- + % envs : :class:`.CtmrgEnv` + % fixed point CTMRG environment. + % + % new_norm : :class:`double` + % corresponding norm. + % + % err : :class:`double` + % final error measure at convergence. + + arguments + alg + peps_top + peps_bot = peps_top + envs = CtmrgEnvironment(peps_top, peps_bot) + end + + if isempty(alg.trunc) + alg.trunc = struct('TruncSpace', space(envs.corners{1,:,:},1)); + end + + err = Inf; + iter = 1; + + t_total = tic; + disp_init(alg); + + old_norm = contract_ctrmg(alg, peps_top, peps_bot, envs); + + eta1 = 1.0; + while (err > alg.tol && iter <= alg.maxiter) || iter <= alg.miniter + + t_iter = tic; + + eta = 0.0; + for i = 1:4 + if isfield(alg.trunc,'doPlot') && alg.trunc.doPlot + subplot(2,2,i) + end + [envs, eta0] = left_move(alg, peps_top, peps_bot, envs); + eta = max(eta, eta0); + envs = rot90(envs); + peps_top = rot90(peps_top); + peps_bot = rot90(peps_bot); + end + + new_norm = contract_ctrmg(alg, peps_top, peps_bot, envs); + + err = abs(old_norm - new_norm); + deta = abs((eta1 - eta) / eta1); + + if alg.doPlot + plot(alg, iter, envs, err); + end + + disp_iter(alg, iter, err, abs(new_norm), eta, deta, toc(t_iter)); + if alg.doSave && mod(iter, alg.saveIterations) == 0 + save_iteration(alg, envs, err, new_norm, iter, eta, toc(t_total)); + end + + if iter > alg.miniter && err < alg.tol + disp_conv(alg, err, abs(new_norm), eta, deta, toc(t_total)); + return + end + + old_norm = new_norm; + eta1 = eta; + iter = iter + 1; + end + + disp_maxiter(alg, err, abs(new_norm), eta, deta, toc(t_total)); + end + + function [envs, eta] = left_move(alg, peps_top, peps_bot, envs) + + eta = 0.0; + h = height(peps_top); + w = width(peps_top); + + for j = 1:w + + [above_projs, below_projs, etaj] = get_projectors(alg, peps_top, peps_bot, envs, j); + eta = max(eta, etaj); + + %use the projectors to grow the corners/edges + for i = 1:h + + envs.corners{1,prev(i,h),j} = contract(envs.corners{1,prev(i,h),prev(j,w)}, [2 1], ... + envs.edges{1,prev(i,h),j}, [1 3 4 -2], ... + above_projs{prev(i,h)}, [-1 4 3 2], ... + 'Rank', [1,1]); + + envs.edges{4,i,j} = contract( envs.edges{4,i,prev(j,w)}, [7 2 4 1], ... + peps_top.A{i,j}.var, [6 2 8 -2 3], ... + conj(peps_bot.A{i,j}.var), [6 4 9 -3 5], ... + below_projs{prev(i,h)}, [1 3 5 -4], ... + above_projs{i}, [-1 9 8 7], ... + 'Rank', [3,1]); + + envs.corners{4,next(i,h),j} = contract(envs.corners{4,next(i,h),prev(j,w)}, [1 2], ... + envs.edges{3,next(i,h),j}, [-1 3 4 1], ... + below_projs{i}, [2 3 4 -2], ... + 'Rank', [1,1]); + + end + + envs.corners(1,:,j) = cellfun(@(x) x ./ norm(x), envs.corners(1,:,j), 'UniformOutput', false); + envs.edges(4,:,j) = cellfun(@(x) x ./ norm(x), envs.edges(4,:,j), 'UniformOutput', false); + envs.corners(4,:,j) = cellfun(@(x) x ./ norm(x), envs.corners(4,:,j), 'UniformOutput', false); + + end + end + + function [above_projs, below_projs, etaj] = get_projectors(alg, peps_top, peps_bot, envs, j) + + etaj = 0; + h = height(peps_top); + w = width(peps_top); + above_projs = cell(h,1); + below_projs = cell(h,1); + + switch alg.projectortype + + case 'cornersboris' + + for i = 1:h + + %SW corner + Q1 = contract( envs.edges{3, mod1(i+2,h), j}, [-1 5 3 1], ... + envs.corners{4, mod1(i+2,h), prev(j,w)}, [1 2], ... + envs.edges{4, next(i,h), prev(j,w)}, [2 6 4 -6], ... + peps_top.A{next(i,h), j}.var, [7 6 5 -2 -5], ... + conj(peps_bot.A{next(i,h), j}.var), [7 4 3 -3 -4], ... + 'Rank', [3,3]); + %NW corner + Q2 = contract( envs.edges{4, i, prev(j,w)}, [-1 3 5 2], ... + envs.corners{1, prev(i,h), prev(j,w)}, [2 1], ... + envs.edges{1, prev(i,h), j}, [1 4 6 -6], ... + peps_top.A{i, j}.var, [7 3 -2 -5 4], ... + conj(peps_bot.A{i, j}.var), [7 5 -3 -4 6], ... + 'Rank', [3,3]); + + trun = [fieldnames(alg.trunc),struct2cell(alg.trunc)]'; + [U, S, V, etaij] = tsvd(contract(Q1,[-1,-2,-3,3,2,1],Q2,[1,2,3,-4,-5,-6],'Rank',[3,3]), trun{:}); + + etaj = max(etaj, etaij); + isqS = inv(sqrtm(S)); + Q = isqS*U'*Q1; + P = Q2*V'*isqS; + + %norm(contract(Q,[-1,1,2,3],P,[3,2,1,-2],'Rank',[1,1])-Tensor.eye(space(Q,1),space(Q,1))) + %norm(P*Q-Tensor.eye(space(Q,2:4)',space(Q,2:4)')) + + above_projs{i} = Q; + below_projs{i} = P; %should be ok for any arrow + + end + + end + + + end + + function total = contract_ctrmg(alg, peps_top, peps_bot, envs) + + total = 1.0; + h = height(peps_top); + w = width(peps_top); + + for i = 1:h + for j = 1:w + total = total * contract( envs.corners{1,prev(i,h),prev(j,w)}, [9 1], ... + envs.edges{1,prev(i,h),j}, [1 4 7 2], ... + envs.corners{2,prev(i,h),next(j,w)}, [2 3], ... + envs.edges{2,i,next(j,w)}, [3 5 8 17], ... + envs.corners{3,next(i,h),next(j,w)}, [17 16], ... + envs.edges{3,next(i,h),j}, [16 14 15 13], ... + envs.corners{4,next(i,h),prev(j,w)}, [13 12], ... + envs.edges{4,i,prev(j,w)}, [12 10 11 9], ... + peps_top.A{i,j}.var, [6 10 14 5 4], ... + conj(peps_bot.A{i,j}.var), [6 11 15 8 7]); + + total = total * contract( envs.corners{1,prev(i,h),prev(j,w)}, [4 1], ... + envs.corners{2,prev(i,h),j}, [1 2], ... + envs.corners{3,i,j}, [2 3], ... + envs.corners{4,i,prev(j,w)}, [3 4]); + + total = total / contract( envs.corners{1,prev(i,h),prev(j,w)}, [2 1], ... + envs.corners{2,prev(i,h),j}, [1 3], ... + envs.edges{2,i,j}, [3 4 5 6], ... + envs.corners{3,next(i,h),j}, [6 7], ... + envs.corners{4,next(i,h),prev(j,w)}, [7 8], ... + envs.edges{4,i,prev(j,w)}, [8 4 5 2]); + + total = total / contract( envs.corners{1,prev(i,h),prev(j,w)}, [1 3], ... + envs.corners{4,i,prev(j,w)}, [2 1], ... + envs.edges{1,prev(i,h),j}, [3 4 5 6], ... + envs.corners{2,prev(i,h),next(j,w)}, [6 7], ... + envs.corners{3,i,next(j,w)}, [7 8], ... + envs.edges{3,i,j}, [8 4 5 2]); + end + end + + end + + end + + %% Display + methods(Access = private) + + function plot(alg, iter, envs, eta) + if ~alg.doPlot, return; end + persistent axhistory axspectrum + + D = height(envs); + W = width(envs); + + if iter == 1 + figure('Name', 'Ctmrg'); + axhistory = subplot(D + 1, W, 1:W); + axspectrum = gobjects(D, W); + for d = 1:D, for w = 1:W + axspectrum(d, w) = subplot(D + 1, W, w + d * W); + end, end + linkaxes(axspectrum, 'y'); + end + + + if isempty(axhistory.Children) + semilogy(axhistory, iter, eta, '.', 'Color', colors(1), ... + 'DisplayName', 'Errors', 'MarkerSize', 10); + hold on + ylim(axhistory, [alg.tol / 5, max([1 eta])]); + yline(axhistory, alg.tol, '-', 'Convergence', ... + 'LineWidth', 2); + hold off + else + axhistory.Children(end).XData(end+1) = iter; + axhistory.Children(end).YData(end+1) = eta; + end + plot_entanglementspectrum(UniformMps({MpsTensor(envs.edges{1,:,:})}), 1:D, 1:W, axspectrum); + drawnow + end + + function disp_init(alg) + if alg.verbosity < Verbosity.conv, return; end + fprintf('---- CTMRG ----\n'); + end + + function disp_iter(alg, iter, err, norm, eta, deta, t) + if alg.verbosity < Verbosity.iter, return; end + fprintf('iteration: %4d\t\terror: %.2e\t\tnorm: %.10e\t\teta: %.2e\t\tdeta: %.2e\t(%s)\n', ... + iter, err, norm, eta, deta, time2str(t)); + end + + function disp_conv(alg, err, norm, eta, deta, t) + if alg.verbosity < Verbosity.conv, return; end + fprintf('CTMRG converged \t\terror: %.2e\t\tnorm: %.10e\t\teta: %.2e\t\tdeta: %.2e\t(%s)\n', ... + err, norm, eta, deta, time2str(t)); + fprintf('---------------\n'); + end + + function disp_maxiter(alg, err, norm, eta, deta, t) + if alg.verbosity < Verbosity.warn, return; end + fprintf('CTMRG max iterations \terror: %.2e\t\tnorm: %.10e\t\teta: %.2e\t\tdeta: %.2e\t(%s)\n', ... + err, norm, eta, deta, time2str(t)); + fprintf('---------------\n'); + end + + function save_iteration(alg, envs, err, norm, iter, eta, t) + fileName = alg.name; + fileData = struct; + fileData.envs = envs; + fileData.err = err; + fileData.norm = norm; + fileData.iteration = iter; + fileData.eta = eta; + fileData.time = t; + save(fileName, '-struct', 'fileData', '-v7.3'); + end + end + +end + diff --git a/src/environments/CtmrgEnvironment.m b/src/environments/CtmrgEnvironment.m new file mode 100644 index 0000000..18bb2da --- /dev/null +++ b/src/environments/CtmrgEnvironment.m @@ -0,0 +1,64 @@ +classdef CtmrgEnvironment + % Data structure for managing CTMRG environments. + % + % Properties + % ---------- + % corners : :class:`cell` of :class:`.Tensor` + % cell array of corner tensors. + % + % edges : :class:`cell` of :class:`.Tensor` + % cell array of edge tensors. + % + % Todo + % ---- + % Document. + + properties + corners + edges + end + + methods + + function obj = CtmrgEnvironment(varargin) + + if nargin == 0, return; end + + if iscell(varargin{1}) + obj.corners = varargin{1}; + obj.edges = varargin{2}; + end + + obj.corners = cellfun(@(x)x ./ norm(x), obj.corners,'UniformOutput',false); + obj.edges = cellfun(@(x)x ./ norm(x), obj.edges,'UniformOutput',false); + end + + function out = rot90(in) + %rotate unit cell + out = CtmrgEnvironment(flip(permute(in.corners,[1,3,2]),3),flip(permute(in.edges,[1,3,2]),3)); + %relabel corners and edges accordingly + in = out; + for i = 1:4 + out.corners(i,:,:) = in.corners(prev(i,4),:,:); + out.edges(i,:,:) = in.edges(prev(i,4),:,:); + end + + end + + function chi = bond_dimensions(obj) + chi = reshape(cellfun(@(x) dims(x,1),obj.corners(1,:,:)), size(obj.corners,2:3)); + end + + function h = height(obj) + % vertical period over which the peps is translation invariant. + h = height(obj.corners{1,:,:}); + end + + function w = width(obj) + % horizontal period over which the peps is translation invariant. + w = width(obj.corners{1,:,:}); + end + + end +end + diff --git a/src/mps/PepsSandwich.m b/src/mps/PepsSandwich.m index 0db3ebd..015670b 100644 --- a/src/mps/PepsSandwich.m +++ b/src/mps/PepsSandwich.m @@ -5,10 +5,10 @@ % Properties % ---------- % top : :class:`.PepsTensor` - % top-layer PEPS tensor, usually interpreted as the 'bra' in the overlap. + % top-layer PEPS tensor, usually interpreted as the 'ket' in the overlap. % % bot : :class:`.PepsTensor` - % bottom-layer PEPS tensor, usually interpreted as the 'ket' in the overlap. + % bottom-layer PEPS tensor, usually interpreted as the 'bra' in the overlap. % % Todo % ---- diff --git a/src/mps/PepsTensor.m b/src/mps/PepsTensor.m index ad21039..9e24d60 100644 --- a/src/mps/PepsTensor.m +++ b/src/mps/PepsTensor.m @@ -209,6 +209,11 @@ function t = randnc(varargin) t = PepsTensor.new(@randnc, varargin{:}); end + + function t = zeros(varargin) + t = PepsTensor.new(@zeros, varargin{:}); + end + end end diff --git a/src/mps/UniformPeps.m b/src/mps/UniformPeps.m new file mode 100644 index 0000000..dd881bb --- /dev/null +++ b/src/mps/UniformPeps.m @@ -0,0 +1,163 @@ +classdef UniformPeps + % Implementation of infinite translation invariant PEPS + % + % Properties + % ---------- + % A : :class:`cell` of :class:`.PepsTensor` + % cell array of PEPS tensors in 2D unit cell. + % + % Todo + % ---- + % Document. + + properties + A cell + end + + + %% Constructors + methods + function peps = UniformPeps(varargin) + % Usage + % ----- + % :code:`peps = UniformPeps(A)` + % + % Arguments + % --------- + % A : :class:`cell` of :class:`.PepsTensor` + % cell array of PEPS tensors in 2D unit cell. + % + % Returns + % ------- + % peps : :class:`.UniformPeps` + % infinite translation-invariant PEPS. + + if nargin == 0, return; end % default empty constructor + + if nargin == 1 + if isa(varargin{1}, 'UniformPeps') % copy constructor + peps.A = varargin{1}.A; + + elseif isa(varargin{1}, 'Tensor') + peps.A{1} = PepsTensor(varargin{1}); + + elseif isa(varargin{1}, 'PepsTensor') + peps.A{1} = varargin{1}; + + elseif iscell(varargin{1}) + assert(isa(varargin{1}{1},'PepsTensor')) + for i = height(varargin{1}):-1:1 + for j = width(varargin{1}):-1:1 + peps.A{i,j} = varargin{1}{i,j}; + end + end + else + error('Invalid constructor for UniformPeps.') + end + + else + error('Invalid constructor for UniformPeps.') + end + end + end + + %% Properties + methods + function h = height(peps) + % vertical period over which the peps is translation invariant. + h = height(peps.A); + end + + function w = width(peps) + % horizontal period over which the peps is translation invariant. + w = width(peps.A); + end + + function s = westvspace(peps, h, w) + % return the virtual space to the left of site (h,w). + if nargin == 1, h = 1:height(peps); w = 1:width(peps); end + s = cellfun(@westvspace, peps.A(h,w)); + end + + function s = southvspace(peps, h, w) + % return the virtual space to the bottom of site (h,w). + if nargin == 1, h = 1:height(peps); w = 1:width(peps); end + s = cellfun(@southvspace, peps.A(h,w)); + end + + function s = eastvspace(peps, h, w) + % return the virtual space to the right of site (h,w). + if nargin == 1, h = 1:height(peps); w = 1:width(peps); end + s = cellfun(@eastvspace, peps.A(h,w)); + end + + function s = northvspace(peps, h, w) + % return the virtual space to the top of site (h,w). + if nargin == 1, h = 1:height(peps); w = 1:width(peps); end + s = cellfun(@northvspace, peps.A(h,w)); + end + + function s = pspace(peps, h, w) + % return the physical space at site (h,w). + if nargin == 1, h = 1:height(peps); w = 1:width(peps); end + s = cellfun(@pspace, peps.A(h,w)); + end + + function peps = rot90(peps) + peps.A = cellfun(@(x)rot90(x),peps.A.','UniformOutput',false); + end + + function peps = rot270(peps) + peps.A = cellfun(@(x)rot270(x),peps.A.','UniformOutput',false); + end + + function type = underlyingType(peps) + type = underlyingType(peps.A{1,1}); + end + end + + + %% Methods + methods + %build horizontalTransferMatrix + + %build verticalTransferMatrix + + %CtmrgEnvironment + + function ctmrgenv = CtmrgEnvironment(peps_top, peps_bot, varargin) + + h = height(peps_top); + w = width(peps_top); + C = cell(4, h, w); + T = cell(4, h, w); + + if nargin == 2 + vspace = repmat(one(space(peps_top.A{1})), h, w); + elseif all(size(varargin{1})==[1,1]) + vspace = repmat(varargin{1}, h, w); + else + vspace = varargin{1}; + end + + for i = 1:h + for j = 1:w + C{1,i,j} = Tensor.randnc(vspace(i,j), vspace(i,j)); + C{2,i,j} = Tensor.randnc(vspace(i,prev(j,w)), vspace(i,prev(j,w))); + C{3,i,j} = Tensor.randnc(vspace(prev(i,h),prev(j,w)), vspace(prev(i,h),prev(j,w))); + C{4,i,j} = Tensor.randnc(vspace(prev(i,h),j), vspace(prev(i,h),j)); + + T{1,i,j} = Tensor.randnc([vspace(i,prev(j,w)), northvspace(peps_top,next(i,h),j)', northvspace(peps_bot,next(i,h),j)], vspace(i,j)); + T{2,i,j} = Tensor.randnc([vspace(prev(i,h),prev(j,w)), eastvspace(peps_top,i,prev(j,w))', eastvspace(peps_bot,i,prev(j,w))], vspace(i,prev(j,w))); + T{3,i,j} = Tensor.randnc([vspace(prev(i,h),j), southvspace(peps_top,prev(i,h),j)', southvspace(peps_bot,prev(i,h),j)], vspace(prev(i,h),prev(j,w))); + T{4,i,j} = Tensor.randnc([vspace(i,j), westvspace(peps_top,i,next(j,w))', westvspace(peps_bot,i,next(j,w))], vspace(prev(i,h),j)); + end + end + + ctmrgenv = CtmrgEnvironment(C, T); + + end + + end +end + diff --git a/src/tensors/Tensor.m b/src/tensors/Tensor.m index b9d1d98..f1a9e5f 100644 --- a/src/tensors/Tensor.m +++ b/src/tensors/Tensor.m @@ -1986,9 +1986,21 @@ end end if isfield(trunc, 'TruncSpace') - error('TBA'); + truncspace = trunc.TruncSpace; + [b, ind] = ismember(truncspace.dimensions.charges, dims.charges); + assert(all(b),"Truncation space contains charges that S does not.") + assert(all(dims.degeneracies(ind) >= truncspace.dimensions.degeneracies), "Truncation space has degeneracies larger than S.") + dims.degeneracies(ind) = truncspace.dimensions.degeneracies; + dims.degeneracies(setxor(ind,1:numel(dims.degeneracies))) = 0; + for i = 1:length(mblocks) + s = diag(Ss{i}); + eta = eta + sum(s(dims.degeneracies(i) + 1:end).^2 * qdim(dims.charges(i))); + Us{i} = Us{i}(:, 1:dims.degeneracies(i)); + Ss{i} = diag(s(1:dims.degeneracies(i))); + Vs{i} = Vs{i}(1:dims.degeneracies(i), :); + end end - + if ~doTrunc W1 = prod(t.codomain); W2 = prod(t.domain); @@ -2014,7 +2026,7 @@ U.var = fill_matrix_data(U.var, Us, dims.charges); S.var = fill_matrix_data(S.var, Ss, dims.charges); V.var = fill_matrix_data(V.var, Vs, dims.charges); - + if nargout <= 1 U = S; end diff --git a/src/utility/plot_schmidtspectrum.m b/src/utility/plot_schmidtspectrum.m new file mode 100644 index 0000000..c63c911 --- /dev/null +++ b/src/utility/plot_schmidtspectrum.m @@ -0,0 +1,59 @@ +function [ax] = plot_schmidtspectrum(S, ax, kwargs) + +arguments + S + ax = [] + kwargs.SymmetrySort = true + kwargs.ExpandQdim = false + kwargs.Marker = '.' +end + + if isempty(ax) + figure; + ax = gca; + end + + if kwargs.Marker == '.' + kwargs.MarkerSize = 10; + else + kwargs.MarkerSize = 6; + end + + %hold(ax, 'off'); + [svals, charges] = matrixblocks(S); + if kwargs.ExpandQdim + for i = 1:length(svals) + svals{i} = reshape(repmat(svals{i}, 1, qdim(charges(i))), [], 1); + end + end + ctr = 0; + labels = arrayfun(@string, charges, 'UniformOutput', false); + lengths = cellfun(@length, svals); + ticks = cumsum(lengths); + if kwargs.SymmetrySort + for i = 1:length(svals) + semilogy(ax, ctr+(1:lengths(i)), svals{i}, 'Marker', kwargs.Marker, 'MarkerSize', kwargs.MarkerSize, 'Color', colors(i)); + if i == 1, hold(ax, 'on'); end + ctr = ctr + lengths(i); + end + set(ax, 'Xtick', ticks, 'fontsize', 10, ... + 'XtickLabelRotation', 60, 'Xgrid', 'on'); + else + [~, p] = sort(vertcat(svals{:}), 'descend'); + p = invperm(p); + for i = 1:length(svals) + semilogy(ax, p(ctr+(1:lengths(i))), svals{i}, 'Marker', kwargs.Marker, 'MarkerSize', kwargs.MarkerSize, 'Color', colors(i)); + if i == 1, hold(ax, 'on'); end + ctr = ctr + lengths(i); + end + + end + legend(ax, arrayfun(@(x) line([0,1],[0,0],'Color', colors(x), 'Marker', '.', 'MarkerSize', 10, 'linestyle', 'none'),1:numel(labels)), labels, 'Location', 'Best') + set(ax, 'TickLabelInterpreter', 'latex'); + xlim(ax, [1 - 1e-8 ticks(end) + 1e-8]); + hold off + + linkaxes(ax, 'y'); + +end + diff --git a/test/TestCtmrg.m b/test/TestCtmrg.m new file mode 100644 index 0000000..40b312f --- /dev/null +++ b/test/TestCtmrg.m @@ -0,0 +1,76 @@ +classdef TestCtmrg < matlab.unittest.TestCase + % Unit tests for infinite matrix product operators. + + properties (TestParameter) + pspace = struct('fermion', GradedSpace.new(fZ2(0, 1), [1 1], false)) + vspace = struct('fermion', GradedSpace.new(fZ2(0, 1), [1 1], false)) + chispace = struct('fermion', GradedSpace.new(fZ2(0, 1), [10 10], false)) + alg = {Ctmrg('tol', 1e-6)} + %other symms TBA + %larger unit cells TBA + end + + methods (Test, ParameterCombination='sequential') + function testEnvironments(tc, pspace, vspace, chispace) + peps = UniformPeps(PepsTensor.randnc(pspace, vspace, vspace)); + envs = CtmrgEnvironment(peps, peps, chispace); + %check matching spaces TBA + end + + function testConvergence(tc, pspace, vspace, chispace, alg) + maxloops = 10; + i = 0; + err = 1; + while err > alg.tol && i < maxloops + i = i+1; + peps = UniformPeps(PepsTensor.randnc(pspace, vspace, vspace)); + envs = CtmrgEnvironment(peps, peps, chispace); + [envs, norm, err] = fixedpoint(alg, peps, peps, envs); + end + assert(i < maxloops, "Convergence seems to fail for simple random PEPS.") + %check left move TBA + %check norm TBA + end + + function testPWave(tc) + %initialize + pspace = GradedSpace.new(fZ2(0, 1), [1 1], false); + vspace = pspace; + chispace = GradedSpace.new(fZ2(0, 1), [10 10], false); + alg = Ctmrg('tol', 1e-6, 'miniter', 1, 'maxiter', 200); + %fill PEPS as PWave from Gaussian + mblocks{1} = [-0.4084 - 0.5139i -0.1033 + 0.1492i -0.7413 - 1.1574i 0.5924 + 0.6552i ... + -0.0660 - 0.1002i -0.1678 - 0.3237i -0.3739 + 0.9139i 0.6497 + 0.8648i]; + mblocks{2} = [ 0.5454 + 0.6080i -0.3003 - 0.2482i 0.3430 - 0.1734i -0.2851 - 0.2658i ... + -0.4272 + 0.5267i 0.4567 + 0.5106i -0.6708 - 0.1306i 0.2184 - 0.1094i]; + A = PepsTensor.zeros(pspace, vspace, vspace); + peps = UniformPeps(PepsTensor(A.var.fill_matrix(mblocks))); + A = peps.A{1}.var; + %run ctmrg + maxloops = 10; + i = 0; + err = 1; + while err > alg.tol && i < maxloops + i = i+1; + envs = CtmrgEnvironment(peps, peps, chispace); + [envs, norm, err] = fixedpoint(alg, peps, peps, envs); + end + assert(i < maxloops, "Convergence seems to fail for PWave example.") + %local density comparison to Gaussian + O = Tensor(pspace, pspace); + O.var.var{2} = 1; + GL = contract(envs.corners{1,1,1}, [1 -4], envs.edges{4,1,1}, [2 -2 -3 1], envs.corners{4,1,1}, [-1 2]); + GR = contract(envs.corners{2,1,1}, [-1 1], envs.edges{2,1,1}, [1 -2 -3 2], envs.corners{3,1,1}, [2 -4]); + n = contract(GL, 1:4, GR, [4,2,3,1]); + GL = GL/sqrt(n); + GR = GR/sqrt(n); + l = contract(envs.edges{3,1,1}, [-1 5 2 1], GL, [1 4 3 8], envs.edges{1,1,1}, [8 9 10 -4], ... + A, [7 4 5 -2 9], conj(A), [6 3 2 -3 10], O, [6 7]); + lw = contract(envs.edges{3,1,1}, [-1 5 2 1], GL, [1 4 3 8], envs.edges{1,1,1}, [8 9 10 -4], ... + A, [6 4 5 -2 9], conj(A), [6 3 2 -3 10]); + rho = contract(l, [2 3 4 1], GR, [1 3 4 2])/contract(lw, [2 3 4 1], GR, [1 3 4 2]); + assert(abs(rho - 0.178787157813086) < sqrt(alg.tol), "CTMRG yields wrong local density for PWave example.") + end + end +end +