diff --git a/ctm/generic/rdm_kagome.py b/ctm/generic/rdm_kagome.py new file mode 100644 index 00000000..8257bb70 --- /dev/null +++ b/ctm/generic/rdm_kagome.py @@ -0,0 +1,1555 @@ +import torch +#from ipeps.ipeps import IPEPS +from ctm.generic.env import ENV +from tn_interface import contract, einsum +from tn_interface import contiguous, view, permute +from tn_interface import conj + + +def _sym_pos_def_matrix(rdm, sym_pos_def=False, verbosity=0, who="unknown"): + rdm_asym= 0.5*(rdm-rdm.conj().t()) + rdm= 0.5*(rdm+rdm.conj().t()) + if verbosity>0: + log.info(f"{who} norm(rdm_sym) {rdm.norm()} norm(rdm_asym) {rdm_asym.norm()}") + if sym_pos_def: + with torch.no_grad(): + D, U= torch.symeig(rdm, eigenvectors=True) + if D.min() < 0: + log.info(f"{who} max(diag(rdm)) {D.max()} min(diag(rdm)) {D.min()}") + D= torch.clamp(D, min=0) + rdm_posdef= U@torch.diag(D)@U.conj().t() + rdm.copy_(rdm_posdef) + rdm = rdm / rdm.diagonal().sum() + return rdm + + +def _sym_pos_def_rdm(rdm, sym_pos_def=False, verbosity=0, who=None): + assert len(rdm.size())%2==0, "invalid rank of RDM" + nsites= len(rdm.size())//2 + + orig_shape= rdm.size() + rdm= rdm.reshape(torch.prod(torch.as_tensor(rdm.size())[:nsites]),-1) + + rdm= _sym_pos_def_matrix(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + rdm= rdm.reshape(orig_shape) + return rdm + + +def rdm1x1(coord, state, env, sym_pos_def=False, verbosity=0): + r""" + :param coord: vertex (x,y) for which reduced density matrix is constructed + :param state: underlying wavefunction + :param env: environment corresponding to ``state`` + :param verbosity: logging verbosity + :type coord: tuple(int,int) + :type state: IPEPS + :type env: ENV + :type verbosity: int + :return: 1-site reduced density matrix with indices :math:`s;s'` + :rtype: torch.tensor + + Computes 1-site reduced density matrix :math:`\rho_{1x1}` centered on vertex ``coord`` by + contracting the following tensor network:: + + C--T-----C + | | | + T--A^+A--T + | | | + C--T-----C + + where the physical indices `s` and `s'` of on-site tensor :math:`A` at vertex ``coord`` + and it's hermitian conjugate :math:`A^\dagger` are left uncontracted + """ + who= "rdm1x1" + # C(-1,-1)--1->0 + # 0 + # 0 + # T(-1,0)--2 + # 1 + rdm = contract(env.C[(coord,(-1,-1))],env.T[(coord,(-1,0))],([0],[0])) + if verbosity>0: + print("rdm=CT "+str(rdm.size())) + # C(-1,-1)--0 + # | + # T(-1,0)--2->1 + # 1 + # 0 + # C(-1,1)--1->2 + rdm = contract(rdm,env.C[(coord,(-1,1))],([1],[0])) + if verbosity>0: + print("rdm=CTC "+str(rdm.size())) + # C(-1,-1)--0 + # | + # T(-1,0)--1 + # | 0->2 + # C(-1,1)--2 1--T(0,1)--2->3 + rdm = contract(rdm,env.T[(coord,(0,1))],([2],[1])) + if verbosity>0: + print("rdm=CTCT "+str(rdm.size())) + # TODO - more efficent contraction with uncontracted-double-layer on-site tensor + # Possibly reshape indices 1,2 of rdm, which are to be contracted with + # on-site tensor and contract bra,ket in two steps instead of creating + # double layer tensor + # / + # --A-- + # /|s + # + # s'|/ + # --A-- + # / + # + dimsA = state.site(coord).size() + a = contiguous(einsum('mefgh,nabcd->eafbgchdmn',state.site(coord),conj(state.site(coord)))) + a = view(a, (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + # C(-1,-1)--0 + # | + # | 0->2 + # T(-1,0)--1 1--a--3 + # | 2\45(s,s') + # | 2 + # C(-1,1)-------T(0,1)--3->1 + rdm = contract(rdm,a,([1,2],[1,2])) + if verbosity>0: + print("rdm=CTCTa "+str(rdm.size())) + # C(-1,-1)--0 0--T(0,-1)--2->0 + # | 1 + # | 2 + # T(-1,0)--------a--3->2 + # | |\45->34(s,s') + # | | + # C(-1,1)--------T(0,1)--1 + rdm = contract(env.T[(coord,(0,-1))],rdm,([0,1],[0,2])) + if verbosity>0: + print("rdm=CTCTaT "+str(rdm.size())) + # C(-1,-1)--T(0,-1)--0 0--C(1,-1) + # | | 1->0 + # | | + # T(-1,0)---a--2 + # | |\34(s,s') + # | | + # C(-1,1)---T(0,1)--0->1 + rdm = contract(env.C[(coord,(1,-1))],rdm,([0],[0])) + if verbosity>0: + print("rdm=CTCTaTC "+str(rdm.size())) + # C(-1,-1)--T(0,-1)-----C(1,-1) + # | | 0 + # | | 0 + # T(-1,0)---a--2 1------T(1,0) + # | |\34->23(s,s') 2->0 + # | | + # C(-1,1)---T(0,1)--1 + rdm = contract(env.T[(coord,(1,0))],rdm,([0,1],[0,2])) + if verbosity>0: + print("rdm=CTCTaTCT "+str(rdm.size())) + # C(-1,-1)--T(0,-1)--------C(1,-1) + # | | | + # | | | + # T(-1,0)---a--------------T(1,0) + # | |\23->12(s,s') 0 + # | | 0 + # C(-1,1)---T(0,1)--1 1----C(1,1) + rdm = contract(rdm,env.C[(coord,(1,1))],([0,1],[0,1])) + if verbosity>0: + print("rdm=CTCTaTCTC "+str(rdm.size())) + + # symmetrize and normalize + rdm= _sym_pos_def_rdm(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + + return rdm + + +# def reduced_density_matrix_kagome(rdm, phys_dim, site_types=('A', 'B', 'C')): +# idp = torch.eye(phys_dim, dtype=rdm.dtype) +# # The final two indices are physical indices! +# dims_rdm = rdm.size() +# reduced_phys_dim = dims_rdm[-1] +# if phys_dim**3 != dims_rdm[-1]: +# raise Exception("Physical dimensions do not agree. 1 site: {}, 1 unit cell: {}".format(phys_dim, dims_rdm[-1])) +# rdm = rdm.view(*dims_rdm[:-2], phys_dim, phys_dim, phys_dim, phys_dim, phys_dim, phys_dim) +# if 'A' in site_types: +# if 'B' in site_types: +# if 'C' in site_types: # do nothing +# rdm = rdm.view(*dims_rdm[:-2], phys_dim ** 3, phys_dim ** 3) +# else: # trace out C +# reduced_phys_dim //= phys_dim +# rdm = torch.einsum('...ijkabc,kc->...ijab', rdm, idp).view(*dims_rdm[:-2], phys_dim ** 2, phys_dim ** 2) +# else: +# reduced_phys_dim //= phys_dim +# if 'C' in site_types: # trace out B +# rdm = torch.einsum('...ijkabc,jb->...ikac', rdm, idp).view(*dims_rdm[:-2], phys_dim ** 2, phys_dim ** 2) +# else: # trace out B & C +# reduced_phys_dim //= phys_dim +# rdm = torch.einsum('...ijkabc,jb,kc->...ia', rdm, idp, idp) +# else: +# reduced_phys_dim //= phys_dim +# if 'B' in site_types: +# if 'C' in site_types: # trace out A +# rdm = torch.einsum('...ijkabc,ia->...jkbc', rdm, idp).view(*dims_rdm[:-2], phys_dim ** 2, phys_dim ** 2) +# else: # trace out A & C +# reduced_phys_dim //= phys_dim +# rdm = torch.einsum('...ijkabc,ia,kc->...jb', rdm, idp, idp) +# else: +# reduced_phys_dim //= phys_dim +# if 'C' in site_types: # trace out A & B +# rdm = torch.einsum('...ijkabc,ia,jb->...kc', rdm, idp, idp) +# else: # trace out A & B & C -> physical dimension = 0 +# reduced_phys_dim //= phys_dim +# rdm = torch.einsum('...ijkabc,ia,jb,kc->...', rdm, idp, idp, idp) +# rdm = rdm.view(*rdm.size(), 1, 1) +# +# return rdm, reduced_phys_dim + + +def build_reduced_density_matrix_kagome(coord, state, site_types=('A', 'B', 'C')): + r""" + :type coord: tuple(int,int) + :type state: IPEPS_KAGOME. State (rank-5 tensor) of the unit cell of kagome lattice + :param site_types: types of kagome sites needed + :rtype: torch.tensor + :return: reduced density matrix for the specified kagome sites + + Build the two-layer reduced density matrix with specified physical degrees of freedom from + the unit cell tensor \'site_tensor\'. + """ + phys_dim = state.get_physical_dim() + dims_site = state.site(coord).size() + reduced_phys_dim = dims_site[0] + if phys_dim**3 != reduced_phys_dim: + raise Exception("Physical dimensions do not agree. 1 site: {}, 1 unit cell: {}".format(phys_dim, reduced_phys_dim)) + tmp_site = state.site(coord) + idp = torch.eye(phys_dim, dtype=state.dtype) / phys_dim + if 'A' in site_types: + if 'B' in site_types: + if 'C' in site_types: # do nothing + rdm = torch.einsum('mefgh,nabcd->eafbgchdmn', tmp_site, conj(tmp_site)) + else: # trace out C + reduced_phys_dim //= phys_dim + tmp_site = state.site(coord).view(phys_dim**2, phys_dim, *dims_site[1:]) + rdm = einsum('miefgh,njabcd,ij->eafbgchdmn', tmp_site, conj(tmp_site), idp) + else: + reduced_phys_dim //= phys_dim + if 'C' in site_types: # trace out B + tmp_site = state.site(coord).view(phys_dim, phys_dim, phys_dim, *dims_site[1:]) + rdm = torch.einsum('ijkefgh,lmnabcd,jm->eafbgchdikln', tmp_site, conj(tmp_site), idp).view(*dims_site[1:], *dims_site[1:], phys_dim**2, phys_dim**2) + else: # trace out B & C + reduced_phys_dim //= phys_dim + tmp_site = state.site(coord).view(phys_dim, phys_dim, phys_dim, *dims_site[1:]) + rdm = einsum('mijefgh,nklabcd,ik,jl->eafbgchdmn', tmp_site, conj(tmp_site), idp, idp) + else: + reduced_phys_dim //= phys_dim + if 'B' in site_types: + if 'C' in site_types: # trace out A + tmp_site = state.site(coord).view(phys_dim, phys_dim**2, *dims_site[1:]) + rdm = einsum('imefgh,jnabcd,ij->eafbgchdmn', tmp_site, conj(tmp_site), idp) + else: # trace out A & C + reduced_phys_dim //= phys_dim + tmp_site = state.site(coord).view(phys_dim, phys_dim, phys_dim, *dims_site[1:]) + rdm = einsum('imjefgh,knlabcd,ik,jl->eafbgchdmn', tmp_site, conj(tmp_site), idp, idp) + else: + reduced_phys_dim //= phys_dim + if 'C' in site_types: # trace out A & B + tmp_site = state.site(coord).view(phys_dim, phys_dim, phys_dim, *dims_site[1:]) + rdm = einsum('ijmefgh,klnabcd,ik,jl->eafbgchdmn', tmp_site, conj(tmp_site), idp, idp) + else: # trace out A & B & C -> physical dimension = 0 + reduced_phys_dim //= phys_dim + tmp_site = state.site(coord).view(phys_dim, phys_dim, phys_dim, *dims_site[1:]) + rdm = einsum('ijkefgh,lmnabcd,il,jm,kn->eafbgchd', tmp_site, conj(tmp_site), idp, idp, idp) + rdm = rdm.view(*rdm.size(), 1, 1) + + return rdm, reduced_phys_dim + + +def rdm1x1_kagome(coord, state, env, sites_to_keep=('A', 'B', 'C'), sym_pos_def=False, verbosity=0): + r""" + :param coord: vertex (x,y) for which reduced density matrix is constructed + :param state: underlying wavefunction + :param env: environment corresponding to ``state`` + :param verbosity: logging verbosity + :param sites_to_keep: physical degrees of freedom to be kept. Default: "ABC" - keep all the DOF + :type coord: tuple(int,int) + :type state: IPEPS_KAGOME + :type env: ENV + :type verbosity: int + :return: 1-site reduced density matrix with indices :math:`s;s'` + :rtype: torch.tensor + + Compute 1-kagome-site reduced density matrix :math:`\rho{1x1}_A` centered on vertex ``coord``. + Inherited from the rdm1x1() method. + """ + who= "rdm1x1_kagome" + # C(-1,-1)--1->0 + # 0 + # 0 + # T(-1,0)--2 + # 1 + rdm = contract(env.C[(coord,(-1,-1))],env.T[(coord,(-1,0))],([0],[0])) + if verbosity>0: + print("rdm=CT "+str(rdm.size())) + # C(-1,-1)--0 + # | + # T(-1,0)--2->1 + # 1 + # 0 + # C(-1,1)--1->2 + rdm = contract(rdm,env.C[(coord,(-1,1))],([1],[0])) + if verbosity>0: + print("rdm=CTC "+str(rdm.size())) + # C(-1,-1)--0 + # | + # T(-1,0)--1 + # | 0->2 + # C(-1,1)--2 1--T(0,1)--2->3 + rdm = contract(rdm,env.T[(coord,(0,1))],([2],[1])) + if verbosity>0: + print("rdm=CTCT "+str(rdm.size())) + # TODO - more efficent contraction with uncontracted-double-layer on-site tensor + # Possibly reshape indices 1,2 of rdm, which are to be contracted with + # on-site tensor and contract bra,ket in two steps instead of creating + # double layer tensor + # / + # --A-- + # /|s + # + # s'|/ + # --A-- + # / + # + dimsA = state.site(coord).size() + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep) + a = view(contiguous(a), (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, rpd, rpd)) + + # C(-1,-1)--0 + # | + # | 0->2 + # T(-1,0)--1 1--a--3 + # | 2\45(s,s') + # | 2 + # C(-1,1)-------T(0,1)--3->1 + rdm = contract(rdm,a,([1,2],[1,2])) + if verbosity>0: + print("rdm=CTCTa "+str(rdm.size())) + # C(-1,-1)--0 0--T(0,-1)--2->0 + # | 1 + # | 2 + # T(-1,0)--------a--3->2 + # | |\45->34(s,s') + # | | + # C(-1,1)--------T(0,1)--1 + rdm = contract(env.T[(coord,(0,-1))],rdm,([0,1],[0,2])) + if verbosity>0: + print("rdm=CTCTaT "+str(rdm.size())) + # C(-1,-1)--T(0,-1)--0 0--C(1,-1) + # | | 1->0 + # | | + # T(-1,0)---a--2 + # | |\34(s,s') + # | | + # C(-1,1)---T(0,1)--0->1 + rdm = contract(env.C[(coord,(1,-1))],rdm,([0],[0])) + if verbosity>0: + print("rdm=CTCTaTC "+str(rdm.size())) + # C(-1,-1)--T(0,-1)-----C(1,-1) + # | | 0 + # | | 0 + # T(-1,0)---a--2 1------T(1,0) + # | |\34->23(s,s') 2->0 + # | | + # C(-1,1)---T(0,1)--1 + rdm = contract(env.T[(coord,(1,0))],rdm,([0,1],[0,2])) + if verbosity>0: + print("rdm=CTCTaTCT "+str(rdm.size())) + # C(-1,-1)--T(0,-1)--------C(1,-1) + # | | | + # | | | + # T(-1,0)---a--------------T(1,0) + # | |\23->12(s,s') 0 + # | | 0 + # C(-1,1)---T(0,1)--1 1----C(1,1) + rdm = contract(rdm,env.C[(coord,(1,1))],([0,1],[0,1])) + if verbosity>0: + print("rdm=CTCTaTCTC "+str(rdm.size())) + + # symmetrize and normalize + rdm= _sym_pos_def_rdm(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + + return rdm + + +def rdm2x1(coord, state, env, sym_pos_def=False, verbosity=0): + r""" + :param coord: vertex (x,y) specifies position of 2x1 subsystem + :param state: underlying wavefunction + :param env: environment corresponding to ``state`` + :param verbosity: logging verbosity + :type coord: tuple(int,int) + :type state: IPEPS + :type env: ENV + :type verbosity: int + :return: 2-site reduced density matrix with indices :math:`s_0s_1;s'_0s'_1` + :rtype: torch.tensor + + Computes 2-site reduced density matrix :math:`\rho_{2x1}` of a horizontal + 2x1 subsystem using following strategy: + + 1. compute four individual corners + 2. construct right and left half of the network + 3. contract right and left halt to obtain final reduced density matrix + + :: + + | | | | | | + T--A^+A(coord)--A^+A(coord+(1,0))--T C2x1_LD(coord)--C2x1(coord+(1,0)) + | | | | + C--T------------T------------------C + + The physical indices `s` and `s'` of on-sites tensors :math:`A` (and :math:`A^\dagger`) + at vertices ``coord``, ``coord+(1,0)`` are left uncontracted + """ + who="rdm2x1" + #----- building C2x2_LU ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord),(-1,-1))] + T1 = env.T[(state.vertexToSite(coord),(0,-1))] + T2 = env.T[(state.vertexToSite(coord),(-1,0))] + dimsA = state.site(coord).size() + a = einsum('mefgh,nabcd->eafbgchdmn',state.site(coord),conj(state.site(coord))) + a = view(contiguous(a),\ + (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + # C--10--T1--2 + # 0 1 + C2x2_LU =contract(C, T1, ([1],[0])) + + # C------T1--2->1 + # 0 1->0 + # 0 + # T2--2->3 + # 1->2 + C2x2_LU =contract(C2x2_LU, T2, ([0],[0])) + + # C-------T1--1->0 + # | 0 + # | 0 + # T2--3 1 a--3 + # 2->1 2\45 + C2x2_LU =contract(C2x2_LU, a, ([0,3],[0,1])) + + # permute 012345->120345 + # reshape (12)(03)45->0123 + # C2x2--1 + # |\23 + # 0 + C2x2_LU= permute(C2x2_LU, (1,2,0,3,4,5)) + C2x2_LU= view(contiguous(C2x2_LU), \ + (T2.size(1)*a.size(2),T1.size(2)*a.size(3),dimsA[0],dimsA[0])) + if verbosity>0: + print("C2X2 LU "+str(coord)+"->"+str(state.vertexToSite(coord))+" (-1,-1): "+str(C2x2_LU.size())) + + #----- building C2x1_LD ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord),(-1,1))] + T2 = env.T[(state.vertexToSite(coord),(0,1))] + + # 0 0->1 + # C--1 1--T2--2 + C2x1_LD=contract(C, T2, ([1],[1])) + + # reshape (01)2->(0)1 + # 0 + # | + # C2x1--1 + C2x1_LD= view(contiguous(C2x1_LD), (C.size(0)*T2.size(0),T2.size(2))) + if verbosity>0: + print("C2X1 LD "+str(coord)+"->"+str(state.vertexToSite(coord))+" (-1,1): "+str(C2x1_LD.size())) + + #----- build left part C2x2_LU--C2x1_LD ------------------------------------ + # C2x2_LU--1 + # |\23 + # 0 + # 0 + # C2x1_LD--1->0 + # TODO is it worthy(performance-wise) to instead overwrite one of C2x2_LU,C2x2_RU ? + left_half= contract(C2x1_LD, C2x2_LU, ([0],[0])) + + #----- building C2x2_RU ---------------------------------------------------- + vec = (1,0) + shitf_coord = state.vertexToSite((coord[0]+vec[0],coord[1]+vec[1])) + C = env.C[(shitf_coord,(1,-1))] + T1 = env.T[(shitf_coord,(1,0))] + T2 = env.T[(shitf_coord,(0,-1))] + dimsA = state.site(shitf_coord).size() + a= einsum('mefgh,nabcd->eafbgchdmn',state.site(shitf_coord),conj(state.site(shitf_coord))) + a= view(contiguous(a), \ + (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + + # 0--C + # 1 + # 0 + # 1--T1 + # 2 + C2x2_RU =contract(C, T1, ([1],[0])) + + # 2<-0--T2--2 0--C + # 3<-1 | + # 0<-1--T1 + # 1<-2 + C2x2_RU =contract(C2x2_RU, T2, ([0],[2])) + + # 1<-2--T2------C + # 3 | + # 45\0 | + # 2<-1--a--3 0--T1 + # 3<-2 0<-1 + C2x2_RU =contract(C2x2_RU, a, ([0,3],[3,0])) + + # permute 012334->120345 + # reshape (12)(03)45->0123 + # 0--C2x2 + # 23/| + # 1 + C2x2_RU= permute(C2x2_RU, (1,2,0,3,4,5)) + C2x2_RU= view(contiguous(C2x2_RU), \ + (T2.size(0)*a.size(1),T1.size(2)*a.size(2), dimsA[0], dimsA[0])) + if verbosity>0: + print("C2X2 RU "+str((coord[0]+vec[0],coord[1]+vec[1]))+"->"+str(shitf_coord)+" (1,-1): "+str(C2x2_RU.size())) + + #----- building C2x1_RD ---------------------------------------------------- + C = env.C[(shitf_coord,(1,1))] + T1 = env.T[(shitf_coord,(0,1))] + + # 1<-0 0 + # 2<-1--T1--2 1--C + C2x1_RD =contract(C, T1, ([1],[2])) + + # reshape (01)2->(0)1 + C2x1_RD = view(contiguous(C2x1_RD), (C.size(0)*T1.size(0),T1.size(1))) + + # 0 + # | + # 1--C2x1 + if verbosity>0: + print("C2X1 RD "+str((coord[0]+vec[0],coord[1]+vec[1]))+"->"+str(shitf_coord)+" (1,1): "+str(C2x1_RD.size())) + + + + #----- build right part C2x2_RU--C2x1_RD ----------------------------------- + # 1<-0--C2x2_RU + # |\23 + # 1 + # 0 + # 0<-1--C2x1_RD + right_half =contract(C2x1_RD, C2x2_RU, ([0],[1])) + + # construct reduced density matrix by contracting left and right halfs + # C2x2_LU--1 1----C2x2_RU + # |\23->01 |\23 + # | | + # C2x1_LD--0 0----C2x1_RD + rdm =contract(left_half,right_half,([0,1],[0,1])) + + # permute into order of s0,s1;s0',s1' where primed indices + # represent "ket" + # 0123->0213 + # symmetrize and normalize + rdm = contiguous(permute(rdm, (0,2,1,3))) + rdm= _sym_pos_def_rdm(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + + return rdm + + +def rdm2x1_kagome(coord, state, env, sites_to_keep_00=('A', 'B', 'C'), sites_to_keep_10=('A', 'B', 'C'), sym_pos_def=False, verbosity=0): + r""" + :param coord: vertex (x,y) specifies position of 2x1 subsystem + :param state: underlying wavefunction + :param env: environment corresponding to ``state`` + :param verbosity: logging verbosity + :param sites_to_keep_00: physical sites needed for the unit cell at coord + (0, 0) + :param sites_to_keep_10: physical sites needed for the unit cell at coord + (1, 0) + :type coord: tuple(int,int) + :type state: IPEPS_KAGOME + :type env: ENV + :type verbosity: int + :return: 2-site reduced density matrix with indices :math:`s_0s_1;s'_0s'_1` + :rtype: torch.tensor + + Computes 2-site reduced density matrix :math:`\rho_{2x1}` of a horizontal + 2x1 subsystem using following strategy: + + 1. compute four individual corners + 2. construct right and left half of the network + 3. contract right and left halt to obtain final reduced density matrix + + :: + + C--T------------T------------------C = C2x2_LU(coord)--C2x2(coord+(1,0)) + | | | | | | + T--A^+A(coord)--A^+A(coord+(1,0))--T C2x1_LD(coord)--C2x1(coord+(1,0)) + | | | | + C--T------------T------------------C + + The physical indices `s` and `s'` of on-sites tensors :math:`A` (and :math:`A^\dagger`) + at vertices ``coord``, ``coord+(1,0)`` are left uncontracted + """ + who = "rdm2x1_kagome" + # ----- building C2x2_LU ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord), (-1, -1))] + T1 = env.T[(state.vertexToSite(coord), (0, -1))] + T2 = env.T[(state.vertexToSite(coord), (-1, 0))] + dimsA = state.site(coord).size() + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep_00) + a = view(contiguous(a), (dimsA[1] ** 2, dimsA[2] ** 2, dimsA[3] ** 2, dimsA[4] ** 2, rpd, rpd)) + # C--10--T1--2 + # 0 1 + C2x2_LU = contract(C, T1, ([1], [0])) + + # C------T1--2->1 + # 0 1->0 + # 0 + # T2--2->3 + # 1->2 + C2x2_LU = contract(C2x2_LU, T2, ([0], [0])) + + # C-------T1--1->0 + # | 0 + # | 0 + # T2--3 1 a--3 + # 2->1 2\45 + C2x2_LU = contract(C2x2_LU, a, ([0, 3], [0, 1])) + + # permute 012345->120345 + # reshape (12)(03)45->0123 + # C2x2--1 + # |\23 + # 0 + C2x2_LU = permute(C2x2_LU, (1, 2, 0, 3, 4, 5)) + # dimsA[0] -> rpd. Modified by Yi. Similar for other functions + C2x2_LU = view(contiguous(C2x2_LU), (T2.size(1) * a.size(2), T1.size(2) * a.size(3), rpd, rpd)) + if verbosity > 0: + print("C2X2 LU " + str(coord) + "->" + str(state.vertexToSite(coord)) + " (-1,-1): " + str(C2x2_LU.size())) + + # ----- building C2x1_LD ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord), (-1, 1))] + T2 = env.T[(state.vertexToSite(coord), (0, 1))] + + # 0 0->1 + # C--1 1--T2--2 + C2x1_LD = contract(C, T2, ([1], [1])) + + # reshape (01)2->(0)1 + # 0 + # | + # C2x1--1 + C2x1_LD = view(contiguous(C2x1_LD), (C.size(0) * T2.size(0), T2.size(2))) + if verbosity > 0: + print("C2X1 LD " + str(coord) + "->" + str(state.vertexToSite(coord)) + " (-1,1): " + str(C2x1_LD.size())) + + # ----- build left part C2x2_LU--C2x1_LD ------------------------------------ + # C2x2_LU--1 + # |\23 + # 0 + # 0 + # C2x1_LD--1->0 + # TODO is it worthy(performance-wise) to instead overwrite one of C2x2_LU,C2x2_RU ? + left_half = contract(C2x1_LD, C2x2_LU, ([0], [0])) + + # ----- building C2x2_RU ---------------------------------------------------- + vec = (1, 0) + shitf_coord = state.vertexToSite((coord[0] + vec[0], coord[1] + vec[1])) + C = env.C[(shitf_coord, (1, -1))] + T1 = env.T[(shitf_coord, (1, 0))] + T2 = env.T[(shitf_coord, (0, -1))] + dimsA = state.site(shitf_coord).size() + + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep_10) + a = view(contiguous(a), (dimsA[1] ** 2, dimsA[2] ** 2, dimsA[3] ** 2, dimsA[4] ** 2, rpd, rpd)) + # 0--C + # 1 + # 0 + # 1--T1 + # 2 + C2x2_RU = contract(C, T1, ([1], [0])) + + # 2<-0--T2--2 0--C + # 3<-1 | + # 0<-1--T1 + # 1<-2 + C2x2_RU = contract(C2x2_RU, T2, ([0], [2])) + + # 1<-2--T2------C + # 3 | + # 45\0 | + # 2<-1--a--3 0--T1 + # 3<-2 0<-1 + C2x2_RU = contract(C2x2_RU, a, ([0, 3], [3, 0])) + + # permute 012334->120345 + # reshape (12)(03)45->0123 + # 0--C2x2 + # 23/| + # 1 + C2x2_RU = permute(C2x2_RU, (1, 2, 0, 3, 4, 5)) + C2x2_RU = view(contiguous(C2x2_RU), \ + (T2.size(0) * a.size(1), T1.size(2) * a.size(2), rpd, rpd)) + if verbosity > 0: + print("C2X2 RU " + str((coord[0] + vec[0], coord[1] + vec[1])) + "->" + str(shitf_coord) + " (1,-1): " + str( + C2x2_RU.size())) + + # ----- building C2x1_RD ---------------------------------------------------- + C = env.C[(shitf_coord, (1, 1))] + T1 = env.T[(shitf_coord, (0, 1))] + + # 1<-0 0 + # 2<-1--T1--2 1--C + C2x1_RD = contract(C, T1, ([1], [2])) + + # reshape (01)2->(0)1 + C2x1_RD = view(contiguous(C2x1_RD), (C.size(0) * T1.size(0), T1.size(1))) + + # 0 + # | + # 1--C2x1 + if verbosity > 0: + print("C2X1 RD " + str((coord[0] + vec[0], coord[1] + vec[1])) + "->" + str(shitf_coord) + " (1,1): " + str( + C2x1_RD.size())) + + # ----- build right part C2x2_RU--C2x1_RD ----------------------------------- + # 1<-0--C2x2_RU + # |\23 + # 1 + # 0 + # 0<-1--C2x1_RD + right_half = contract(C2x1_RD, C2x2_RU, ([0], [1])) + + # construct reduced density matrix by contracting left and right halfs + # C2x2_LU--1 1----C2x2_RU + # |\23->01 |\23 + # | | + # C2x1_LD--0 0----C2x1_RD + rdm = contract(left_half, right_half, ([0, 1], [0, 1])) + + # permute into order of s0,s1;s0',s1' where primed indices + # represent "ket" + # 0123->0213 + # symmetrize and normalize + rdm = contiguous(permute(rdm, (0, 2, 1, 3))) + rdm = _sym_pos_def_rdm(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + + return rdm + + +def rdm1x2(coord, state, env, sym_pos_def=False, verbosity=0): + r""" + :param coord: vertex (x,y) specifies position of 1x2 subsystem + :param state: underlying wavefunction + :param env: environment corresponding to ``state`` + :param verbosity: logging verbosity + :type coord: tuple(int,int) + :type state: IPEPS + :type env: ENV + :type verbosity: int + :return: 2-site reduced density matrix with indices :math:`s_0s_1;s'_0s'_1` + :rtype: torch.tensor + + Computes 2-site reduced density matrix :math:`\rho_{1x2}` of a vertical + 1x2 subsystem using following strategy: + + 1. compute four individual corners + 2. construct upper and lower half of the network + 3. contract upper and lower halt to obtain final reduced density matrix + + :: + + C--T------------------C = C2x2_LU(coord)--------C1x2(coord) + | | | | | + T--A^+A(coord)--------T C2x2_LD(coord+(0,1))--C1x2(coord+0,1)) + | | | + T--A^+A(coord+(0,1))--T + | | | + C--T------------------C + + The physical indices `s` and `s'` of on-sites tensors :math:`A` (and :math:`A^\dagger`) + at vertices ``coord``, ``coord+(0,1)`` are left uncontracted + """ + who="rdm1x2" + #----- building C2x2_LU ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord),(-1,-1))] + T1 = env.T[(state.vertexToSite(coord),(0,-1))] + T2 = env.T[(state.vertexToSite(coord),(-1,0))] + dimsA = state.site(coord).size() + a= einsum('mefgh,nabcd->eafbgchdmn',state.site(coord), conj(state.site(coord))) + a= view(contiguous(a), \ + (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + + # C--10--T1--2 + # 0 1 + C2x2_LU =contract(C, T1, ([1],[0])) + + # C------T1--2->1 + # 0 1->0 + # 0 + # T2--2->3 + # 1->2 + C2x2_LU =contract(C2x2_LU, T2, ([0],[0])) + + # C-------T1--1->0 + # | 0 + # | 0 + # T2--3 1 a--3 + # 2->1 2\45 + C2x2_LU =contract(C2x2_LU, a, ([0,3],[0,1])) + + # permute 012345->120345 + # reshape (12)(03)45->0123 + # C2x2--1 + # |\23 + # 0 + C2x2_LU= permute(C2x2_LU, (1,2,0,3,4,5)) + C2x2_LU= view(contiguous(C2x2_LU), \ + (T2.size(1)*a.size(2),T1.size(2)*a.size(3),dimsA[0],dimsA[0])) + if verbosity>0: + print("C2X2 LU "+str(coord)+"->"+str(state.vertexToSite(coord))+" (-1,-1): "+str(C2x2_LU.size())) + + #----- building C1x2_RU ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord),(1,-1))] + T1 = env.T[(state.vertexToSite(coord),(1,0))] + + # 0--C + # 1 + # 0 + # 1--T1 + # 2 + C1x2_RU =contract(C, T1, ([1],[0])) + + # reshape (01)2->(0)1 + # 0--C1x2 + # 23/| + # 1 + C1x2_RU= view(contiguous(C1x2_RU), (C.size(0)*T1.size(1),T1.size(2))) + if verbosity>0: + print("C1X2 RU "+str(coord)+"->"+str(state.vertexToSite(coord))+" (1,-1): "+str(C1x2_RU.size())) + + #----- build upper part C2x2_LU--C1x2_RU ----------------------------------- + # C2x2_LU--1 0--C1x2_RU + # |\23 | + # 0->1 1->0 + upper_half =contract(C1x2_RU, C2x2_LU, ([0],[1])) + + #----- building C2x2_LD ---------------------------------------------------- + vec = (0,1) + shitf_coord = state.vertexToSite((coord[0]+vec[0],coord[1]+vec[1])) + C = env.C[(shitf_coord,(-1,1))] + T1 = env.T[(shitf_coord,(-1,0))] + T2 = env.T[(shitf_coord,(0,1))] + dimsA = state.site(shitf_coord).size() + a= einsum('mefgh,nabcd->eafbgchdmn',state.site(shitf_coord),conj(state.site(shitf_coord))) + a= view(contiguous(a), \ + (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + + # 0->1 + # T1--2 + # 1 + # 0 + # C--1->0 + C2x2_LD =contract(C, T1, ([0],[1])) + + # 1->0 + # T1--2->1 + # | + # | 0->2 + # C--0 1--T2--2->3 + C2x2_LD =contract(C2x2_LD, T2, ([0],[1])) + + # 0 0->2 + # T1--1 1--a--3 + # | 2\45 + # | 2 + # C--------T2--3->1 + C2x2_LD =contract(C2x2_LD, a, ([1,2],[1,2])) + + # permute 012345->021345 + # reshape (02)(13)45->0123 + # 0 + # |/23 + # C2x2--1 + C2x2_LD= permute(C2x2_LD, (0,2,1,3,4,5)) + C2x2_LD= view(contiguous(C2x2_LD), \ + (T1.size(0)*a.size(0),T2.size(2)*a.size(3), dimsA[0], dimsA[0])) + if verbosity>0: + print("C2X2 LD "+str((coord[0]+vec[0],coord[1]+vec[1]))+"->"+str(shitf_coord)+" (-1,1): "+str(C2x2_LD.size())) + + #----- building C2x2_RD ---------------------------------------------------- + C = env.C[(shitf_coord,(1,1))] + T2 = env.T[(shitf_coord,(1,0))] + + # 0 + # 1--T2 + # 2 + # 0 + # 2<-1--C + C1x2_RD =contract(T2, C, ([2],[0])) + + # permute 012->021 + # reshape 0(12)->0(1) + C1x2_RD = view(contiguous(permute(C1x2_RD,(0,2,1))), \ + (T2.size()[0],C.size()[1]*T2.size()[1])) + + # 0 + # | + # 1--C1x2 + if verbosity>0: + print("C1X2 RD "+str((coord[0]+vec[0],coord[1]+vec[1]))+"->"+str(shitf_coord)+" (1,1): "+str(C1x2_RD.size())) + + #----- build lower part C2x2_LD--C1x2_RD ----------------------------------- + # 0->1 0 + # |/23 | + # C2x2_LD--1 1--C1x2_RD + lower_half =contract(C1x2_RD, C2x2_LD, ([1],[1])) + + # construct reduced density matrix by contracting lower and upper halfs + # C2x2_LU------C1x2_RU + # |\23->01 | + # 1 0 + # 1 0 + # |/23 | + # C2x2_LD------C1x2_RD + rdm =contract(upper_half,lower_half,([0,1],[0,1])) + + # permute into order of s0,s1;s0',s1' where primed indices + # represent "ket" + # 0123->0213 + # symmetrize and normalize + rdm = contiguous(permute(rdm, (0,2,1,3))) + rdm= _sym_pos_def_rdm(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + + return rdm + + +def rdm1x2_kagome(coord, state, env, sites_to_keep_00=('A', 'B', 'C'), sites_to_keep_01=('A', 'B', 'C'), sym_pos_def=False, verbosity=0): + r""" + :param coord: vertex (x,y) specifies position of 1x2 subsystem + :param state: underlying wavefunction + :param env: environment corresponding to ``state`` + :param verbosity: logging verbosity + :param sites_to_keep_00: physical sites needed for the unit cell at coord + (0, 0) + :param sites_to_keep_01: physical sites needed for the unit cell at coord + (0, 1) + :type coord: tuple(int,int) + :type state: IPEPS + :type env: ENV + :type verbosity: int + :return: 2-site reduced density matrix with indices :math:`s_0s_1;s'_0s'_1` + :rtype: torch.tensor + + Computes 2-site reduced density matrix :math:`\rho_{1x2}` of a vertical + 1x2 subsystem using following strategy: + """ + who = "rdm1x2_kagome" + # ----- building C2x2_LU ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord), (-1, -1))] + T1 = env.T[(state.vertexToSite(coord), (0, -1))] + T2 = env.T[(state.vertexToSite(coord), (-1, 0))] + dimsA = state.site(coord).size() + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep_00) + a = view(contiguous(a), (dimsA[1] ** 2, dimsA[2] ** 2, dimsA[3] ** 2, dimsA[4] ** 2, rpd, rpd)) + + # C--10--T1--2 + # 0 1 + C2x2_LU = contract(C, T1, ([1], [0])) + + # C------T1--2->1 + # 0 1->0 + # 0 + # T2--2->3 + # 1->2 + C2x2_LU = contract(C2x2_LU, T2, ([0], [0])) + + # C-------T1--1->0 + # | 0 + # | 0 + # T2--3 1 a--3 + # 2->1 2\45 + C2x2_LU = contract(C2x2_LU, a, ([0, 3], [0, 1])) + + # permute 012345->120345 + # reshape (12)(03)45->0123 + # C2x2--1 + # |\23 + # 0 + C2x2_LU = permute(C2x2_LU, (1, 2, 0, 3, 4, 5)) + C2x2_LU = view(contiguous(C2x2_LU), (T2.size(1) * a.size(2), T1.size(2) * a.size(3), rpd, rpd)) + if verbosity > 0: + print("C2X2 LU " + str(coord) + "->" + str(state.vertexToSite(coord)) + " (-1,-1): " + str(C2x2_LU.size())) + + # ----- building C1x2_RU ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord), (1, -1))] + T1 = env.T[(state.vertexToSite(coord), (1, 0))] + + # 0--C + # 1 + # 0 + # 1--T1 + # 2 + C1x2_RU = contract(C, T1, ([1], [0])) + + # reshape (01)2->(0)1 + # 0--C1x2 + # 23/| + # 1 + C1x2_RU = view(contiguous(C1x2_RU), (C.size(0) * T1.size(1), T1.size(2))) + if verbosity > 0: + print("C1X2 RU " + str(coord) + "->" + str(state.vertexToSite(coord)) + " (1,-1): " + str(C1x2_RU.size())) + + # ----- build upper part C2x2_LU--C1x2_RU ----------------------------------- + # C2x2_LU--1 0--C1x2_RU + # |\23 | + # 0->1 1->0 + upper_half = contract(C1x2_RU, C2x2_LU, ([0], [1])) + + # ----- building C2x2_LD ---------------------------------------------------- + vec = (0, 1) + shitf_coord = state.vertexToSite((coord[0] + vec[0], coord[1] + vec[1])) + C = env.C[(shitf_coord, (-1, 1))] + T1 = env.T[(shitf_coord, (-1, 0))] + T2 = env.T[(shitf_coord, (0, 1))] + dimsA = state.site(shitf_coord).size() + + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep_01) + a = view(contiguous(a), (dimsA[1] ** 2, dimsA[2] ** 2, dimsA[3] ** 2, dimsA[4] ** 2, rpd, rpd)) + # 0->1 + # T1--2 + # 1 + # 0 + # C--1->0 + C2x2_LD = contract(C, T1, ([0], [1])) + + # 1->0 + # T1--2->1 + # | + # | 0->2 + # C--0 1--T2--2->3 + C2x2_LD = contract(C2x2_LD, T2, ([0], [1])) + + # 0 0->2 + # T1--1 1--a--3 + # | 2\45 + # | 2 + # C--------T2--3->1 + C2x2_LD = contract(C2x2_LD, a, ([1, 2], [1, 2])) + + # permute 012345->021345 + # reshape (02)(13)45->0123 + # 0 + # |/23 + # C2x2--1 + C2x2_LD = permute(C2x2_LD, (0, 2, 1, 3, 4, 5)) + C2x2_LD = view(contiguous(C2x2_LD), (T1.size(0) * a.size(0), T2.size(2) * a.size(3), rpd, rpd)) + if verbosity > 0: + print("C2X2 LD " + str((coord[0] + vec[0], coord[1] + vec[1])) + "->" + str(shitf_coord) + " (-1,1): " + str( + C2x2_LD.size())) + + # ----- building C2x2_RD ---------------------------------------------------- + C = env.C[(shitf_coord, (1, 1))] + T2 = env.T[(shitf_coord, (1, 0))] + + # 0 + # 1--T2 + # 2 + # 0 + # 2<-1--C + C1x2_RD = contract(T2, C, ([2], [0])) + + # permute 012->021 + # reshape 0(12)->0(1) + C1x2_RD = view(contiguous(permute(C1x2_RD, (0, 2, 1))), \ + (T2.size()[0], C.size()[1] * T2.size()[1])) + + # 0 + # | + # 1--C1x2 + if verbosity > 0: + print("C1X2 RD " + str((coord[0] + vec[0], coord[1] + vec[1])) + "->" + str(shitf_coord) + " (1,1): " + str( + C1x2_RD.size())) + + # ----- build lower part C2x2_LD--C1x2_RD ----------------------------------- + # 0->1 0 + # |/23 | + # C2x2_LD--1 1--C1x2_RD + lower_half = contract(C1x2_RD, C2x2_LD, ([1], [1])) + + # construct reduced density matrix by contracting lower and upper halfs + # C2x2_LU------C1x2_RU + # |\23->01 | + # 1 0 + # 1 0 + # |/23 | + # C2x2_LD------C1x2_RD + rdm = contract(upper_half, lower_half, ([0, 1], [0, 1])) + + # permute into order of s0,s1;s0',s1' where primed indices + # represent "ket" + # 0123->0213 + # symmetrize and normalize + rdm = contiguous(permute(rdm, (0, 2, 1, 3))) + rdm = _sym_pos_def_rdm(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + + return rdm + + +def rdm2x2(coord, state, env, sym_pos_def=False, verbosity=0): + r""" + :param coord: vertex (x,y) specifies upper left site of 2x2 subsystem + :param state: underlying wavefunction + :param env: environment corresponding to ``state`` + :param verbosity: logging verbosity + :type coord: tuple(int,int) + :type state: IPEPS + :type env: ENV + :type verbosity: int + :return: 4-site reduced density matrix with indices :math:`s_0s_1s_2s_3;s'_0s'_1s'_2s'_3` + :rtype: torch.tensor + + Computes 4-site reduced density matrix :math:`\rho_{2x2}` of 2x2 subsystem specified + by the vertex ``coord`` of its upper left corner using strategy: + + 1. compute four individual corners + 2. construct upper and lower half of the network + 3. contract upper and lower half to obtain final reduced density matrix + + :: + + C--T------------------T------------------C = C2x2_LU(coord)--------C2x2(coord+(1,0)) + | | | | | | + T--A^+A(coord)--------A^+A(coord+(1,0))--T C2x2_LD(coord+(0,1))--C2x2(coord+(1,1)) + | | | | + T--A^+A(coord+(0,1))--A^+A(coord+(1,1))--T + | | | | + C--T------------------T------------------C + + The physical indices `s` and `s'` of on-sites tensors :math:`A` (and :math:`A^\dagger`) + at vertices ``coord``, ``coord+(1,0)``, ``coord+(0,1)``, and ``coord+(1,1)`` are + left uncontracted and given in the same order:: + + s0 s1 + s2 s3 + + """ + who= "rdm2x2" + #----- building C2x2_LU ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord),(-1,-1))] + T1 = env.T[(state.vertexToSite(coord),(0,-1))] + T2 = env.T[(state.vertexToSite(coord),(-1,0))] + dimsA = state.site(coord).size() + a = contiguous(einsum('mefgh,nabcd->eafbgchdmn',state.site(coord),conj(state.site(coord)))) + a = view(a, (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + + # C--10--T1--2 + # 0 1 + C2x2_LU = contract(C, T1, ([1],[0])) + + # C------T1--2->1 + # 0 1->0 + # 0 + # T2--2->3 + # 1->2 + C2x2_LU = contract(C2x2_LU, T2, ([0],[0])) + + # C-------T1--1->0 + # | 0 + # | 0 + # T2--3 1 a--3 + # 2->1 2\45 + C2x2_LU = contract(C2x2_LU, a, ([0,3],[0,1])) + + # permute 012345->120345 + # reshape (12)(03)45->0123 + # C2x2--1 + # |\23 + # 0 + C2x2_LU = contiguous(permute(C2x2_LU,(1,2,0,3,4,5))) + C2x2_LU = view(C2x2_LU, (T2.size(1)*a.size(2),T1.size(2)*a.size(3),dimsA[0],dimsA[0])) + if verbosity>0: + print("C2X2 LU "+str(coord)+"->"+str(state.vertexToSite(coord))+" (-1,-1): "+str(C2x2_LU.size())) + + #----- building C2x2_RU ---------------------------------------------------- + vec = (1,0) + shitf_coord = state.vertexToSite((coord[0]+vec[0],coord[1]+vec[1])) + C = env.C[(shitf_coord,(1,-1))] + T1 = env.T[(shitf_coord,(1,0))] + T2 = env.T[(shitf_coord,(0,-1))] + dimsA = state.site(shitf_coord).size() + a = contiguous(einsum('mefgh,nabcd->eafbgchdmn',state.site(shitf_coord),conj(state.site(shitf_coord)))) + a = view(a, (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + + # 0--C + # 1 + # 0 + # 1--T1 + # 2 + C2x2_RU = contract(C, T1, ([1],[0])) + + # 2<-0--T2--2 0--C + # 3<-1 | + # 0<-1--T1 + # 1<-2 + C2x2_RU = contract(C2x2_RU, T2, ([0],[2])) + + # 1<-2--T2------C + # 3 | + # 45\0 | + # 2<-1--a--3 0--T1 + # 3<-2 0<-1 + C2x2_RU = contract(C2x2_RU, a, ([0,3],[3,0])) + + # permute 012334->120345 + # reshape (12)(03)45->0123 + # 0--C2x2 + # 23/| + # 1 + C2x2_RU = contiguous(permute(C2x2_RU, (1,2,0,3,4,5))) + C2x2_RU = view(C2x2_RU, (T2.size(0)*a.size(1),T1.size(2)*a.size(2), dimsA[0], dimsA[0])) + if verbosity>0: + print("C2X2 RU "+str((coord[0]+vec[0],coord[1]+vec[1]))+"->"+str(shitf_coord)+" (1,-1): "+str(C2x2_RU.size())) + + #----- build upper part C2x2_LU--C2x2_RU ----------------------------------- + # C2x2_LU--1 0--C2x2_RU C2x2_LU------C2x2_RU + # |\23->12 |\23->45 & permute |\12->23 |\45 + # 0 1->3 0 3->1 + # TODO is it worthy(performance-wise) to instead overwrite one of C2x2_LU,C2x2_RU ? + upper_half = contract(C2x2_LU, C2x2_RU, ([1],[0])) + upper_half = permute(upper_half, (0,3,1,2,4,5)) + + #----- building C2x2_RD ---------------------------------------------------- + vec = (1,1) + shitf_coord = state.vertexToSite((coord[0]+vec[0],coord[1]+vec[1])) + C = env.C[(shitf_coord,(1,1))] + T1 = env.T[(shitf_coord,(0,1))] + T2 = env.T[(shitf_coord,(1,0))] + dimsA = state.site(shitf_coord).size() + a = contiguous(einsum('mefgh,nabcd->eafbgchdmn',state.site(shitf_coord),conj(state.site(shitf_coord)))) + a = view(a, (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + + # 1<-0 0 + # 2<-1--T1--2 1--C + C2x2_RD = contract(C, T1, ([1],[2])) + + # 2<-0 + # 3<-1--T2 + # 2 + # 0<-1 0 + # 1<-2--T1---C + C2x2_RD = contract(C2x2_RD, T2, ([0],[2])) + + # 2<-0 1<-2 + # 3<-1--a--3 3--T2 + # 2\45 | + # 0 | + # 0<-1--T1------C + C2x2_RD = contract(C2x2_RD, a, ([0,3],[2,3])) + + # permute 012345->120345 + # reshape (12)(03)45->0123 + C2x2_RD = contiguous(permute(C2x2_RD, (1,2,0,3,4,5))) + C2x2_RD = view(C2x2_RD, (T2.size(0)*a.size(0),T1.size(1)*a.size(1), dimsA[0], dimsA[0])) + + # 0 + # |/23 + # 1--C2x2 + if verbosity>0: + print("C2X2 RD "+str((coord[0]+vec[0],coord[1]+vec[1]))+"->"+str(shitf_coord)+" (1,1): "+str(C2x2_RD.size())) + + #----- building C2x2_LD ---------------------------------------------------- + vec = (0,1) + shitf_coord = state.vertexToSite((coord[0]+vec[0],coord[1]+vec[1])) + C = env.C[(shitf_coord,(-1,1))] + T1 = env.T[(shitf_coord,(-1,0))] + T2 = env.T[(shitf_coord,(0,1))] + dimsA = state.site(shitf_coord).size() + a = contiguous(einsum('mefgh,nabcd->eafbgchdmn',state.site(shitf_coord),conj(state.site(shitf_coord)))) + a = view(a, (dimsA[1]**2, dimsA[2]**2, dimsA[3]**2, dimsA[4]**2, dimsA[0], dimsA[0])) + + # 0->1 + # T1--2 + # 1 + # 0 + # C--1->0 + C2x2_LD = contract(C, T1, ([0],[1])) + + # 1->0 + # T1--2->1 + # | + # | 0->2 + # C--0 1--T2--2->3 + C2x2_LD = contract(C2x2_LD, T2, ([0],[1])) + + # 0 0->2 + # T1--1 1--a--3 + # | 2\45 + # | 2 + # C--------T2--3->1 + C2x2_LD = contract(C2x2_LD, a, ([1,2],[1,2])) + + # permute 012345->021345 + # reshape (02)(13)45->0123 + # 0 + # |/23 + # C2x2--1 + C2x2_LD = contiguous(permute(C2x2_LD, (0,2,1,3,4,5))) + C2x2_LD = view(C2x2_LD, (T1.size(0)*a.size(0),T2.size(2)*a.size(3), dimsA[0], dimsA[0])) + if verbosity>0: + print("C2X2 LD "+str((coord[0]+vec[0],coord[1]+vec[1]))+"->"+str(shitf_coord)+" (-1,1): "+str(C2x2_LD.size())) + + #----- build lower part C2x2_LD--C2x2_RD ----------------------------------- + # 0 0->3 0 3->1 + # |/23->12 |/23->45 & permute |/12->23 |/45 + # C2x2_LD--1 1--C2x2_RD C2x2_LD------C2x2_RD + # TODO is it worthy(performance-wise) to instead overwrite one of C2x2_LD,C2x2_RD ? + lower_half = contract(C2x2_LD, C2x2_RD, ([1],[1])) + lower_half = permute(lower_half, (0,3,1,2,4,5)) + + # construct reduced density matrix by contracting lower and upper halfs + # C2x2_LU------C2x2_RU + # |\23->01 |\45->23 + # 0 1 + # 0 1 + # |/23->45 |/45->67 + # C2x2_LD------C2x2_RD + rdm = contract(upper_half,lower_half,([0,1],[0,1])) + + # permute into order of s0,s1,s2,s3;s0',s1',s2',s3' where primed indices + # represent "ket" + # 01234567->02461357 + # symmetrize and normalize + rdm= contiguous(permute(rdm, (0,2,4,6,1,3,5,7))) + rdm= _sym_pos_def_rdm(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + + return rdm + + +def rdm2x2_kagome(coord, state, env, sites_to_keep_00=('A', 'B', 'C'), sites_to_keep_10=('A', 'B', 'C'), sites_to_keep_01=('A', 'B', 'C'), sites_to_keep_11=('A', 'B', 'C'), sym_pos_def=False, verbosity=0): + r""" + :param coord: vertex (x,y) specifies upper left site of 2x2 subsystem + :param state: underlying wavefunction + :param env: environment corresponding to ``state`` + :param verbosity: logging verbosity + :param sites_to_keep_00: physical sites needed for the unit cell at coord + (0, 0) + :param sites_to_keep_10: physical sites needed for the unit cell at coord + (1, 0) + :param sites_to_keep_01: physical sites needed for the unit cell at coord + (0, 1) + :param sites_to_keep_11: physical sites needed for the unit cell at coord + (1, 1) + :type coord: tuple(int,int) + :type state: IPEPS_KAGOME + :type env: ENV + :type verbosity: int + :return: 4-site reduced density matrix with indices :math:`s_0s_1s_2s_3;s'_0s'_1s'_2s'_3` + :rtype: torch.tensor + + Computes 4-site reduced density matrix :math:`\rho_{2x2}` of 2x2 subsystem specified + by the vertex ``coord`` of its upper left corner using strategy: + + 1. compute four individual corners + 2. construct upper and lower half of the network + 3. contract upper and lower half to obtain final reduced density matrix + + :: + + C--T------------------T------------------C = C2x2_LU(coord)--------C2x2(coord+(1,0)) + | | | | | | + T--A^+A(coord)--------A^+A(coord+(1,0))--T C2x2_LD(coord+(0,1))--C2x2(coord+(1,1)) + | | | | + T--A^+A(coord+(0,1))--A^+A(coord+(1,1))--T + | | | | + C--T------------------T------------------C + + The physical indices `s` and `s'` of on-sites tensors :math:`A` (and :math:`A^\dagger`) + at vertices ``coord``, ``coord+(1,0)``, ``coord+(0,1)``, and ``coord+(1,1)`` are + left uncontracted and given in the same order:: + + s0 s1 + s2 s3 + + """ + who = "rdm2x2_kagome" + # ----- building C2x2_LU ---------------------------------------------------- + C = env.C[(state.vertexToSite(coord), (-1, -1))] + T1 = env.T[(state.vertexToSite(coord), (0, -1))] + T2 = env.T[(state.vertexToSite(coord), (-1, 0))] + dimsA = state.site(coord).size() + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep_00) + a = view(contiguous(a), (dimsA[1] ** 2, dimsA[2] ** 2, dimsA[3] ** 2, dimsA[4] ** 2, rpd, rpd)) + # C--10--T1--2 + # 0 1 + C2x2_LU = contract(C, T1, ([1], [0])) + + # C------T1--2->1 + # 0 1->0 + # 0 + # T2--2->3 + # 1->2 + C2x2_LU = contract(C2x2_LU, T2, ([0], [0])) + + # C-------T1--1->0 + # | 0 + # | 0 + # T2--3 1 a--3 + # 2->1 2\45 + C2x2_LU = contract(C2x2_LU, a, ([0, 3], [0, 1])) + + # permute 012345->120345 + # reshape (12)(03)45->0123 + # C2x2--1 + # |\23 + # 0 + C2x2_LU = contiguous(permute(C2x2_LU, (1, 2, 0, 3, 4, 5))) + C2x2_LU = view(C2x2_LU, (T2.size(1) * a.size(2), T1.size(2) * a.size(3), rpd, rpd)) + if verbosity > 0: + print("C2X2 LU " + str(coord) + "->" + str(state.vertexToSite(coord)) + " (-1,-1): " + str(C2x2_LU.size())) + + # ----- building C2x2_RU ---------------------------------------------------- + vec = (1, 0) + shitf_coord = state.vertexToSite((coord[0] + vec[0], coord[1] + vec[1])) + C = env.C[(shitf_coord, (1, -1))] + T1 = env.T[(shitf_coord, (1, 0))] + T2 = env.T[(shitf_coord, (0, -1))] + dimsA = state.site(shitf_coord).size() + + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep_10) + a = view(contiguous(a), (dimsA[1] ** 2, dimsA[2] ** 2, dimsA[3] ** 2, dimsA[4] ** 2, rpd, rpd)) + # 0--C + # 1 + # 0 + # 1--T1 + # 2 + C2x2_RU = contract(C, T1, ([1], [0])) + + # 2<-0--T2--2 0--C + # 3<-1 | + # 0<-1--T1 + # 1<-2 + C2x2_RU = contract(C2x2_RU, T2, ([0], [2])) + + # 1<-2--T2------C + # 3 | + # 45\0 | + # 2<-1--a--3 0--T1 + # 3<-2 0<-1 + C2x2_RU = contract(C2x2_RU, a, ([0, 3], [3, 0])) + + # permute 012334->120345 + # reshape (12)(03)45->0123 + # 0--C2x2 + # 23/| + # 1 + C2x2_RU = contiguous(permute(C2x2_RU, (1, 2, 0, 3, 4, 5))) + C2x2_RU = view(C2x2_RU, (T2.size(0) * a.size(1), T1.size(2) * a.size(2), rpd, rpd)) + if verbosity > 0: + print("C2X2 RU " + str((coord[0] + vec[0], coord[1] + vec[1])) + "->" + str(shitf_coord) + " (1,-1): " + str( + C2x2_RU.size())) + + # ----- build upper part C2x2_LU--C2x2_RU ----------------------------------- + # C2x2_LU--1 0--C2x2_RU C2x2_LU------C2x2_RU + # |\23->12 |\23->45 & permute |\12->23 |\45 + # 0 1->3 0 3->1 + # TODO is it worthy(performance-wise) to instead overwrite one of C2x2_LU,C2x2_RU ? + upper_half = contract(C2x2_LU, C2x2_RU, ([1], [0])) + upper_half = permute(upper_half, (0, 3, 1, 2, 4, 5)) + + # ----- building C2x2_RD ---------------------------------------------------- + vec = (1, 1) + shitf_coord = state.vertexToSite((coord[0] + vec[0], coord[1] + vec[1])) + C = env.C[(shitf_coord, (1, 1))] + T1 = env.T[(shitf_coord, (0, 1))] + T2 = env.T[(shitf_coord, (1, 0))] + dimsA = state.site(shitf_coord).size() + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep_11) + a = view(contiguous(a), (dimsA[1] ** 2, dimsA[2] ** 2, dimsA[3] ** 2, dimsA[4] ** 2, rpd, rpd)) + # 1<-0 0 + # 2<-1--T1--2 1--C + C2x2_RD = contract(C, T1, ([1], [2])) + + # 2<-0 + # 3<-1--T2 + # 2 + # 0<-1 0 + # 1<-2--T1---C + C2x2_RD = contract(C2x2_RD, T2, ([0], [2])) + + # 2<-0 1<-2 + # 3<-1--a--3 3--T2 + # 2\45 | + # 0 | + # 0<-1--T1------C + C2x2_RD = contract(C2x2_RD, a, ([0, 3], [2, 3])) + + # permute 012345->120345 + # reshape (12)(03)45->0123 + C2x2_RD = contiguous(permute(C2x2_RD, (1, 2, 0, 3, 4, 5))) + C2x2_RD = view(C2x2_RD, (T2.size(0) * a.size(0), T1.size(1) * a.size(1), rpd, rpd)) + + # 0 + # |/23 + # 1--C2x2 + if verbosity > 0: + print("C2X2 RD " + str((coord[0] + vec[0], coord[1] + vec[1])) + "->" + str(shitf_coord) + " (1,1): " + str( + C2x2_RD.size())) + + # ----- building C2x2_LD ---------------------------------------------------- + vec = (0, 1) + shitf_coord = state.vertexToSite((coord[0] + vec[0], coord[1] + vec[1])) + C = env.C[(shitf_coord, (-1, 1))] + T1 = env.T[(shitf_coord, (-1, 0))] + T2 = env.T[(shitf_coord, (0, 1))] + dimsA = state.site(shitf_coord).size() + a, rpd = build_reduced_density_matrix_kagome(coord, state, site_types=sites_to_keep_01) + a = view(contiguous(a), (dimsA[1] ** 2, dimsA[2] ** 2, dimsA[3] ** 2, dimsA[4] ** 2, rpd, rpd)) + # 0->1 + # T1--2 + # 1 + # 0 + # C--1->0 + C2x2_LD = contract(C, T1, ([0], [1])) + + # 1->0 + # T1--2->1 + # | + # | 0->2 + # C--0 1--T2--2->3 + C2x2_LD = contract(C2x2_LD, T2, ([0], [1])) + + # 0 0->2 + # T1--1 1--a--3 + # | 2\45 + # | 2 + # C--------T2--3->1 + C2x2_LD = contract(C2x2_LD, a, ([1, 2], [1, 2])) + + # permute 012345->021345 + # reshape (02)(13)45->0123 + # 0 + # |/23 + # C2x2--1 + C2x2_LD = contiguous(permute(C2x2_LD, (0, 2, 1, 3, 4, 5))) + C2x2_LD = view(C2x2_LD, (T1.size(0) * a.size(0), T2.size(2) * a.size(3), rpd, rpd)) + if verbosity > 0: + print("C2X2 LD " + str((coord[0] + vec[0], coord[1] + vec[1])) + "->" + str(shitf_coord) + " (-1,1): " + str( + C2x2_LD.size())) + + # ----- build lower part C2x2_LD--C2x2_RD ----------------------------------- + # 0 0->3 0 3->1 + # |/23->12 |/23->45 & permute |/12->23 |/45 + # C2x2_LD--1 1--C2x2_RD C2x2_LD------C2x2_RD + # TODO is it worthy(performance-wise) to instead overwrite one of C2x2_LD,C2x2_RD ? + lower_half = contract(C2x2_LD, C2x2_RD, ([1], [1])) + lower_half = permute(lower_half, (0, 3, 1, 2, 4, 5)) + + # construct reduced density matrix by contracting lower and upper halfs + # C2x2_LU------C2x2_RU + # |\23->01 |\45->23 + # 0 1 + # 0 1 + # |/23->45 |/45->67 + # C2x2_LD------C2x2_RD + rdm = contract(upper_half, lower_half, ([0, 1], [0, 1])) + + # permute into order of s0,s1,s2,s3;s0',s1',s2',s3' where primed indices + # represent "ket" + # 01234567->02461357 + # symmetrize and normalize + rdm = contiguous(permute(rdm, (0, 2, 4, 6, 1, 3, 5, 7))) + rdm = _sym_pos_def_rdm(rdm, sym_pos_def=sym_pos_def, verbosity=verbosity, who=who) + + return rdm diff --git a/examples/ctmrg_kagome.py b/examples/ctmrg_kagome.py new file mode 100644 index 00000000..c38dca23 --- /dev/null +++ b/examples/ctmrg_kagome.py @@ -0,0 +1,291 @@ +import os +os.environ['MKL_THREADING_LAYER'] = 'GNU' +import context +import torch +import argparse +import config as cfg +from ipeps.ipeps_kagome import * +from ctm.generic.env import * +from ctm.generic import ctmrg +# from ctm.generic import transferops +from models import kagome +import unittest +import numpy as np +import pickle + +# parse command line args and build necessary configuration objects +parser= cfg.get_args_parser() +# additional model-dependent arguments +parser.add_argument("--theta", type=float, default=0.5, help="parametrization between 2- and 3-site terms. theta * pi") +parser.add_argument("--phi", type=float, default=0., help="parametrization between normal and chiral terms. phi * pi") +parser.add_argument("--tiling", default="1SITE", help="tiling of the lattice") +# additional observables-related arguments +parser.add_argument("--corrf_r", type=int, default=1, help="maximal correlation function distance") +parser.add_argument("--top_n", type=int, default=8, help="number of leading eigenvalues"+ + "of transfer operator to compute") +parser.add_argument("--input_prefix", type=str, default="theta_0_phi_0_bonddim_3_chi_16", help="parameters of input state") +parser.add_argument("--output_path", type=str, default="/scratch/yx51/kagome", help="path of output") +parser.add_argument("--restrictions", type=bool, default=False, help="restrictions on the 5 site tensors") +args, unknown_args = parser.parse_known_args() + + +def main(): + cfg.configure(args) + cfg.print_config() + torch.set_num_threads(args.omp_cores) + torch.manual_seed(args.seed) + + param_j = np.round(np.cos(np.pi*args.theta), decimals=12) + param_k = np.round(np.sin(np.pi*args.theta) * np.cos(np.pi*args.phi), decimals=12) + param_h = np.round(np.sin(np.pi*args.theta) * np.sin(np.pi*args.phi), decimals=12) + print("J = {}; K = {}; H = {}".format(param_j, param_k, param_h)) + model = kagome.KAGOME(phys_dim=3, j=param_j, k=param_k, h=param_h) + # initialize an ipeps + # 1) define lattice-tiling function, that maps arbitrary vertex of square lattice + # coord into one of coordinates within unit-cell of iPEPS ansatz + if args.tiling == "BIPARTITE": + def lattice_to_site(coord): + vx = (coord[0] + abs(coord[0]) * 2) % 2 + vy = abs(coord[1]) + return ((vx + vy) % 2, 0) + elif args.tiling == "1SITE": + def lattice_to_site(coord): + return (0, 0) + elif args.tiling == "2SITE": + def lattice_to_site(coord): + vx = (coord[0] + abs(coord[0]) * 2) % 2 + # vy = (coord[1] + abs(coord[1]) * 1) % 1 + return (vx, 0) + elif args.tiling == "4SITE": + def lattice_to_site(coord): + vx = (coord[0] + abs(coord[0]) * 2) % 2 + vy = (coord[1] + abs(coord[1]) * 2) % 2 + return (vx, vy) + # elif args.tiling == "8SITE": + # def lattice_to_site(coord): + # shift_x = coord[0] + 2 * (coord[1] // 2) + # vx = shift_x % 4 + # vy = coord[1] % 2 + # return (vx, vy) + else: + raise ValueError("Invalid tiling: " + str(args.tiling) + " Supported options: " \ + + "1SITE, BIPARTITE, 2SITE, 4SITE, 8SITE") + + # initialize an ipeps + if args.instate != None: + state = read_ipeps_kagome(args.instate, restrictions=args.restrictions, vertexToSite=lattice_to_site) + if args.bond_dim > max(state.get_aux_bond_dims()): + # extend the auxiliary dimensions + state = extend_bond_dim_kagome(state, args.bond_dim) + state.add_noise(args.instate_noise) + elif args.ipeps_init_type == 'RANDOM': + bond_dim = args.bond_dim + A = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + B = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + C = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + RD = torch.rand((bond_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + RU = torch.rand((bond_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + + # normalization of initial random tensors + A = A/torch.max(torch.abs(A)) + B = B/torch.max(torch.abs(B)) + C = C/torch.max(torch.abs(C)) + RD = RD/torch.max(torch.abs(RD)) + RU = RU/torch.max(torch.abs(RU)) + kagome_sites = {(0, 0, 0): A, (0, 0, 1): B, (0, 0, 2): C, (0, 0, 3): RD, (0, 0, 4): RU} + if args.tiling in ["BIPARTITE", "2SITE", "4SITE"]: + A2 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + B2 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + C2 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + RD2 = torch.rand((bond_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + RU2 = torch.rand((bond_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + kagome_sites[(1, 0, 0)] = A2/torch.max(torch.abs(A2)) + kagome_sites[(1, 0, 1)] = B2/torch.max(torch.abs(B2)) + kagome_sites[(1, 0, 2)] = C2/torch.max(torch.abs(C2)) + kagome_sites[(1, 0, 3)] = RD2/torch.max(torch.abs(RD2)) + kagome_sites[(1, 0, 4)] = RU2/torch.max(torch.abs(RU2)) + if args.tiling in ["4SITE"]: + A3 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + B3 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + C3 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + RD3 = torch.rand((bond_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + RU3 = torch.rand((bond_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + kagome_sites[(0, 1, 0)] = A3 / torch.max(torch.abs(A3)) + kagome_sites[(0, 1, 1)] = B3 / torch.max(torch.abs(B3)) + kagome_sites[(0, 1, 2)] = C3 / torch.max(torch.abs(C3)) + kagome_sites[(0, 1, 3)] = RD3 / torch.max(torch.abs(RD3)) + kagome_sites[(0, 1, 4)] = RU3 / torch.max(torch.abs(RU3)) + + A4 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + B4 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + C4 = torch.rand((model.phys_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + RD4 = torch.rand((bond_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + RU4 = torch.rand((bond_dim, bond_dim, bond_dim), dtype=cfg.global_args.torch_dtype, + device=cfg.global_args.device) + kagome_sites[(1, 1, 0)] = A4 / torch.max(torch.abs(A4)) + kagome_sites[(1, 1, 1)] = B4 / torch.max(torch.abs(B4)) + kagome_sites[(1, 1, 2)] = C4 / torch.max(torch.abs(C4)) + kagome_sites[(1, 1, 3)] = RD4 / torch.max(torch.abs(RD4)) + kagome_sites[(1, 1, 4)] = RU4 / torch.max(torch.abs(RU4)) + + state = IPEPS_KAGOME(kagome_sites, restrictions=args.restrictions, vertexToSite=lattice_to_site) + else: + raise ValueError("Missing trial state: -instate=None and -ipeps_init_type= " \ + + str(args.ipeps_init_type) + " is not supported") + + if not state.dtype==model.dtype: + cfg.global_args.dtype= state.dtype + print(f"dtype of initial state {state.dtype} and model {model.dtype} do not match.") + print(f"Setting default dtype to {cfg.global_args.dtype} and reinitializing "\ + +" the model") + model= kagome.KAGOME(phys_dim=3, j=param_j, k=param_k, h=param_h) + + print(state) + + # 2) select the "energy" function + if args.tiling == "BIPARTITE" or args.tiling == "2SITE": + energy_f = model.energy_1site + eval_obs_f= model.eval_obs + elif args.tiling == "1SITE": + energy_f= model.energy_1site + # TODO include eval_obs with rotation on B-sublattice + eval_obs_f= model.eval_obs + elif args.tiling == "4SITE": + energy_f = model.energy_1site + eval_obs_f= model.eval_obs + else: + raise ValueError("Invalid tiling: "+str(args.tiling)+" Supported options: "\ + +"BIPARTITE, 2SITE, 4SITE") + + def ctmrg_conv_energy(state, env, history, ctm_args=cfg.ctm_args): + with torch.no_grad(): + if not history: + history = [] + e_curr = energy_f(state, env) + obs_values, obs_labels = model.eval_obs(state, env) + history.append([e_curr.item()] + obs_values) + print(", ".join([f"{len(history)}", f"{e_curr}"] + [f"{v}" for v in obs_values])) + + if len(history) > 1 and abs(history[-1][0] - history[-2][0]) < ctm_args.ctm_conv_tol: + return True, history + return False, history + + ctm_env_init = ENV(args.chi, state) + init_env(state, ctm_env_init) + print(ctm_env_init) + + e_curr0 = energy_f(state, ctm_env_init) + obs_values0, obs_labels = model.eval_obs(state, ctm_env_init) + + print(", ".join(["epoch", "energy"] + obs_labels)) + print(", ".join([f"{-1}", f"{e_curr0}"] + [f"{v}" for v in obs_values0])) + + ctm_env_init, *ctm_log = ctmrg.run(state, ctm_env_init, conv_check=ctmrg_conv_energy) + + e_final = energy_f(state, ctm_env_init) + obs_values, obs_labels = model.eval_obs(state, ctm_env_init) + print(", ".join(["epoch", "energy"] + obs_labels)) + print(", ".join([f"{-1}", f"{e_final}"] + [f"{v}" for v in obs_values])) + energy_dn, energy_up = model.energy_1triangle(state, ctm_env_init) + + obs = dict() + obs["energy"] = e_final + obs["energy_dn"] = energy_dn + obs["energy_up"] = energy_up + generators = model.eval_generators(state, ctm_env_init) + c1 = model.eval_C1(state, ctm_env_init) + for label, value in generators.items(): + obs[label] = value + for label, value in c1.items(): + obs[label] = value + c2s = model.eval_C2(state, ctm_env_init) + for label, value in c2s.items(): + obs[label] = value + + def inner_prod(gens1, gens2): + dot_prod = gens1[0] * gens2[0] + 0.75 * gens1[7] * gens2[7] + 0.5 * (gens1[1] * gens2[2] + + gens1[3] * gens2[4] + gens1[5] * gens2[6] + + gens1[1] * gens2[2] + gens1[3] * gens2[4] + + gens1[5] * gens2[6]) + return dot_prod + + obs["f_spin_A"] = inner_prod(generators["generators_A"], generators["generators_A"]) + obs["f_spin_B"] = inner_prod(generators["generators_B"], generators["generators_B"]) + obs["f_spin_C"] = inner_prod(generators["generators_C"], generators["generators_C"]) + obs["avg_bonds_dn"] = obs_values[0] + obs["avg_bonds_up"] = obs_values[1] + obs["chirality_dn"] = obs_values[2] + obs["chirality_up"] = obs_values[3] + print(obs) + # if args.restrictions: + # filename_obs = "./data/onsite_obs_restricted_theta_{}_phi_{}_bonddim_{}_chi_{}.json".format(int(args.theta*100), int(args.phi*100), args.bond_dim, args.chi) + # else: + # filename_obs = "./data/onsite_obs_theta_{}_phi_{}_bonddim_{}_chi_{}.json".format(int(args.theta * 100), int(args.phi * 100), args.bond_dim, args.chi) + # filename_obs = "{}/obs_{}.pkl".format(args.output_path, args.input_prefix) + # with open(filename_obs, "wb") as fp: + # pickle.dump(obs, fp) + # + # with open(filename_obs, "rb") as fp: + # tmp = pickle.load(fp) + # print(tmp) + + # corrSS = model.eval_corrf_SS((0, 0), (1, 0), state, ctm_env_init, args.corrf_r) + # print("\n\nSS[(0,0),(1,0)] r " + " ".join([label for label in corrSS.keys()])) + # for i in range(args.corrf_r): + # print(f"{i} " + " ".join([f"{corrSS[label][i]}" for label in corrSS.keys()])) + # + # corrSS = model.eval_corrf_SS((0, 0), (0, 1), state, ctm_env_init, args.corrf_r) + # print("\n\nSS[(0,0),(0,1)] r " + " ".join([label for label in corrSS.keys()])) + # for i in range(args.corrf_r): + # print(f"{i} " + " ".join([f"{corrSS[label][i]}" for label in corrSS.keys()])) + # + # corrSSSS = model.eval_corrf_SSSS((0, 0), (1, 0), state, ctm_env_init, args.corrf_r) + # print("\n\nSSSS[(0,0),(1,0)] r " + " ".join([label for label in corrSSSS.keys()])) + # for i in range(args.corrf_r): + # print(f"{i} " + " ".join([f"{corrSSSS[label][i]}" for label in corrSSSS.keys()])) + # + # # environment diagnostics + # print("\n") + # for c_loc, c_ten in ctm_env_init.C.items(): + # u, s, v = torch.svd(c_ten, compute_uv=False) + # print(f"spectrum C[{c_loc}]") + # for i in range(args.chi): + # print(f"{i} {s[i]}") + # + # transfer operator spectrum + # site_dir_list = [((0, 0), (1, 0)), ((0, 0), (0, 1)), ((1, 1), (1, 0)), ((1, 1), (0, 1))] + # for sdp in site_dir_list: + # print(f"\n\nspectrum(T)[{sdp[0]},{sdp[1]}]") + # l = transferops.get_Top_spec(args.top_n, *sdp, state, ctm_env_init) + # for i in range(l.size()[0]): + # print(f"{i} {l[i, 0]} {l[i, 1]}") + # + # # environment diagnostics + # for c_loc,c_ten in ctm_env_init.C.items(): + # u,s,v= torch.svd(c_ten, compute_uv=False) + # print(f"\n\nspectrum C[{c_loc}]") + # for i in range(args.chi): + # print(f"{i} {s[i]}") + + +if __name__ == '__main__': + if len(unknown_args) > 0: + print("args not recognized: " + str(unknown_args)) + raise Exception("Unknown command line arguments") + main() + + diff --git a/examples/optim_kagome.py b/examples/optim_kagome.py new file mode 100644 index 00000000..5215eff9 --- /dev/null +++ b/examples/optim_kagome.py @@ -0,0 +1,216 @@ +import os +os.environ['MKL_THREADING_LAYER'] = 'GNU' +import context +import torch +import argparse +import config as cfg +from ipeps.ipeps import * +from ctm.generic.env import * +from ctm.generic import ctmrg +# from ctm.generic import ctmrg_sl as ctmrg +from models import bbh_j1j2 + +from optim.ad_optim_lbfgs_mod import optimize_state +import unittest +import logging + +log = logging.getLogger(__name__) + +# parse command line args and build necessary configuration objects +parser = cfg.get_args_parser() +# additional model-dependent arguments +parser.add_argument("--theta", type=float, default=0., help="theta") +parser.add_argument("--ratio", type=float, default=1., help="y/x ratio") +parser.add_argument("--j1", type=float, default=1., help="nn bilinear coupling") +parser.add_argument("--k1", type=float, default=0., help="nn biquadratic coupling") +parser.add_argument("--j2", type=float, default=0., help="nnn bilinear coupling") +parser.add_argument("--k2", type=float, default=0., help="nnn biquadratic coupling") +parser.add_argument("--tiling", default="BIPARTITE", help="tiling of the lattice") +args, unknown_args = parser.parse_known_args() + + +def main(): + cfg.configure(args) + cfg.print_config() + torch.set_num_threads(args.omp_cores) + torch.manual_seed(args.seed) + + model = bbh_j1j2.BBH_J1J2(j1=args.j1, k1=args.k1, j2=args.j2, k2=args.k2) + + # initialize an ipeps + # 1) define lattice-tiling function, that maps arbitrary vertex of square lattice + # coord into one of coordinates within unit-cell of iPEPS ansatz + if args.tiling == "BIPARTITE": + def lattice_to_site(coord): + vx = (coord[0] + abs(coord[0]) * 2) % 2 + vy = abs(coord[1]) + return ((vx + vy) % 2, 0) + elif args.tiling == "2SITE": + def lattice_to_site(coord): + vx = (coord[0] + abs(coord[0]) * 2) % 2 + # vy = (coord[1] + abs(coord[1]) * 1) % 1 + return (vx, 0) + elif args.tiling == "4SITE": + def lattice_to_site(coord): + vx = (coord[0] + abs(coord[0]) * 2) % 2 + vy = (coord[1] + abs(coord[1]) * 2) % 2 + return (vx, vy) + # elif args.tiling == "8SITE": + # def lattice_to_site(coord): + # shift_x = coord[0] + 2 * (coord[1] // 2) + # vx = shift_x % 4 + # vy = coord[1] % 2 + # return (vx, vy) + else: + raise ValueError("Invalid tiling: " + str(args.tiling) + " Supported options: " \ + + "BIPARTITE, 2SITE, 4SITE") + + if args.instate != None: + state = read_ipeps(args.instate, vertexToSite=lattice_to_site) + if args.bond_dim > max(state.get_aux_bond_dims()): + # extend the auxiliary dimensions + state = extend_bond_dim(state, args.bond_dim) + state.add_noise(args.instate_noise) + elif args.opt_resume is not None: + if args.tiling == "BIPARTITE" or args.tiling == "2SITE": + state = IPEPS(dict(), lX=2, lY=1) + elif args.tiling == "4SITE": + state = IPEPS(dict(), lX=2, lY=2) + # elif args.tiling == "8SITE": + # state = IPEPS(dict(), lX=4, lY=2) + state.load_checkpoint(args.opt_resume) + elif args.ipeps_init_type == 'RANDOM': + bond_dim = args.bond_dim + + A = torch.rand((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + B = torch.rand((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + # A = torch.zeros((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + # dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + # A[:, 0, :, 0, :] = torch.rand_like(A[:, 0, :, 0, :], dtype=A.torch_dtype, device=A.device) + # A = (1 - args.ratio) * A + args.ratio * torch.rand_like(A, dtype=A.torch_dtype, device=A.device) + # + # B = torch.zeros((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), + # dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + # B[:, 0, :, 0, :] = torch.rand_like(B[:, 0, :, 0, :], dtype=B.torch_dtype, device=B.device) + # B = (1 - args.ratio) * B + args.ratio * torch.rand_like(B, dtype=B.torch_dtype, device=B.device) + + # normalization of initial random tensors + A = A / torch.max(torch.abs(A)) + B = B / torch.max(torch.abs(B)) + sites = {(0, 0): A, (1, 0): B} + if args.tiling == "4SITE" or args.tiling == "8SITE": + C = torch.rand((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + D = torch.rand((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + dtype=cfg.global_args.torch_dtype, device=cfg.global_args.device) + sites[(0, 1)] = C / torch.max(torch.abs(C)) + sites[(1, 1)] = D / torch.max(torch.abs(D)) + + # if args.tiling == "8SITE": + # E = torch.rand((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + # dtype=cfg.global_args.dtype, device=cfg.global_args.device) + # F = torch.rand((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + # dtype=cfg.global_args.dtype, device=cfg.global_args.device) + # G = torch.rand((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + # dtype=cfg.global_args.dtype, device=cfg.global_args.device) + # H = torch.rand((model.phys_dim, bond_dim, bond_dim, bond_dim, bond_dim), \ + # dtype=cfg.global_args.dtype, device=cfg.global_args.device) + # sites[(2, 0)] = E / torch.max(torch.abs(E)) + # sites[(3, 0)] = F / torch.max(torch.abs(F)) + # sites[(2, 1)] = G / torch.max(torch.abs(G)) + # sites[(3, 1)] = H / torch.max(torch.abs(H)) + state = IPEPS(sites, vertexToSite=lattice_to_site) + else: + raise ValueError("Missing trial state: --instate=None and --ipeps_init_type= " \ + + str(args.ipeps_init_type) + " is not supported") + + if not state.dtype==model.dtype: + cfg.global_args.dtype= state.dtype + print(f"dtype of initial state {state.dtype} and model {model.dtype} do not match.") + print(f"Setting default dtype to {cfg.global_args.dtype} and reinitializing "\ + +" the model") + model= bbh_j1j2.BBH_J1J2(j1=args.j1, k1=args.k1, j2=args.j2, k2=args.k2) + + print(state) + + # 2) select the "energy" function + if args.tiling == "BIPARTITE" or args.tiling == "2SITE": + energy_f = model.energy_2x2 + elif args.tiling == "4SITE": + energy_f = model.energy_2x2 + # elif args.tiling == "8SITE": + # energy_f = model.energy_2x2_8site + else: + raise ValueError("Invalid tiling: " + str(args.tiling) + " Supported options: " \ + + "BIPARTITE, 2SITE, 4SITE, 8SITE") + + @torch.no_grad() + def ctmrg_conv_energy(state, env, history, ctm_args=cfg.ctm_args): + if not history: + history = [] + e_curr = energy_f(state, env) + history.append(e_curr.item()) + + if (len(history) > 1 and abs(history[-1] - history[-2]) < ctm_args.ctm_conv_tol) \ + or len(history) >= ctm_args.ctm_max_iter: + log.info({"history_length": len(history), "history": history}) + return True, history + return False, history + + ctm_env = ENV(args.chi, state) + init_env(state, ctm_env) + + ctm_env, *ctm_log = ctmrg.run(state, ctm_env, conv_check=ctmrg_conv_energy) + loss0 = energy_f(state, ctm_env) + obs_values, obs_labels = model.eval_obs(state, ctm_env) + print(", ".join(["epoch", "energy"] + obs_labels)) + print(", ".join([f"{-1}", f"{loss0}"] + [f"{v}" for v in obs_values])) + + def loss_fn(state, ctm_env_in, opt_context): + ctm_args = opt_context["ctm_args"] + opt_args = opt_context["opt_args"] + + # possibly re-initialize the environment + if opt_args.opt_ctm_reinit: + init_env(state, ctm_env_in) + + # 1) compute environment by CTMRG + ctm_env_out, *ctm_log = ctmrg.run(state, ctm_env_in, \ + conv_check=ctmrg_conv_energy, ctm_args=ctm_args) + + # 2) evaluate loss with the converged environment + loss = energy_f(state, ctm_env_out) + + return (loss, ctm_env_out, *ctm_log) + + @torch.no_grad() + def obs_fn(state, ctm_env, opt_context): + if ("line_search" in opt_context.keys() and not opt_context["line_search"]) \ + or not "line_search" in opt_context.keys(): + epoch = len(opt_context["loss_history"]["loss"]) + loss = opt_context["loss_history"]["loss"][-1] + obs_values, obs_labels = model.eval_obs(state, ctm_env) + print(", ".join([f"{epoch}", f"{loss}"] + [f"{v}" for v in obs_values])) + log.info("Norm(sites): " + ", ".join([f"{t.norm()}" for c, t in state.sites.items()])) + + # optimize + optimize_state(state, ctm_env, loss_fn, obs_fn=obs_fn) + + # compute final observables for the best variational state + outputstatefile = args.out_prefix + "_state.json" + state = read_ipeps(outputstatefile, vertexToSite=state.vertexToSite) + ctm_env = ENV(args.chi, state) + init_env(state, ctm_env) + ctm_env, *ctm_log = ctmrg.run(state, ctm_env, conv_check=ctmrg_conv_energy) + opt_energy = energy_f(state, ctm_env) + obs_values, obs_labels = model.eval_obs(state, ctm_env) + print(", ".join([f"{args.opt_max_iter}", f"{opt_energy}"] + [f"{v}" for v in obs_values])) + + +if __name__ == '__main__': + if len(unknown_args) > 0: + print("args not recognized: " + str(unknown_args)) + raise Exception("Unknown command line arguments") + main() diff --git a/groups/su3.py b/groups/su3.py new file mode 100644 index 00000000..fe8a2747 --- /dev/null +++ b/groups/su3.py @@ -0,0 +1,201 @@ +import torch +from math import factorial, sqrt +from tn_interface import einsum +import numpy as np + +class SU3_DEFINING(): + def __init__(self, p=1, q=0, dtype=torch.complex128, device='cpu'): + r""" + :param (p, q): labels of highest weight for su(3) representations. (1, 0) - defining representation + :param dtype: data type of matrix representation of operators + :param device: device on which the torch.tensor objects are stored + :type J: int + :type dtype: torch.dtype + :type device: int + + Build the defining representation "3" of su(3) Lie algebra using the Cartan-Weyl basis, $\lambda_i$. + + The quadratic Casimir operator of su(3) can be expressed in terms of the C-W basis, defined as follow. + + .. math:: + \begin{align*} + T^\pm &= \frac{1}{2} (\lambda_1 \pm i\lambda_2) = (F_1 \pm iF_2)\\ + T^z &= \frac{1}{2} \lambda_3 = F_3\\ + V^\pm &= \frac{1}{2} (\lambda_4 \pm i\lambda_5) = (F_4 \pm iF_5)\\ + U^\pm &= \frac{1}{2} (\lambda_6 \pm i\lambda_7) = (F_6 \pm iF_7)\\ + Y &= \frac{1}{\sqrt{3}} \lambda_8 = \frac{2}{\sqrt{3}} F_8 + \end{align*} + + \begin{align*} + C_1 = \sum_{k}{F_k F_k} &= \frac{1}{2} (T^+ T^- + T^- T^+ + V^+ V^- + V^- V^+ + U^+ U^- + U^- U^+) \\ + &+ T^z T^z + \frac{3}{4} Y Y + \end{align*} + """ + self.p = p + self.q = q + self.dtype = dtype + self.device = device + + def I(self): + r""" + :return: Identity operator of irrep + :rtype: torch.tensor + """ + return get_op("I", dtype=self.dtype, device=self.device) + + def TZ(self): + r""" + :return: :math:`T^z` operator of irrep + :rtype: torch.tensor + """ + return get_op("tz", dtype=self.dtype, device=self.device) + + def Y(self): + r""" + :return: :math:`Y` operator of irrep + :rtype: torch.tensor + """ + return get_op("y", dtype=self.dtype, device=self.device) + + def TP(self): + r""" + :return: :math:`T^+` operator of irrep + :rtype: torch.tensor + """ + return get_op("tp", dtype=self.dtype, device=self.device) + + def TM(self): + r""" + :return: :math:`T^-` operator of irrep + :rtype: torch.tensor + """ + return get_op("tm", dtype=self.dtype, device=self.device) + + def VP(self): + r""" + :return: :math:`V^+` operator of irrep + :rtype: torch.tensor + """ + return get_op("vp", dtype=self.dtype, device=self.device) + + def VM(self): + r""" + :return: :math:`V^-` operator of irrep + :rtype: torch.tensor + """ + return get_op("vm", dtype=self.dtype, device=self.device) + + def UP(self): + r""" + :return: :math:`U^+` operator of irrep + :rtype: torch.tensor + """ + return get_op("up", dtype=self.dtype, device=self.device) + + def UM(self): + r""" + :return: :math:`U^-` operator of irrep + :rtype: torch.tensor + """ + return get_op("um", dtype=self.dtype, device=self.device) + + def C1(self): + r""" + :return: The quadratic Casimir of su(3) as rank-4 for tensor + :rtype: torch.tensor + """ + expr_kron = 'ij,ab->iajb' + # spin-spin interaction \sum_k{\vec{F}_{1,k}\vec{S}_{2,k}} between F-spins on sites 1 and 2 + C1 = einsum(expr_kron, self.TZ(), self.TZ()) + 0.75 * einsum(expr_kron, self.Y(), self.Y())\ + + 0.5 * (einsum(expr_kron, self.TP(), self.TM()) + einsum(expr_kron, self.TM(), self.TP()) + + einsum(expr_kron, self.VP(), self.VM()) + einsum(expr_kron, self.VM(), self.VP()) + + einsum(expr_kron, self.UP(), self.UM()) + einsum(expr_kron, self.UM(), self.UP())) + return C1 + + def C2(self): + r""" + :return: The cubic Casimir of su(3) as rank-6 for tensor + :rtype: torch.tensor + """ + expr_kron = 'ia,jb,kc->ijkabc' + Fs = dict() + Fs["f1"] = 0.5 * (self.TP() + self.TM()) + Fs["f2"] = - 0.5j * (self.TP() - self.TM()) + Fs["f3"] = self.TZ() + Fs["f4"] = 0.5 * (self.VP() + self.VM()) + Fs["f5"] = - 0.5j * (self.VP() - self.VM()) + Fs["f6"] = 0.5 * (self.UP() + self.UM()) + Fs["f7"] = - 0.5j * (self.UP() - self.UM()) + Fs["f8"] = np.sqrt(3.0) / 2 * self.Y() + C2 = torch.zeros((3, 3, 3, 3, 3, 3), dtype=torch.complex128, device='cpu') + # C2 = None + for i in range(8): + for j in range(8): + for k in range(8): + d = 2 * torch.trace((Fs[f"f{i+1}"]@Fs[f"f{j+1}"]+Fs[f"f{j+1}"]@Fs[f"f{i+1}"])@Fs[f"f{k+1}"]) + C2 += d * einsum(expr_kron, Fs[f"f{i+1}"], Fs[f"f{j+1}"], Fs[f"f{k+1}"]) + + return C2 + + +def get_op(op, dtype=torch.complex128, device='cpu', dbg=False): + if op == "I": + if dbg: + print(">>>>> Constructing 1sO: Id <<<<<") + return torch.eye(3, dtype=dtype, device=device) + + elif op == "tz": + if dbg: + print(">>>>> Constructing 1sO: T^z <<<<<") + res = torch.zeros((3, 3), dtype=dtype, device=device) + res[0, 0] = 0.5 + res[1, 1] = -0.5 + return res + elif op == "y": + if dbg: + print(">>>>> Constructing 1sO: Y <<<<<") + res = torch.zeros((3, 3), dtype=dtype, device=device) + res[0, 0] = 1.0 / 3.0 + res[1, 1] = 1.0 / 3.0 + res[2, 2] = - 2.0 / 3.0 + return res + elif op == "tp": + if dbg: + print(">>>>> Constructing 1sO: T^+ <<<<<") + res = torch.zeros((3, 3), dtype=dtype, device=device) + res[0, 1] = 1.0 + return res + elif op == "tm": + if dbg: + print(">>>>> Constructing 1sO: T^- <<<<<") + res = torch.zeros((3, 3), dtype=dtype, device=device) + res[1, 0] = 1.0 + return res + elif op == "vp": + if dbg: + print(">>>>> Constructing 1sO: V^+ <<<<<") + res = torch.zeros((3, 3), dtype=dtype, device=device) + res[0, 2] = 1.0 + return res + elif op == "vm": + if dbg: + print(">>>>> Constructing 1sO: V^- <<<<<") + res = torch.zeros((3, 3), dtype=dtype, device=device) + res[2, 0] = 1.0 + return res + elif op == "up": + if dbg: + print(">>>>> Constructing 1sO: U^+ <<<<<") + res = torch.zeros((3, 3), dtype=dtype, device=device) + res[1, 2] = 1.0 + return res + elif op == "um": + if dbg: + print(">>>>> Constructing 1sO: U^- <<<<<") + res = torch.zeros((3, 3), dtype=dtype, device=device) + res[2, 1] = 1.0 + return res + else: + raise Exception("Unsupported operator requested: " + op) + +#TODO: CG series of su(3), i.e., the expansion of the tensor product of two irrep into direct sum of irreps diff --git a/ipeps/ipeps_kagome.py b/ipeps/ipeps_kagome.py new file mode 100644 index 00000000..07886486 --- /dev/null +++ b/ipeps/ipeps_kagome.py @@ -0,0 +1,361 @@ +import torch +from collections import OrderedDict +import json +import math +import config as cfg +import ipeps.ipeps as ipeps +from ipeps.tensor_io import * + +class IPEPS_KAGOME(ipeps.IPEPS): + def __init__(self, kagome_tensors, restrictions=False, vertexToSite=None, lX=None, lY=None, \ + peps_args=cfg.peps_args, global_args=cfg.global_args): + r""" + :param kagome_tensors: list of selected KAGOME-type tensors A, B, C, R1, R2; A', B', C', R1', R2'; ... + :param vertexToSite: function mapping arbitrary vertex of a square lattice + into a vertex within elementary unit cell + :param lX: length of the elementary unit cell in X direction + :param lY: length of the elementary unit cell in Y direction + :param peps_args: ipeps configuration + :param global_args: global configuration + :type kagome_tensors: list[tuple(dict(str,str,str), torch.tensor)] + :type vertexToSite: function(tuple(int,int))->tuple(int,int) + :type lX: int + :type lY: int + :type peps_args: PEPSARGS + :type global_args: GLOBALARGS + + Member ``sites`` is a dictionary of non-equivalent on-site tensors + indexed by tuple of coordinates (x,y) within the elementary unit cell. + The index-position convetion for on-site tensors is defined as follows:: + + u s + |/ + l--A--r <=> A[s,u,l,d,r] + | + d + + where s denotes physical index, and u,l,d,r label four principal directions + up, left, down, right in anti-clockwise order starting from up + + Member ``vertexToSite`` is a mapping function from vertex on a square lattice + passed in as tuple(x,y) to a corresponding tuple(x,y) within elementary unit cell. + + On-site tensor of an IPEPS object ``wfc`` at vertex (x,y) is conveniently accessed + through the member function ``site``, which internally uses ``vertexToSite`` mapping:: + + coord= (0,0) + A_00= wfc.site(coord) + + By combining the appropriate ``vertexToSite`` mapping function with elementary unit + cell specified through ``sites`` various tilings of a square lattice can be achieved:: + + # Example 1: 1-site translational iPEPS + + sites={(0,0): A} + def vertexToSite(coord): + return (0,0) + wfc= IPEPS(sites,vertexToSite) + + # resulting tiling: + # y\x -2 -1 0 1 2 + # -2 A A A A A + # -1 A A A A A + # 0 A A A A A + # 1 A A A A A + + # Example 2: 2-site bipartite iPEPS + + sites={(0,0): A, (1,0): B} + def vertexToSite(coord): + x = (coord[0] + abs(coord[0]) * 2) % 2 + y = abs(coord[1]) + return ((x + y) % 2, 0) + wfc= IPEPS(sites,vertexToSite) + + # resulting tiling: + # y\x -2 -1 0 1 2 + # -2 A B A B A + # -1 B A B A B + # 0 A B A B A + # 1 B A B A B + + + where in the last example we used default setting for ``vertexToSite``, which + maps square lattice into elementary unit cell of size `lX` x `lY` assuming + periodic boundary conditions (PBC) along both X and Y directions. + + Kagome sites: + coord_kagome - ( 0, 0, x=0(A),1(B),2(C),3(RD),4(RU) ) + + """ + self.lX = lX + self.lY = lY + if vertexToSite is not None: + self.vertexToSite = vertexToSite + else: + def vertexToSite(coord): + x = coord[0] + y = coord[1] + return ((x + abs(x)*self.lX) % self.lX, (y + abs(y)*self.lY) % self.lY) + self.vertexToSite = vertexToSite + + self.restrictions = restrictions + self.kagome_tensors = kagome_tensors + self.phys_dim = self.get_physical_dim() + self.bond_dims = self.get_aux_bond_dims() + sites = self.build_unit_cell_tensors_kagome() + + super().__init__(sites, vertexToSite=vertexToSite, peps_args=peps_args,\ + global_args=global_args) + + def __str__(self): + print(f"lX x lY: {self.lX} x {self.lY}") + for nid, coord_kagome, kagome_site in [(t[0], *t[1]) for t in enumerate(self.kagome_tensors.items())]: + print(f"A{nid} {coord_kagome}: {kagome_site.size()}") + + # show the unit cell + unit_cell = """-- A^a_kl B^b_mn -- +-- \\ / + R1_lno + | + C^c_op + | + R2_pqr + / \\ +""" + print(unit_cell) + + return "" + + def get_parameters(self): + return self.kagome_tensors.values() + + def get_checkpoint(self): + return self.kagome_tensors + + def load_checkpoint(self, checkpoint_file): + checkpoint = torch.load(checkpoint_file) + self.kagome_tensors = checkpoint["parameters"] + for kagome_t in self.kagome_tensors.values(): + kagome_t.requires_grad_(False) + self.sites = self.build_unit_cell_tensors_kagome() + + def kagome_vertex_to_vertex(self, coord_kagome): + x = coord_kagome[0] + y = coord_kagome[1] + # z = coord_kagome[2] + return (x, y) + + def build_unit_cell_tensors_kagome(self): + sites = dict() + coords = [] + for coord_kagome in self.kagome_tensors.keys(): + tmp_coord = self.kagome_vertex_to_vertex(coord_kagome) + if torch.all(torch.tensor(coords != tmp_coord, dtype=torch.bool)): + coords.append(tmp_coord) + if self.restrictions: + for coord in coords: + self.kagome_tensors[(coord[0], coord[1], 1)] = self.kagome_tensors[(coord[0], coord[1], 0)] + self.kagome_tensors[(coord[0], coord[1], 2)] = self.kagome_tensors[(coord[0], coord[1], 0)] + self.kagome_tensors[(coord[0], coord[1], 4)] = self.kagome_tensors[(coord[0], coord[1], 3)] + for coord in coords: + tmp_tensor_transpose = (self.kagome_tensors[(coord[0], coord[1], 0)].transpose(1, 2)).contiguous() + t_a = (self.kagome_tensors[(coord[0], coord[1], 0)] + tmp_tensor_transpose) / 2 + t_b = t_a + t_c = t_a + t_r1 = self.kagome_tensors[(coord[0], coord[1], 3)] + t_r2 = t_r1 + sites[coord] = torch.einsum('akl,lno,bmn,cpo,qrp->abckrqm', t_a, t_r1, t_b, t_c, t_r2).flatten(start_dim=0, end_dim=2) + else: + for coord in coords: + t_a = self.kagome_tensors[(coord[0], coord[1], 0)] + t_b = self.kagome_tensors[(coord[0], coord[1], 1)] + t_c = self.kagome_tensors[(coord[0], coord[1], 2)] + t_r1 = self.kagome_tensors[(coord[0], coord[1], 3)] + t_r2 = self.kagome_tensors[(coord[0], coord[1], 4)] + sites[coord] = torch.einsum('akl,lno,bmn,cop,pqr->abckrqm', t_a, t_r1, t_b, t_c, t_r2).flatten(start_dim=0, end_dim=2) + + return sites + + def add_noise(self, noise): + r""" + :param noise: magnitude of the noise + :type noise: float + + Take IPEPS and add random uniform noise with magnitude ``noise`` to all on-site tensors + """ + if self.dtype == torch.float64: + for coord_kagome in self.kagome_tensors.keys(): + rand_t = torch.rand(self.kagome_tensors[coord_kagome].size(), dtype=self.dtype, device=self.device) + self.kagome_tensors[coord_kagome] = self.kagome_tensors[coord_kagome] + noise * rand_t + elif self.dtype == torch.complex128: + for coord_kagome in self.kagome_tensors.keys(): + rand_t_real = torch.rand(self.kagome_tensors[coord_kagome].size(), dtype=self.dtype, device=self.device) + rand_t_img = torch.rand(self.kagome_tensors[coord_kagome].size(), dtype=self.dtype, device=self.device) + self.kagome_tensors[coord_kagome] = self.kagome_tensors[coord_kagome] + noise * rand_t_real + self.kagome_tensors[coord_kagome] = self.kagome_tensors[coord_kagome] + noise * rand_t_img * 1j + else: + raise Exception("Unsuppoted data type. Optional: \"float64\", \"complex128\".") + + def get_physical_dim(self): + phys_dim = None + for coord_kagome, t in self.kagome_tensors.items(): + if phys_dim == None: + if coord_kagome[2] == 0 or coord_kagome[2] == 1 or coord_kagome[2] == 2: + phys_dim = t.size()[0] + return phys_dim + + def get_aux_bond_dims(self): + bond_dims = None + coords = [] + for coord_kagome in self.kagome_tensors.keys(): + tmp_coord = self.kagome_vertex_to_vertex(coord_kagome) + if torch.all(torch.tensor(coords != tmp_coord, dtype=torch.bool)): + coords.append(tmp_coord) + for coord in coords: + t_a = self.kagome_tensors[(coord[0], coord[1], 0)] + t_b = self.kagome_tensors[(coord[0], coord[1], 1)] + t_r2 = self.kagome_tensors[(coord[0], coord[1], 4)] + if bond_dims == None: + bond_dims = [t_b.size()[1], t_r2.size()[1], t_r2.size()[2], t_a.size()[1]] + + return bond_dims + + def write_to_file(self, outputfile, tol=1.0e-14, normalize=False): + write_ipeps_kagome(self, outputfile, tol=tol, normalize=normalize) + + +def extend_bond_dim_kagome(state, new_d): + coords = [] + new_state = state + for coord_kagome in state.kagome_tensors.keys(): + tmp_coord = state.kagome_vertex_to_vertex(coord_kagome) + if torch.all(torch.tensor(coords != tmp_coord, dtype=torch.bool)): + coords.append(tmp_coord) + for coord in coords: + for i in [0, 1, 2, 3, 4]: + if i < 3: + t = state.kagome_tensors[(coord[0], coord[1], i)] + dims = t.size() + new_t = torch.zeros((dims[0], new_d, new_d), dtype=state.dtype, device=state.device) + new_t[:, :dims[1], :dims[2]] = t + new_state.kagome_tensors[(coord[0], coord[1], i)] = new_t + else: + t = state.kagome_tensors[(coord[0], coord[1], i)] + dims = t.size() + new_t = torch.zeros((new_d, new_d, new_d), dtype=state.dtype, device=state.device) + new_t[:dims[0], :dims[1], :dims[2]] = t + new_state.kagome_tensors[(coord[0], coord[1], i)] = new_t + new_state.build_unit_cell_tensors_kagome() + + return new_state + + +def read_ipeps_kagome(jsonfile, restrictions=False, vertexToSite=None, peps_args=cfg.peps_args,\ + global_args=cfg.global_args): + r""" + :param jsonfile: input file describing iPEPS in json format + :param vertexToSite: function mapping arbitrary vertex of a square lattice + into a vertex within elementary unit cell + :param aux_seq: array specifying order of auxiliary indices of on-site tensors stored + in `jsonfile` + :param peps_args: ipeps configuration + :param global_args: global configuration + :type jsonfile: str or Path object + :type vertexToSite: function(tuple(int,int))->tuple(int,int) + :type peps_args: PEPSARGS + :type global_args: GLOBALARGS + :return: wavefunction + :rtype: IPEPS + + + A simple PBC ``vertexToSite`` function is used by default + """ + # asq = [x+1 for x in aux_seq] + kagome_tensors = OrderedDict() + with open(jsonfile) as j: + raw_state = json.load(j) + # read the list of considered kagome-type tensors + for ts in raw_state["map"]: + coord_kagome = (ts["x"],ts["y"],ts["z"]) + t = None + for s in raw_state["kagome_sites"]: + if s["siteId"] == ts["siteId"]: + t = s + if t == None: + raise Exception("Tensor with siteId: "+ts["sideId"]+" NOT FOUND in \"sites\"") + + if "format" in t.keys(): + if t["format"] == "1D": + X = torch.from_numpy(read_bare_json_tensor_np(t)) + else: + # default + X = torch.from_numpy(read_bare_json_tensor_np_legacy(t)) + + kagome_tensors[coord_kagome] = X.to(device=global_args.device) + + # Unless given, construct a function mapping from + # any site of square-lattice back to unit-cell + if vertexToSite == None: + # check for legacy keys + lX = 0 + lY = 0 + lX = raw_state["sizeM"] if "sizeM" in raw_state else raw_state["lX"] + lY = raw_state["sizeN"] if "sizeN" in raw_state else raw_state["lY"] + + def vertexToSite(coord): + x = coord[0] + y = coord[1] + return ( (x + abs(x)*lX)%lX, (y + abs(y)*lY)%lY ) + + state = IPEPS_KAGOME(kagome_tensors=kagome_tensors, restrictions=restrictions, vertexToSite=vertexToSite, \ + lX=lX, lY=lY, peps_args=peps_args, global_args=global_args) + else: + state = IPEPS_KAGOME(kagome_tensors=kagome_tensors, restrictions=restrictions, vertexToSite=vertexToSite, \ + peps_args=peps_args, global_args=global_args) + return state + +def write_ipeps_kagome(state, outputfile, tol=1.0e-14, normalize=False): + r""" + :param state: wavefunction to write out in json format + :param outputfile: target file + :param aux_seq: array specifying order in which the auxiliary indices of on-site tensors + will be stored in the `outputfile` + :param tol: minimum magnitude of tensor elements which are written out + :param normalize: if True, on-site tensors are normalized before writing + :type state: IPEPS + :type ouputfile: str or Path object + :type tol: float + :type normalize: bool + + TODO implement cutoff on elements with magnitude below tol + """ + json_state = dict({"lX": state.lX, "lY": state.lY, "kagome_sites": []}) + + def vertexToSite(coord): + x = coord[0] + y = coord[1] + return ((x + abs(x) * state.lX) % state.lX, (y + abs(y) * state.lY) % state.lY) + + site_ids = [] + site_map = [] + # The default input/output order + tensor_names = ['A', 'B', 'C', 'RD', 'RU'] + for coord_kagome, kagome_tensor in state.kagome_tensors.items(): + if normalize: + kagome_tensor = kagome_tensor / kagome_tensor.abs().max() + tensor_name = tensor_names[coord_kagome[2]%5] + nid = coord_kagome[0] + coord_kagome[1] * state.lX + site_ids.append(f"{tensor_name}{nid}") + site_map.append(dict({"siteId": site_ids[-1], "x": coord_kagome[0], "y": coord_kagome[1], "z": coord_kagome[2]})) + + json_tensor = serialize_bare_tensor_legacy(kagome_tensor) + json_tensor["siteId"]=site_ids[-1] + json_state["kagome_sites"].append(json_tensor) + + json_state["siteIds"] = site_ids + json_state["map"] = site_map + + with open(outputfile, 'w') as f: + json.dump(json_state, f, indent=4, separators=(',', ': ')) + diff --git a/models/kagome.py b/models/kagome.py new file mode 100644 index 00000000..731b9f8b --- /dev/null +++ b/models/kagome.py @@ -0,0 +1,242 @@ +import torch +import groups.su2 as su2 +import groups.su3 as su3 +import config as cfg +from ctm.generic.env import ENV +from ctm.generic import corrf +from ctm.generic import rdm_kagome # modified by Yi, 06/15/21 +from math import sqrt +from tn_interface import einsum, mm +from tn_interface import view, permute, contiguous +import itertools +import numpy as np + + +def _cast_to_real(t): + return t.real if t.is_complex() else t + + +class KAGOME(): + def __init__(self, phys_dim=3, j=0.0, k=1.0, h=0.0, global_args=cfg.global_args): + self.dtype = global_args.torch_dtype + self.device = global_args.device + self.phys_dim = phys_dim + self.j = j # 2-site permutation coupling constant + self.k = k # 3-site ring exchange coupling constant + self.h = h # 3-site ring exchange coupling constantNN coupling between A & C in the same unit cell + + # For now, only one-site + self.obs_ops = self.get_obs_ops() + self.h2, self.h3_l, self.h3_r = self.get_h() + + def get_obs_ops(self): + obs_ops = dict() + irrep = su2.SU2(self.phys_dim, dtype=self.dtype, device=self.device) + obs_ops["sz"] = irrep.SZ() + obs_ops["sp"] = irrep.SP() + obs_ops["sm"] = irrep.SM() + obs_ops["SS"] = irrep.SS() + irrep_su3_def = su3.SU3_DEFINING(dtype=self.dtype, device=self.device) + obs_ops["tz"] = irrep_su3_def.TZ() + obs_ops["tp"] = irrep_su3_def.TP() + obs_ops["tm"] = irrep_su3_def.TM() + obs_ops["vp"] = irrep_su3_def.VP() + obs_ops["vm"] = irrep_su3_def.VM() + obs_ops["up"] = irrep_su3_def.UP() + obs_ops["um"] = irrep_su3_def.UM() + obs_ops["y"] = irrep_su3_def.Y() + # obs_ops["C1"] = irrep_su3_def.C1() + # obs_ops["C2"] = irrep_su3_def.C2() + + return obs_ops + + def get_h(self): + pd = self.phys_dim + irrep = su2.SU2(pd, dtype=self.dtype, device=self.device) + # identity operator on two spin-S spins + idp2x2 = torch.eye(pd ** 2, dtype=self.dtype, device=self.device) + SS = irrep.SS() + SS = SS.view(pd ** 2, pd ** 2) + perm_2site = SS + SS @ SS - idp2x2 + # Reshape back into rank-4 tensor for later use with reduced density matrices + perm_2site = perm_2site.view(pd, pd, pd, pd).contiguous() + + h2 = perm_2site + h3_l = torch.einsum('ijal,lkbc->ijkabc', perm_2site, perm_2site) + h3_r = torch.einsum('ijal,klbc->ikjabc', perm_2site, perm_2site) + + return h2, h3_l, h3_r + + def energy_1site(self, state, env): + pd = self.phys_dim + energy = 0.0 + idp = torch.eye(pd, dtype=self.dtype, device=self.device) + for coord, site in state.sites.items(): + # intra-cell 2-site exchange terms + norm = rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ia,jb,kc->ijkabc', idp, idp, idp)) + energy += self.j * rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ijab,kc->ijkabc', self.h2, idp)) / norm + energy += self.j * rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('jkbc,ia->ijkabc', self.h2, idp)) / norm + energy += self.j * rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ikac,jb->ijkabc', self.h2, idp)) / norm + + # intra-cell 3-site ring exchange terms + energy += (self.k + self.h * 1j) * rdm_kagome.trace1x1_kagome(coord, state, env, self.h3_l) / norm + energy += (self.k - self.h * 1j) * rdm_kagome.trace1x1_kagome(coord, state, env, self.h3_r) / norm + + # inter-cell 2-site exchange terms + rdm2x2_ab = rdm_kagome.rdm2x2_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=(), + sites_to_keep_01=(), sites_to_keep_11=('A')) + energy += self.j * torch.einsum('ilad,ilad', rdm2x2_ab, self.h2) + rdm1x2_bc = rdm_kagome.rdm2x1_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=('C')) + energy += self.j * torch.einsum('ijab,ijab', rdm1x2_bc, self.h2) + rdm2x1_ac = rdm_kagome.rdm1x2_kagome(coord, state, env, sites_to_keep_00=('C'), sites_to_keep_01=('A')) + energy += self.j * torch.einsum('ijab,ijab', rdm2x1_ac, self.h2) + + # inter-cell 3-site ring exchange terms + rdm2x2_ring = rdm_kagome.rdm2x2_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=('C'), + sites_to_keep_01=(), sites_to_keep_11=('A')) + energy += (self.k + self.h * 1j) * torch.einsum('ijlabd,lijdab', rdm2x2_ring, self.h3_l) + energy += (self.k - self.h * 1j) * torch.einsum('ijlabd,lijdab', rdm2x2_ring, self.h3_r) + + energy_per_site = energy / (len(state.sites.items()) * 3.0) + energy_per_site = _cast_to_real(energy_per_site) + return energy_per_site + + def eval_obs(self, state, env): + pd = self.phys_dim + chirality = self.h3_l - self.h3_r + idp = torch.eye(pd, dtype=self.dtype, device=self.device) + obs = dict() + with torch.no_grad(): + for coord, site in state.sites.items(): + norm = rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ia,jb,kc->ijkabc', idp, idp, idp)) + obs["chirality_dn"] = rdm_kagome.trace1x1_kagome(coord, state, env, chirality) / norm + obs["chirality_dn"] = _cast_to_real(obs["chirality_dn"] * 1j) + rdm2x2_ring = rdm_kagome.rdm2x2_kagome(coord, state, env, sites_to_keep_00=('B'), + sites_to_keep_10=('C'), sites_to_keep_01=(), + sites_to_keep_11=('A')) + obs["chirality_up"] = torch.einsum('ijlabd,lijdab', rdm2x2_ring, chirality) + obs["chirality_up"] = _cast_to_real(obs["chirality_up"] * 1j) + + obs["avg_bonds_dn"] = rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ijab,kc->ijkabc', self.h2, idp)) / norm + obs["avg_bonds_dn"] += rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('jkbc,ia->ijkabc', self.h2, idp)) / norm + obs["avg_bonds_dn"] += rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ikac,jb->ijkabc', self.h2, idp)) / norm + obs["avg_bonds_dn"] = _cast_to_real(obs["avg_bonds_dn"]) / 3.0 + + rdm2x2_ab = rdm_kagome.rdm2x2_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=(), + sites_to_keep_01=(), sites_to_keep_11=('A')) + obs["avg_bonds_up"] = torch.einsum('ilad,ilad', rdm2x2_ab, self.h2) + rdm1x2_bc = rdm_kagome.rdm2x1_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=('C')) + obs["avg_bonds_up"] += torch.einsum('ijab,ijab', rdm1x2_bc, self.h2) + rdm2x1_ac = rdm_kagome.rdm1x2_kagome(coord, state, env, sites_to_keep_00=('C'), sites_to_keep_01=('A')) + obs["avg_bonds_up"] += torch.einsum('ijab,ijab', rdm2x1_ac, self.h2) + obs["avg_bonds_up"] = _cast_to_real(obs["avg_bonds_up"]) / 3.0 + + # prepare list with labels and values + obs_labels = ["avg_bonds_dn", "avg_bonds_up", "chirality_dn", "chirality_up"] + obs_values = [obs[label] for label in obs_labels] + return obs_values, obs_labels + + def energy_1triangle(self, state, env): + energy_dn = 0.0 + energy_up = 0.0 + pd = self.phys_dim + idp = torch.eye(pd, dtype=self.dtype, device=self.device) + for coord, site in state.sites.items(): + # intra-cell 2-site exchange terms + norm = rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ia,jb,kc->ijkabc', idp, idp, idp)) + energy_dn += self.j * rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ijab,kc->ijkabc', self.h2, idp)) / norm + energy_dn += self.j * rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('jkbc,ia->ijkabc', self.h2, idp)) / norm + energy_dn += self.j * rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ikac,jb->ijkabc', self.h2, idp)) / norm + # intra-cell 3-site ring exchange terms + energy_dn += (self.k + self.h * 1j) * rdm_kagome.trace1x1_kagome(coord, state, env, self.h3_l) / norm + energy_dn += (self.k - self.h * 1j) * rdm_kagome.trace1x1_kagome(coord, state, env, self.h3_r) / norm + + # inter-cell 2-site exchange terms + rdm2x2_ab = rdm_kagome.rdm2x2_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=(), + sites_to_keep_01=(), sites_to_keep_11=('A')) + energy_up += self.j * torch.einsum('ilad,ilad', rdm2x2_ab, self.h2) + rdm1x2_bc = rdm_kagome.rdm2x1_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=('C')) + energy_up += self.j * torch.einsum('ijab,ijab', rdm1x2_bc, self.h2) + rdm2x1_ac = rdm_kagome.rdm1x2_kagome(coord, state, env, sites_to_keep_00=('C'), sites_to_keep_01=('A')) + energy_up += self.j * torch.einsum('ijab,ijab', rdm2x1_ac, self.h2) + + # inter-cell 3-site ring exchange terms + rdm2x2_up = rdm_kagome.rdm2x2_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=('C'), + sites_to_keep_01=(), sites_to_keep_11=('A')) + energy_up += (self.k + self.h * 1j) * torch.einsum('ijlabd,lijdab', rdm2x2_up, self.h3_l) + energy_up += (self.k - self.h * 1j) * torch.einsum('ijlabd,lijdab', rdm2x2_up, self.h3_r) + energy_dn = energy_dn / (len(state.sites.items()) * 3.0) + energy_up = energy_up / (len(state.sites.items()) * 3.0) + return energy_dn, energy_up + + def eval_generators(self, state, env): + pd = self.phys_dim + gens = dict({"generators_A": torch.zeros(8, dtype=self.dtype), "generators_B": torch.zeros(8, dtype=self.dtype), + "generators_C": torch.zeros(8, dtype=self.dtype)}) + with torch.no_grad(): + for coord, site in state.sites.items(): + for site_type in ['A', 'B', 'C']: + rdm1x1 = rdm_kagome.rdm1x1_kagome(coord, state, env, sites_to_keep=(site_type)).view(pd, pd) + tmp = dict() + for label in ["tz", "tp", "tm", "vp", "vm", "up", "um", "y"]: + op = self.obs_ops[label] + tmp[f"{label}{coord}"] = einsum('ij,ji', rdm1x1, op) + gens[f"generators_{site_type}"] += np.array( + [tmp[f"tz{coord}"], tmp[f"tp{coord}"], tmp[f"tm{coord}"], + tmp[f"vp{coord}"], tmp[f"vm{coord}"], + tmp[f"up{coord}"], tmp[f"um{coord}"], tmp[f"y{coord}"]]) + for site_type in ['A', 'B', 'C']: + gens[f"generators_{site_type}"] /= len(state.sites.items()) + return gens + + def eval_C2(self, state, env): + pd = self.phys_dim + idp = torch.eye(pd, dtype=self.dtype, device=self.device) + irrep_su3_def = su3.SU3_DEFINING(dtype=self.dtype, device=self.device) + c2 = irrep_su3_def.C2() + c2_list = dict({"C2_dn": 0., "C2_up": 0.}) + with torch.no_grad(): + for coord, site in state.sites.items(): + norm = rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ia,jb,kc->ijkabc', idp, idp, idp)) + c2_list["C2_dn"] += rdm_kagome.trace1x1_kagome(coord, state, env, c2) / norm + + rdm2x2_up = rdm_kagome.rdm2x2_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=('C'), + sites_to_keep_01=(), sites_to_keep_11=('A')) + c2_list["C2_up"] += torch.einsum('ijlabd,lijdab', rdm2x2_up, c2) + + c2_list["C2_dn"] /= len(state.sites.items()) + c2_list["C2_up"] /= len(state.sites.items()) + return c2_list + + def eval_C1(self, state, env): + pd = self.phys_dim + idp = torch.eye(pd, dtype=self.dtype, device=self.device) + irrep_su3_def = su3.SU3_DEFINING(dtype=self.dtype, device=self.device) + c1 = irrep_su3_def.C1() + c1_dict = dict({"C1_AB_dn": 0., "C1_BC_dn": 0., "C1_AC_dn": 0., "C1_AB_up": 0., "C1_BC_up": 0., "C1_AC_up": 0.}) + with torch.no_grad(): + for coord, site in state.sites.items(): + norm = rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ia,jb,kc->ijkabc', idp, idp, idp)) + c1_dict["C1_AB_dn"] += rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ijab,kc->ijkabc', c1, idp)) / norm + c1_dict["C1_BC_dn"] += rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('jkbc,ia->ijkabc', c1, idp)) / norm + c1_dict["C1_AC_dn"] += rdm_kagome.trace1x1_kagome(coord, state, env, torch.einsum('ikac,jb->ijkabc', c1, idp)) / norm + + rdm2x2_ab = rdm_kagome.rdm2x2_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=(), + sites_to_keep_01=(), sites_to_keep_11=('A')) + c1_dict["C1_AB_up"] += torch.einsum('ilad,ilad', rdm2x2_ab, c1) + rdm1x2_bc = rdm_kagome.rdm2x1_kagome(coord, state, env, sites_to_keep_00=('B'), sites_to_keep_10=('C')) + c1_dict["C1_BC_up"] += torch.einsum('ijab,ijab', rdm1x2_bc, c1) + rdm2x1_ac = rdm_kagome.rdm1x2_kagome(coord, state, env, sites_to_keep_00=('C'), sites_to_keep_01=('A')) + c1_dict["C1_AC_up"] += torch.einsum('ijab,ijab', rdm2x1_ac, c1) + + c1_dict["C1_AB_dn"] /= len(state.sites.items()) + c1_dict["C1_BC_dn"] /= len(state.sites.items()) + c1_dict["C1_AC_dn"] /= len(state.sites.items()) + c1_dict["C1_AB_up"] /= len(state.sites.items()) + c1_dict["C1_BC_up"] /= len(state.sites.items()) + c1_dict["C1_AC_up"] /= len(state.sites.items()) + c1_dict["total_C1_dn"] = c1_dict["C1_AB_dn"] + c1_dict["C1_BC_dn"] + c1_dict["C1_AC_dn"] + c1_dict["total_C1_up"] = c1_dict["C1_AB_up"] + c1_dict["C1_BC_up"] + c1_dict["C1_AC_up"] + + return c1_dict +