From cdf67bca67c7a19e5e988c5914452073ab572851 Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Sun, 19 May 2024 11:56:44 +0800 Subject: [PATCH] heisenberg mpo --- Project.toml | 1 + pycode/doVUMPS.py | 228 ++++++++++++++++++++++++++++++++++++++++++++++ pycode/ncon.py | 197 +++++++++++++++++++++++++++++++++++++++ pycode/vumps.py | 97 ++++++++++++++++++++ src/SimpleTDVP.jl | 4 + src/dmrg.jl | 52 +++++++++++ src/mpo.jl | 12 +++ src/mps.jl | 26 ++++++ test/mpo.jl | 10 +- test/mps.jl | 10 ++ 10 files changed, 636 insertions(+), 1 deletion(-) create mode 100644 pycode/doVUMPS.py create mode 100644 pycode/ncon.py create mode 100644 pycode/vumps.py diff --git a/Project.toml b/Project.toml index b8d14f2..7b10050 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" LuxorGraphPlot = "1f49bdf2-22a7-4bc4-978b-948dc219fbbc" OMEinsum = "ebe7aa44-baf0-506c-a96f-8464559b3922" +Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c" [compat] KrylovKit = "0.7" diff --git a/pycode/doVUMPS.py b/pycode/doVUMPS.py new file mode 100644 index 0000000..2914f56 --- /dev/null +++ b/pycode/doVUMPS.py @@ -0,0 +1,228 @@ +""" +doVUMPS.py +--------------------------------------------------------------------- +Module for optimizing an infinite MPS using the VUMPS algorithm. + +by Glen Evenbly (c) for www.tensors.net, (v1.0) - last modified 07/2020 +""" + +import numpy as np +from numpy import linalg as LA +from ncon import ncon +from scipy.sparse.linalg import LinearOperator, eigsh, gmres +from scipy.linalg import polar +from typing import List, Optional + + +def doVUMPS(AL: np.ndarray, + C: np.ndarray, + AR: np.ndarray, + h0: np.ndarray, + HL: Optional[np.ndarray] = None, + HR: Optional[np.ndarray] = None, + m: Optional[int] = None, + update_mode: Optional[str] = 'polar', + num_iter: Optional[int] = 1000, + ev_tol: Optional[float] = 1e-12, + en_exact: Optional[float] = 0.0): + """ + Implementation the VUMPS algorithm to optimize an infinite MPS with a 1-site + unit cell. + Args: + AL: the rank-3 MPS tensor in the left orthogonal gauge. + C: vector of the singular weights of the MPS. + AR: the rank-3 MPS tensor in the right orthogonal gauge. + h0: two-site Hamiltonian coupling. + HL: the left block Hamiltonian (if previously computed). + HR: the right block Hamiltonian (if previously computed). + m: the MPS bond dimension. If the input MPS tensors are of smaller + dimension, then they will be exapnded to match `m`. + update_mode: set updates via 'polar' decomposition (more stable) of the + 'svd' (less stable). + num_iter: number of iterations to perform. + ev_tol: tolerance for convergence of eighs. + en_exact: the exact energy (if known). + Returns: + np.ndarray: MPS tensor AL + np.ndarray: MPS singular weights C + np.ndarray: MPS tensor AR + np.ndarray: block Hamiltonian HL + np.ndarray: block Hamiltonian HR + """ + + # Initialise tensors + d = h0.shape[0] + h = h0.copy() + if HL is None: + HL = np.zeros([m, m]) + if HR is None: + HR = np.zeros([m, m]) + + if m is None: + m = AL.shape[2] + else: + if m > AL.shape[2]: + # Expand tensors to new dimension + AL = Orthogonalize(TensorExpand(AL, [m, d, m]), 2) + C = TensorExpand(C, [m]) + AR = Orthogonalize(TensorExpand(AR, [m, d, m]), 2) + HL = TensorExpand(HL, [m, m]) + HR = TensorExpand(HR, [m, m]) + AC = ncon([AL, np.diag(C)], [[-1, -2, 1], [1, -3]]) + + # Begin variational update iterations + for p in range(num_iter): + """ Evaluate the energy """ + tensors = [AL, AL, h0, AL.conj(), AL.conj()] + connects = [[7, 1, 6], [6, 2, -1], [3, 4, 1, 2], [7, 3, 5], [5, 4, -2]] + con_order = [7, 1, 3, 6, 2, 5, 4] + hL0 = ncon(tensors, connects, con_order) + energyL = np.trace(np.diag(C**2) @ hL0) + + tensors = [AR, AR, h0, AR.conj(), AR.conj()] + connects = [[6, 4, 1], [1, 3, -1], [5, 7, 3, 4], [6, 7, 2], [2, 5, -2]] + con_order = [6, 4, 7, 1, 3, 2, 5] + hR0 = ncon(tensors, connects, con_order) + energyR = np.trace(np.diag(C**2) @ hR0) + + energy = 0.25 * (energyL + energyR) + en_error = energy - en_exact + print('iteration: %d of %d, dim: %d, energy: %4.4f, ' + 'en-error: %2.2e' % (p, num_iter, m, energy, en_error)) + + """ Contract the MPS from the left """ + # Shift left edge Hamiltonian terms to set energy to zero + tensors = [AL, AL, h, AL.conj(), AL.conj()] + connects = [[7, 1, 6], [6, 2, -1], [3, 4, 1, 2], [7, 3, 5], [5, 4, -2]] + con_order = [7, 1, 3, 6, 2, 5, 4] + hL = ncon(tensors, connects, con_order) + + en_density = np.dot(np.diag(hL), C**2) + hL -= en_density * np.eye(m) + h -= en_density * np.eye(d**2).reshape(d, d, d, d) + HL -= np.dot(np.diag(HL), C**2) * np.eye(m) + + # Solve left edge Hamiltonian eigenvector + def LeftEdge(HL): + m = AL.shape[0] + d = AL.shape[1] + HL_T = AL.reshape(m * d, m).T @ (HL.reshape(m, m) @ AL.reshape(m, d * m) + ).reshape(m * d, m) + return HL.flatten() - HL_T.flatten() + LeftOp = LinearOperator((m**2, m**2), matvec=LeftEdge, dtype=np.float64) + HL_temp, is_conv = gmres(LeftOp, hL.flatten(), x0=HL.flatten(), tol=1e-10, + restart=None, maxiter=5, atol=None) + HL_temp = HL_temp.reshape(m, m) + HL = 0.5 * (HL_temp + HL_temp.T) + + """ Contract the MPS from the right """ + # Shift right edge Hamiltonian terms to set energy to zero + tensors = [AR, AR, h, AR.conj(), AR.conj()] + connects = [[6, 4, 1], [1, 3, -1], [5, 7, 3, 4], [6, 7, 2], [2, 5, -2]] + con_order = [6, 4, 7, 1, 3, 2, 5] + hR = ncon(tensors, connects, con_order) + en_density = np.dot(np.diag(hR), C**2) + + hR -= en_density * np.eye(m) + h -= np.dot(np.diag(hR), C**2) * np.eye(d**2).reshape(d, d, d, d) + HR -= np.dot(np.diag(HR), C**2) * np.eye(m) + + # Solve right edge Hamiltonian eigenvector + def RightEdge(HR): + m = AR.shape[0] + d = AR.shape[1] + HR_T = AR.reshape(m * d, m).T @ (HR.reshape(m, m) @ AR.reshape(m, d * m) + ).reshape(m * d, m) + return HR.flatten() - HR_T.flatten() + RightOp = LinearOperator((m**2, m**2), matvec=RightEdge, dtype=np.float64) + HR_temp, is_conv = gmres(RightOp, hR.flatten(), x0=HR.flatten(), tol=1e-10, + restart=None, maxiter=5, atol=None) + HR_temp = HR_temp.reshape(m, m) + HR = 0.5 * (HR_temp + HR_temp.T) + + """ Update the MPS singular values """ + # define function for applying the block Hamiltonian + def MidWeights(C_mat): + m = AL.shape[2] + C_mat = C_mat.reshape(m, m) + tensors = [AL, h, AL.conj(), AR, AR.conj(), C_mat] + connects = [[3, 1, 8], [2, 7, 1, 5], [3, 2, -1], [6, 5, 4], [6, 7, -2], + [8, 4]] + con_order = [4, 8, 1, 5, 6, 7, 3, 2] + return (ncon(tensors, connects, con_order) + (HL @ C_mat) + + (C_mat @ HR)).flatten() + + WeightOp = LinearOperator((m**2, m**2), matvec=MidWeights, dtype=np.float64) + C_temp = eigsh(WeightOp, k=1, which='SA', v0=np.diag(C).flatten(), + ncv=None, maxiter=None, tol=ev_tol)[1] + + # Change to diagonal gauge + ut, C, vt = LA.svd(C_temp.reshape(m, m)) + AL = ncon([ut.T, AL, ut], [[-1, 1], [1, -2, 2], [2, -3]]) + AR = ncon([vt, AR, vt.T], [[-1, 1], [1, -2, 2], [2, -3]]) + HL = ut.T @ HL @ ut + HR = vt @ HR @ vt.T + + # Contract left/middle Hamiltonian term + tensors = [AL, h0, AL.conj()] + connects = [[3, 1, -1], [2, -4, 1, -2], [3, 2, -3]] + con_order = [2, 3, 1] + hL_mid = ncon(tensors, connects, con_order).reshape(d * m, d * m) + + # Contract right/middle Hamiltonian term + tensors = [AR, h0, AR.conj()] + connects = [[2, 1, -2], [-3, 3, -1, 1], [2, 3, -4]] + con_order = [2, 1, 3] + hR_mid = ncon(tensors, connects, con_order).reshape(d * m, d * m) + + """ Update the MPS tensors """ + # define function for applying the block Hamiltonian + def MidTensor(AC): + m = AL.shape[2] + d = AL.shape[1] + return ((HL @ AC.reshape(m, d * m)).flatten() + + (hL_mid.reshape(m * d, m * d) @ AC.reshape(m * d, m)).flatten() + + (AC.reshape(m * d, m) @ HR).flatten() + + (AC.reshape(m, d * m) @ hR_mid.reshape(d * m, d * m)).flatten()) + + TensorOp = LinearOperator((d * m**2, d * m**2), + matvec=MidTensor, dtype=np.float64) + AC = (eigsh(TensorOp, k=1, which='SA', v0=AC.flatten(), + ncv=None, maxiter=None, tol=ev_tol)[1]).reshape(m, d, m) + + if update_mode == 'polar': + AL = (polar(AC.reshape(m * d, m))[0]).reshape(m, d, m) + AR = (polar(AC.reshape(m, d * m), side='left')[0] + ).reshape(m, d, m).transpose(2, 1, 0) + elif update_mode == 'svd': + ut, _, vt = LA.svd(AC.reshape(m * d, m) @ np.diag(C), full_matrices=False) + AL = (ut @ vt).reshape(m, d, m) + ut, _, vt = LA.svd(np.diag(C) @ AC.reshape(m, d * m), full_matrices=False) + AR = (ut @ vt).reshape(m, d, m).transpose(2, 1, 0) + + return AL, C, AR, HL, HR + + +def TensorExpand(A: np.ndarray, chivec: List): + """ expand tensor dimension by padding with zeros """ + + if [*A.shape] == chivec: + return A + else: + for k in range(len(chivec)): + if A.shape[k] != chivec[k]: + indloc = list(range(-1, -len(chivec) - 1, -1)) + indloc[k] = 1 + A = ncon([A, np.eye(A.shape[k], chivec[k])], [indloc, [1, -k - 1]]) + + return A + + +def Orthogonalize(A: np.ndarray, pivot: int): + """ orthogonalize an array with respect to a pivot """ + + A_sh = A.shape + ut, st, vht = LA.svd( + A.reshape(np.prod(A_sh[:pivot]), np.prod(A_sh[pivot:])), + full_matrices=False) + return (ut @ vht).reshape(A_sh) diff --git a/pycode/ncon.py b/pycode/ncon.py new file mode 100644 index 0000000..b35d9cd --- /dev/null +++ b/pycode/ncon.py @@ -0,0 +1,197 @@ +# ncon network contractor +# by Glen Evenbly (c) for www.tensors.net, (v1.2) - last modified 6/2020 + +# ncon.py +import numpy as np +from typing import List, Union, Tuple, Optional + + +def ncon(tensors: List[np.ndarray], + connects: List[Union[List[int], Tuple[int]]], + con_order: Optional[Union[List[int], str]] = None, + check_network: Optional[bool] = True, + which_env: Optional[int] = 0): + """ + Network CONtractor: contracts a tensor network of N tensors via a sequence + of (N-1) tensordot operations. More detailed instructions and examples can + be found at: https://arxiv.org/abs/1402.0939. + Args: + tensors: list of the tensors in the network. + connects: length-N list of lists (or tuples) specifying the network + connections. The jth entry of the ith list in connects labels the edge + connected to the jth index of the ith tensor. Labels should be positive + integers for internal indices and negative integers for free indices. + con_order: optional argument to specify the order for contracting the + positive indices. Defaults to ascending order if omitted. Can also be + set at "greedy" or "full" to call a solver to automatically determine + the order. + check_network: if true then the input network is checked for consistency; + this can catch many common user mistakes for defining networks. + which_env: if provided, ncon will produce the environment of the requested + tensor (i.e. the network given by removing the specified tensor from + the original network). Only valid for networks with no open indices. + Returns: + Union[np.ndarray,float]: the result of the network contraction; an + np.ndarray if the network contained open indices, otherwise a scalar. + """ + num_tensors = len(tensors) + tensor_list = [tensors[ele] for ele in range(num_tensors)] + connect_list = [np.array(connects[ele]) for ele in range(num_tensors)] + + # generate contraction order if necessary + flat_connect = np.concatenate(connect_list) + if con_order is None: + con_order = np.unique(flat_connect[flat_connect > 0]) + else: + con_order = np.array(con_order) + + # check inputs if enabled + if check_network: + dims_list = [list(tensor.shape) for tensor in tensor_list] + check_inputs(connect_list, flat_connect, dims_list, con_order) + + # do all partial traces + for ele in range(len(tensor_list)): + num_cont = len(connect_list[ele]) - len(np.unique(connect_list[ele])) + if num_cont > 0: + tensor_list[ele], connect_list[ele], cont_ind = partial_trace( + tensor_list[ele], connect_list[ele]) + con_order = np.delete( + con_order, + np.intersect1d(con_order, cont_ind, return_indices=True)[1]) + + # do all binary contractions + while len(con_order) > 0: + # identify tensors to be contracted + cont_ind = con_order[0] + locs = [ + ele for ele in range(len(connect_list)) + if sum(connect_list[ele] == cont_ind) > 0 + ] + + # do binary contraction + cont_many, A_cont, B_cont = np.intersect1d( + connect_list[locs[0]], + connect_list[locs[1]], + assume_unique=True, + return_indices=True) + if np.size(tensor_list[locs[0]]) < np.size(tensor_list[locs[1]]): + ind_order = np.argsort(A_cont) + else: + ind_order = np.argsort(B_cont) + + tensor_list.append( + np.tensordot( + tensor_list[locs[0]], + tensor_list[locs[1]], + axes=(A_cont[ind_order], B_cont[ind_order]))) + connect_list.append( + np.append( + np.delete(connect_list[locs[0]], A_cont), + np.delete(connect_list[locs[1]], B_cont))) + + # remove contracted tensors from list and update con_order + del tensor_list[locs[1]] + del tensor_list[locs[0]] + del connect_list[locs[1]] + del connect_list[locs[0]] + con_order = np.delete( + con_order, + np.intersect1d(con_order, cont_many, return_indices=True)[1]) + + # do all outer products + while len(tensor_list) > 1: + s1 = tensor_list[-2].shape + s2 = tensor_list[-1].shape + tensor_list[-2] = np.outer(tensor_list[-2].reshape(np.prod(s1)), + tensor_list[-1].reshape(np.prod(s2))).reshape( + np.append(s1, s2)) + connect_list[-2] = np.append(connect_list[-2], connect_list[-1]) + del tensor_list[-1] + del connect_list[-1] + + # do final permutation + if len(connect_list[0]) > 0: + return np.transpose(tensor_list[0], np.argsort(-connect_list[0])) + else: + return tensor_list[0].item() + + +def partial_trace(A, A_label): + """ Partial trace on tensor A over repeated labels in A_label """ + + num_cont = len(A_label) - len(np.unique(A_label)) + if num_cont > 0: + dup_list = [] + for ele in np.unique(A_label): + if sum(A_label == ele) > 1: + dup_list.append([np.where(A_label == ele)[0]]) + + cont_ind = np.array(dup_list).reshape(2 * num_cont, order='F') + free_ind = np.delete(np.arange(len(A_label)), cont_ind) + + cont_dim = np.prod(np.array(A.shape)[cont_ind[:num_cont]]) + free_dim = np.array(A.shape)[free_ind] + + B_label = np.delete(A_label, cont_ind) + cont_label = np.unique(A_label[cont_ind]) + B = np.zeros(np.prod(free_dim)) + A = A.transpose(np.append(free_ind, cont_ind)).reshape( + np.prod(free_dim), cont_dim, cont_dim) + for ip in range(cont_dim): + B = B + A[:, ip, ip] + + return B.reshape(free_dim), B_label, cont_label + + else: + return A, A_label, [] + + +def check_inputs(connect_list, flat_connect, dims_list, con_order): + """ Check consistancy of NCON inputs""" + + pos_ind = flat_connect[flat_connect > 0] + neg_ind = flat_connect[flat_connect < 0] + + # check that lengths of lists match + if len(dims_list) != len(connect_list): + raise ValueError( + ('mismatch between %i tensors given but %i index sublists given') % + (len(dims_list), len(connect_list))) + + # check that tensors have the right number of indices + for ele in range(len(dims_list)): + if len(dims_list[ele]) != len(connect_list[ele]): + raise ValueError(( + 'number of indices does not match number of labels on tensor %i: ' + '%i-indices versus %i-labels') + % (ele, len(dims_list[ele]), len(connect_list[ele]))) + + # check that contraction order is valid + if not np.array_equal(np.sort(con_order), np.unique(pos_ind)): + raise ValueError(('NCON error: invalid contraction order')) + + # check that negative indices are valid + for ind in np.arange(-1, -len(neg_ind) - 1, -1): + if sum(neg_ind == ind) == 0: + raise ValueError(('NCON error: no index labelled %i') % (ind)) + elif sum(neg_ind == ind) > 1: + raise ValueError(('NCON error: more than one index labelled %i') % (ind)) + + # check that positive indices are valid and contracted tensor dimensions match + flat_dims = np.array([item for sublist in dims_list for item in sublist]) + for ind in np.unique(pos_ind): + if sum(pos_ind == ind) == 1: + raise ValueError(('NCON error: only one index labelled %i') % (ind)) + elif sum(pos_ind == ind) > 2: + raise ValueError( + ('NCON error: more than two indices labelled %i') % (ind)) + + cont_dims = flat_dims[flat_connect == ind] + if cont_dims[0] != cont_dims[1]: + raise ValueError( + ('NCON error: tensor dimension mismatch on index labelled %i: ' + 'dim-%i versus dim-%i') % (ind, cont_dims[0], cont_dims[1])) + + return True + diff --git a/pycode/vumps.py b/pycode/vumps.py new file mode 100644 index 0000000..7e6a86e --- /dev/null +++ b/pycode/vumps.py @@ -0,0 +1,97 @@ +""" +mainVUMPS.py +--------------------------------------------------------------------- +Script file for initializing the Hamiltonian and MPS tensors before passing to +the VUMPS routine to optimize for the ground state. + +by Glen Evenbly (c) for www.tensors.net, (v1.0) - last modified 07/2020 +""" + +import numpy as np +from numpy import linalg as LA +from scipy.sparse.linalg import LinearOperator, eigs +from ncon import ncon +from doVUMPS import doVUMPS + + +""" set algorithm options """ +m = 32 # MPS bond dimension +update_mode = 'polar' # set update decomposition ('polar' or 'svd') +ev_tol = 1e-5 # tolerance for convergence of eighs +model = 'heisenberg' # set model ('heisenberg', 'XX' or 'ising') +num_iter = 20 # number of iterations to perform + + +""" define the Hamiltonian """ +d = 4 +sX = np.array([[0, 1.0], [1.0, 0]]) +sY = np.array([[0, -1.0 * 1j], [1.0 * 1j, 0]]) +sZ = np.array([[1.0, 0], [0, -1.0]]) +sI = np.array([[1.0, 0], [0, 1.0]]) +if model == 'heisenberg': + h_temp = (np.real(np.kron(sX, sX) + np.kron(sY, sY) + np.kron(sZ, sZ))) + h0 = (0.5 * np.kron(h_temp, np.eye(4)) + 0.5 * np.kron(np.eye(4), h_temp) + + np.kron(np.eye(2), np.kron(h_temp, np.eye(2)))).reshape(d, d, d, d) + en_exact = (1 - 4 * np.log(2)) +elif model == 'XX': + h_temp = (np.real(np.kron(sX, sX) + np.kron(sY, sY))) + h0 = (0.5 * np.kron(h_temp, np.eye(4)) + 0.5 * np.kron(np.eye(4), h_temp) + + np.kron(np.eye(2), np.kron(h_temp, np.eye(2)))).reshape(d, d, d, d) + en_exact = -4 / np.pi +elif model == 'ising': + h_temp = (-np.real(np.kron(sX, sX) + 0.5 * np.kron(sZ, sI) + + 0.5 * np.kron(sI, sZ))) + h0 = (0.5 * np.kron(h_temp, np.eye(4)) + 0.5 * np.kron(np.eye(4), h_temp) + + np.kron(np.eye(2), np.kron(h_temp, np.eye(2)))).reshape(d, d, d, d) + en_exact = -4 / np.pi + +""" initialize the MPS tensors """ +C = np.random.rand(m) +C = C / LA.norm(C) +AL = (LA.svd(np.random.rand(m * d, m), full_matrices=False)[0]).reshape(m, d, m) +AR = (LA.svd(np.random.rand(m * d, m), full_matrices=False)[0]).reshape(m, d, m) +HL = np.zeros([m, m]) +HR = np.zeros([m, m]) + + +""" run the optimization algorithm """ +AL, C, AR, HL, HR = doVUMPS(AL, C, AR, h0, HL=HL, HR=HR, m=m, + update_mode=update_mode, num_iter=num_iter, + ev_tol=ev_tol, en_exact=en_exact) + + +""" increase bond dim, reduce tolerance, then continue updates """ +m = 64 +ev_tol = 1e-7 +AL, C, AR, HL, HR = doVUMPS(AL, C, AR, h0, HL=HL, HR=HR, m=m, + update_mode=update_mode, num_iter=num_iter, + ev_tol=ev_tol, en_exact=en_exact) + + +""" compute the local 2-site reduced density matrix """ +# define linear operator for contraction from the right +def RightDensity(rhoR): + m = AL.shape[2] + tensors = [AL, AL.conj(), rhoR.reshape(m, m)] + connects = [[-2, 1, 2], [-1, 1, 3], [3, 2]] + con_order = [2, 1, 3] + return ncon(tensors, connects, con_order).flatten() + + +RightDensityOp = LinearOperator((m**2, m**2), matvec=RightDensity, + dtype=np.float64) +# solve eigenvalue problem for r.h.s density matrix +rhoR_temp = eigs(RightDensityOp, k=1, which='LM', v0=np.diag(C**2).flatten(), + ncv=None, maxiter=None, tol=0)[1] +rhoR_temp = np.real(rhoR_temp.reshape(m, m)) +rhoR_temp = rhoR_temp + rhoR_temp.T.conj() +rhoR = rhoR_temp / np.trace(rhoR_temp) + +# contract for the 2-site local reduced density matrix +tensors = [AL, AL, AL.conj(), AL.conj(), rhoR] +connects = [[3, -1, 4], [4, -2, 1], [3, -3, 5], [5, -4, 2], [1, 2]] +rho_two = ncon(tensors, connects) + +en_final = ncon([h0, rho_two], [[1, 2, 3, 4], [1, 2, 3, 4]]) / 2 +en_error = en_final - en_exact +print('final energy: %4.4f, en-error: %2.2e' % (en_final, en_error)) diff --git a/src/SimpleTDVP.jl b/src/SimpleTDVP.jl index 511f7f5..665abcd 100644 --- a/src/SimpleTDVP.jl +++ b/src/SimpleTDVP.jl @@ -2,11 +2,15 @@ module SimpleTDVP using LinearAlgebra, OMEinsum, LinearOperators, KrylovKit using LuxorGraphPlot: Node, Connection, Luxor +using Yao +import Yao: mat import LuxorGraphPlot export MPS, nsite, nflavor, vec2mps, IndexStore, newindex!, code_mps2vec, rand_mps, vec +export product_mps, ghz_mps export left_move!, right_move!, canonical_move!, is_canonicalized, canonical_center, to_left_canonical!, to_right_canonical! export mat, MPO, code_mpo2mat, mat2mpo, rand_mpo +export heisenberg_mpo export dot, sandwich, compress!, num_of_elements export dmrg!, dmrg export TensorLayout diff --git a/src/dmrg.jl b/src/dmrg.jl index 251f5c1..9ac8f28 100644 --- a/src/dmrg.jl +++ b/src/dmrg.jl @@ -59,6 +59,58 @@ function dmrg!(mps::MPS{T}, mpo::MPO{T}; nsweeps, atol=1e-10, Dmax=typemax(Int)) return DMRGResult(eng, mps, truncation_error, convergence_error) end +function idmrg!(mps::MPS{T}, mpo::MPO{T}; nsweeps, atol=1e-10, Dmax=typemax(Int)) where T + @assert nsite(mps) == nsite(mpo) "The number of sites in MPS and MPO are different" + @assert nflavor(mps) == nflavor(mpo) "The number of flavors in MPS and MPO are different" + to_right_canonical!(mps) + left_env = empty(mps.tensors) + right_env = [fill!(similar(mps.tensors[end], 1, 1, 1), one(T))] + for i=nsite(mps):-1:3 # update the right environment 3:nsite(mps) + push!(right_env, update_right_env(right_env[end], mps.tensors[i], mpo.tensors[i])) + end + convergence_error = eng = eng_prev = real(T)(Inf) + truncation_error = zero(real(T)) + for k = 1:nsweeps + for i = 1:nsite(mps)-1 # NOTE: site `1` is unneseccary to update, but is convenient to have + # update site i, i+1 + if i == 1 + L = fill!(similar(mps.tensors[1], 1, 1, 1), one(T)) + else + L = update_left_env(left_env[end], mps.tensors[i-1], mpo.tensors[i-1]) + end + push!(left_env, L) + r = DMRGRuntime(mps.tensors[i], mps.tensors[i+1], mpo.tensors[i], mpo.tensors[i+1], L, pop!(right_env)) + eng, A, S, B, err = update_sites(r; Dmax, atol) + B .*= reshape(S, :, 1, 1) + truncation_error += err + @info "Sweep = $k (right moving), site = ($i, $(i+1)), energy = $eng, bond size = $(length(S)), error = $err" + mps.tensors[i], mps.tensors[i+1] = A, B + right_move!(mps, Dmax=Dmax, atol=atol) + end + #push!(right_env, fill!(similar(mps.tensors[end], 1, 1, 1), one(T))) + for i = nsite(mps):-1:2 # NOTE: site `nsite(mps)` is unneseccary to update, but is convenient to have + # update site i-1, i + if i == nsite(mps) + R = fill!(similar(mps.tensors[end], 1, 1, 1), one(T)) + else + R = update_right_env(right_env[end], mps.tensors[i+1], mpo.tensors[i+1]) + end + push!(right_env, R) + r = DMRGRuntime(mps.tensors[i-1], mps.tensors[i], mpo.tensors[i-1], mpo.tensors[i], pop!(left_env), R) + eng, A, S, B, err = update_sites(r; Dmax, atol) + A .*= reshape(S, 1, 1, :) + truncation_error += err + @info "Sweep = $k (left moving), site = ($(i-1), $i), energy = $eng, bond size = $(length(S)), error = $err" + mps.tensors[i-1], mps.tensors[i] = A, B + left_move!(mps, Dmax=Dmax, atol=atol) + end + convergence_error = eng_prev - eng + eng_prev = eng + end + return DMRGResult(eng, mps, truncation_error, convergence_error) +end + + struct DMRGRuntime{T, ST<:AbstractArray{T, 3}, OT<:AbstractArray{T, 4}, TL<:AbstractArray{T, 3}, TR<:AbstractArray{T, 3}} A::ST B::ST diff --git a/src/mpo.jl b/src/mpo.jl index cb19d51..bd3e8ba 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -109,4 +109,16 @@ end function sandwich(bra::MPS, op::MPO, ket::MPS) code = code_sandwich(bra, op, ket) return code(conj.(bra.tensors)..., op.tensors..., ket.tensors...)[] +end + +function heisenberg_mpo(::Type{T}, n::Int) where T + @assert n > 1 + tensor1 = zeros(T, 1, 2, 2, 5) + tensor2 = zeros(T, 5, 2, 2, 5) + tensor3 = zeros(T, 5, 2, 2, 1) + tensor1[1, :, :, 1] = tensor2[1, :, :, 1] = tensor2[5, :, :, 5] = tensor3[5, :, :, 1] = Matrix{T}(I2) + tensor1[1, :, :, 2] = tensor2[2, :, :, 5] = tensor2[1, :, :, 2] = tensor3[2, :, :, 1] = Matrix{T}(X) + tensor1[1, :, :, 3] = tensor2[3, :, :, 5] = tensor2[1, :, :, 3] = tensor3[3, :, :, 1] = Matrix{T}(Y) + tensor1[1, :, :, 4] = tensor2[4, :, :, 5] = tensor2[1, :, :, 4] = tensor3[4, :, :, 1] = Matrix{T}(Z) + MPO([tensor1, fill(tensor2, n-2)..., tensor3]) end \ No newline at end of file diff --git a/src/mps.jl b/src/mps.jl index 84bd455..6bef2b3 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -101,4 +101,30 @@ end function LinearAlgebra.dot(mps1::MPS, mps2::MPS) code = code_dot(mps1, mps2) return code(conj.(mps1.tensors)..., mps2.tensors...)[] +end + +function product_mps(states::AbstractVector...) + return MPS([reshape(state, 1, size(state, 1), 1) for state in states]) +end + +function ghz_mps(::Type{T}, n::Int) where T + @assert n >= 2 + state1 = zeros(T, 1, 2, 2) + state1[1, 1, 1] = state1[1, 2, 2] = 1/sqrt(2) + state2 = zeros(T, 2, 2, 2) + state2[1, 1, 1] = state2[2, 2, 2] = 1 + state3 = zeros(T, 2, 2, 1) + state3[1, 1, 1] = state3[2, 2, 1] = 1 + return MPS([state1, fill(state2, n-2)..., state3]) +end + +function aklt_mps(::Type{T}, n::Int) where T + @assert n >= 2 + state1 = zeros(T, 1, 3, 3) + state1[1, 1, 1] = state1[1, 2, 2] = state1[1, 3, 3] = 1/sqrt(3) + state2 = zeros(T, 3, 3, 3) + state2[1, 1, 1] = state2[2, 2, 2] = state2[3, 3, 3] = 1 + state3 = zeros(T, 3, 3, 1) + state3[1, 1, 1] = state3[2, 2, 1] = state3[3, 3, 1] = 1 + return MPS([state1, fill(state2, n-2)..., state3]) end \ No newline at end of file diff --git a/test/mpo.jl b/test/mpo.jl index ee7e327..a57be6e 100644 --- a/test/mpo.jl +++ b/test/mpo.jl @@ -1,4 +1,4 @@ -using Test, SimpleTDVP +using Test, SimpleTDVP, Yao @testset "MPO conversion" begin m = randn(ComplexF64, 2^5, 2^5) @@ -35,4 +35,12 @@ end @test num_of_elements(mpo2) == 1*2*2*4 + 4*2*2*10 + 10*2*2*10 + 10*2*2*4 + 4*2*2*1 @test mat(mpo2) ≈ mat(mpo) @test SimpleTDVP.check_canonical(mpo2) +end + +@testset "heisenberg" begin + n = 7 + mpo = heisenberg_mpo(ComplexF64, n) + @test nsite(mpo) == n + h = EasyBuild.heisenberg(n; periodic=false) + @test mat(mpo) ≈ mat(h) end \ No newline at end of file diff --git a/test/mps.jl b/test/mps.jl index 8d046a8..ba06027 100644 --- a/test/mps.jl +++ b/test/mps.jl @@ -57,4 +57,14 @@ end @test num_of_elements(mps2) == 1*2*2 + 2*2*4 + 4*2*7 + 7*2*7 + 7*2*4 + 4*2*2 + 2*2*1 @test vec(mps2) ≈ vec(mps) @test SimpleTDVP.check_canonical(mps2) +end + +@testset "special mps" begin + mps = ghz_mps(ComplexF64, 3) + @test vec(mps) ≈ [1, 0, 0, 0, 0, 0, 0, 1]/sqrt(2) + @test SimpleTDVP.check_canonical(mps) + tensors = [randn(ComplexF64, 2) for _=1:3] + mps = product_mps(tensors...) + @test vec(mps) ≈ kron(tensors[end:-1:1]...) + @test SimpleTDVP.check_canonical(mps) end \ No newline at end of file