tntMpsMpo.c

/*
Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
Date:    $LastChangedDate$
(c) University of Oxford 2014
*/

/* Include the internal header for the TNT library */
#include "tntMpsInternal.h"


void tntMpsMpoConnect(tntNetwork mps, 
                      tntNetwork mpo) 
{
    
    tntNode A, O; /*  the current MPS node and MPO node */
    unsigned L; /* Length of system */
    
    /* First check that the MPS and the MPO have the same number of terms */
    L = tntMpsLength(mps);
    
    if (tntMpsLength(mpo) != L) {
        tntErrorPrint("Cannot connect the MPS to the MPO as there are not the same number of nodes in the two networks|" /* NO_COVERAGE */
                         "MPS contains %d nodes, MPO contains %d nodes",L,tntMpsLength(mpo)); /* NO_COVERAGE */
    } /* NO_COVERAGE */
    
    /* Make a copy of the MPO network */
    mpo = tntNetworkCopy(mpo);
    
    /* Find the MPO node for site 0 */
    O = tntNodeFindFirst(mpo);
    
    /* Convert the network to a node group, leaving terminating connections intact */
    tntNetworkToNodeGroup(&mpo,0);
    
    /* Find the MPS node corresponding to site 0 */
    A  = tntNodeFindFirst(mps);
        
    /* Now loop through all sites, attaching the copy of the MPO term to the MPS term */
    while (1) {
        
        /* Join a copy of the MPO to the MPS */
        tntNodeJoin(O,"U",A,"D");
        
        /* Find the next group of sites if they exist, or exit the loop if the last site has been reached */
        if (A == tntNodeFindLast(mps)) {
            break;
        } else {
            /* Find the next nodes along */
            A = tntNodeFindConn(A, "R");
            O = tntNodeFindConn(O, "R");
        }
    }
}

/* To add once refs working  * The routine performs the zipping algorithm described in \cite Stoudenmire2010 . */
double tntMpsMpoContract(tntNetwork mpsmpo, 
                         int chi) 
{
    
    tntNode A, O;       /* The current node of the MPS and the MPO */
    tntNode C, Cc;      /* The node that is the result of contracting the MPS and its MPO node, and a copy of it taken to SVD */
    tntNode U, Udag;    /* First unitary that is a result of SVDing tnCc, and its conjugate */
    unsigned L, j;      /* The length of the MPS network, site counter */
    double err_tot=0.0, err; /* total truncation error, and truncation error per SVD. */
    
    /* Find the length of the network */
    L = tntMpsLength(mpsmpo);
    
    /* ------ contracting from the left to the right, truncating to truncation tolerace, to result in a single MPS ---- */
    
    /* Find the nodes corresponding to site 0 */
    A = tntNodeFindFirst(mpsmpo);
    O = tntNodeFindConn(A, "D");
    
    /* Contract these two nodes */
    C = tntNodeContract(A, O, NULL, "LR=MS");
    
    /* Fuse the extra leg coming from the MPO with the MPS leg */
    tntNodeFuse(C,"LM","L","");
    
    for (j = 0; j < L-1; j++) {
        
        /* Take a copy of tnC - this copy will not be connected to any other nodes in the network */
        Cc = tntNodeCopy(C);
        
        /* Take an untruncated SVD of Cc, assigning the result back to the copy. */
        Cc = tntNodeSVD(Cc, "LD", "R", "L", "R", "L", -1, &err);
        err_tot += err;
        
        /* Make an unconnected copy of the first unitary - this will become the new MPS node to replace tnA when inserted back in the network */
        U = tntNodeCopy(Cc);
        
        /* the rest of the nodes are not needed -> discard them */
        tntNodeGroupFree(&Cc);
        
        /* Legs L and D of U define one dimension (lets call rows) of the Unitary matrix, leg R of U defines the other dimension (lets call columns) of the unitary matrix */
        
        /* Take a conjugate copy of tnU */
        Udag = tntNodeCopy(U,1);
        
        /* Now map the legs of tnUdag so that leg2 -> leg 1 and leg 1 defines rows, and leg 1 -> leg 2 and leg 2 and 3 define columns. 
         * tnUdag is now hermitian conjuate of tnU, and when leg 2 of tnU contracts with leg 1 of tnUdag the result will be the identity */
        tntNodeMapLegs(Udag, "LR=RL");
        
        /* Now insert --tnU--tnUdag-- = the identity (when the physical legs are fused with the outer legs), back into the MPS network */
        /* Take a different action depending on whether we are at the start of the network, or in the middle of the network */
        if (0 == j) {
            /* if we inserting are at the beginning of the network, we can put insert the nodes directly, without temporarily (and illegally) connecting legs with non-matching dimensions */
            
            /* Insert the node at the start, which connects leg 2 of tnUdag to leg 1 of tnW */
            tntNodeInsertAtStart(Udag, "L", "R", mpsmpo);
            
            /* Note that we defined legs 1 and 3 of tnC as rows, and legs 2 and 3 of tnUdag as columns - for a full insertion we also need to contract their physical legs */
            tntNodeJoin(Udag, "D", C, "D");
            
            /* Now insert U at the start of the network. U becomes the new MPS node for this set (the previous MPS node is part of C) */
            tntNodeInsertAtStart(U, "L", "R", mpsmpo);
        } else {
            /* if we are inserting in the middle of the network, we can't simply insert the nodes in, as this would cause a mismatch between dimensions of connecting legs */
            /* This operation is thus forbidden and would cause an error in the program */
            
            /* Instead, split the connections between tnC and the previous MPS node, and then join up nodes one by one */
            /* first identify the previous node */
            A = tntNodeFindConn(C, "L");
            
            /* Split connections between tnA and tnC */
            tntNodeSplit(A, C);
            
            /* connect tnA to tnU */
            tntNodeJoin(A, "R", U, "L");
            
            /* connect tnU to tnUdag */
            tntNodeJoin(U, "R", Udag, "L");
            
            /* Connect tnUdag to tnC, remembering they are joined on 2 legs */
            tntNodeJoin(Udag, "R", C, "L");
            tntNodeJoin(Udag, "D", C, "D");
            
            /* the network is now complete again */
        }
        
        /* find the next MPS node along, and it's corresponding MPO node */
        A = tntNodeFindConn(C, "R");
        O = tntNodeFindConn(A, "D");
        
        /* Now contract all these nodes */
        C = tntNodeListContract("LRDS",Udag,C,A,O);
        
        /* If tnA was the last node, tnC simply forms the final MPS node.
         * Otherwise continue the loop to form further MPS nodes */
    }
    
    /* Fuse the extra leg coming from the MPO with the MPS leg */
    tntNodeFuse(C,"RS","R","");
    
    /* Now sweep back right to left, this time truncating to the maximum internal dimension chi */
    /* Note setting desired orthogonality centre to zero enures a full right to left sweep */
    err = tntMpsTruncate(mpsmpo,0,chi);
    err_tot += err;
    
    /* return the total truncation error */
    return err_tot;

}


double tntMpsMpoProduct(tntNetwork mps, 
                        tntNetwork mpo, 
                        int chi) 
{

    /* Create the full network */
    tntMpsMpoConnect(mps,mpo);
    
    /* Contract the full network and return the error */
    return tntMpsMpoContract(mps,chi);

}