How to calculate string operator in the S=1 Heisenberg antiferromagnetic chain in ITensor C++ version ?

Dear ITensor teams,
Thanks for the reply. I have a question about ITensor C++ version: How to calculate string operator in the S=1 Heisenberg antiferromagnetic chain ? e.g. <\psi_{0}|S^{z}{1}exp(i \pi \sum{l=2}^{L-1}S^{z}{l})S^{z}{L}|\psi_{0}>, here |\psi_{0}> is the ground state of spin-1 antiferromagnetic chain.
We know the benchmark result is : \lim_{|L->\infty|}<\psi_{0}|S^{z}{1}exp(i \pi \sum{l=2}^{L-1}S^{z}{l})S^{z}{L}|\psi_{0}>=-4/9
The following i write a naive code but the results is wrong (for L=24, the result = -0.28 )

auto sites = SpinOne(L,{"ConserveQNs",false});
auto ampo = AutoMPO(sites);
    // OBC
    for(int i = 1; i< L; ++i)
            {
                ampo += 0.5*JJ,"S+",i,"S-",i+1;
                ampo += 0.5*JJ,"S-",i,"S+",i+1;
                ampo += JJ,"Sz",i,"Sz",i+1;
            }
    auto H = toMPO(ampo);

    auto psi0 = randomMPS(sites);
    auto sweeps = Sweeps(100); //number of sweeps is 5
    sweeps.maxdim() = 128;
    sweeps.cutoff() = 1E-8;

    // Begin the DMRG calculation
    auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet",quiet,"EnergyErrgoal",1E-5});

    //Measure a string correlation function
    //at sites i and j
    //Below we will assume j > i
    auto i = 1;
    auto j = L;
    //Make the operators you want to measure
    auto op_i = op(sites,"Sz",i);
    auto op_j = op(sites,"Sz",j);

    //'gauge' the MPS to site i
    //any 'position' between i and j, inclusive, would work here
    psi.position(i); 

    //Create the bra/dual version of the MPS psi
    auto psidag = dag(psi);

    //Prime the link indices to make them distinct from
    //the original ket links
    psidag.prime("Link");

    //index linking i-1 to i:
    auto li_1 = leftLinkIndex(psi,i);

    auto C = prime(psi(i),li_1)*op_i;
    C *= prime(psidag(i),"Site");
for(int k = i+1; k < j; ++k)
    {
    C *= psi(k);
    C *= expHermitian(PI*op(sites,"Sz",k),1_i); // construct exponential operator 
    C *= prime(psidag(k),"Site");
    }
//index linking j to j+1:
auto lj = rightLinkIndex(psi,j);

C *= prime(psi(j),lj)*op_j;
C *= prime(psidag(j),"Site");

auto result = eltC(C); 

printfln("The String operator is = ",real(result));

Maybe there are some wrong in the code, please help me figure it out, many thanks!
Best regards,
Sugar