Hi Miles,

Sure, here are the code snippets

```
function Ham_constructor(sites, coeff, bnd_str, on_site =false, D_vec=nothing)
i=1
target = i+1
# Input operator terms which define a Hamiltonian
ampo = OpSum()
while i <= (length(sites))
if (bnd_str == "pbc" && i == length(sites))
target = 1
elseif (bnd_str == "obc" && i == length(sites))
break
end
ampo += (0.5*coeff[1], "Sz * Sz", i)
ampo += (0.5*coeff[1], "Sy * Sy", i)
ampo += (0.5*coeff[1], "Sx * Sx", i)
ampo += (coeff[2], "Sz", i, "Sz", target)
ampo += (coeff[3], "Sx", i, "Sx", target)
ampo += (coeff[3], "Sy", i, "Sy", target)
ampo += (coeff[4], "Sz * Sz", i, "Sz * Sz", target)
ampo += (coeff[5], "Sx * Sx", i, "Sx * Sx", target)
ampo += (coeff[5], "Sx * Sy", i, "Sx * Sy", target)
ampo += (coeff[5], "Sy * Sx", i, "Sy * Sx", target)
ampo += (coeff[5], "Sy * Sy", i, "Sy * Sy", target)
ampo += (coeff[6], "Sz * Sx", i, "Sz * Sx", target)
ampo += (coeff[6], "Sz * Sy", i, "Sz * Sy", target)
ampo += (coeff[6], "Sx * Sz", i, "Sx * Sz", target)
ampo += (coeff[6], "Sy * Sz", i, "Sy * Sz", target)
ampo += (coeff[7], "Sx * Sx", i, "Sy * Sy", target)
ampo += (-coeff[7], "Sx * Sy", i, "Sy * Sx", target)
ampo += (-coeff[7], "Sy * Sx", i, "Sx * Sy", target)
ampo += (coeff[7], "Sy * Sy", i, "Sx * Sx", target)
if (on_site == true)
println(D_vec[i])
ampo += (D_vec[i], "Sz * Sz", i)
ampo += (-D_vec[i]*(1/3), "Sz * Sz", i)
ampo += (-D_vec[i]*(1/3), "Sy * Sy", i)
ampo += (-D_vec[i]*(1/3), "Sx * Sx", i)
end
i+=1
target = i+1
end
# Convert these terms to an MPO tensor network
H = MPO(ampo, sites)
return H
end
function DMRG_driver(sites, H, state_reqd, init_state=nothing, gr_state=nothing, weights=nothing)
let
# Create an initial random matrix product state
#Random.seed!(1234)
#psi0 = randomMPS(sites, 60)
state = [isodd(i) ? "Up" : "Dn" for i=1:length(sites)]
#state = ["Z0" for i=1:length(sites)]
psi0 = productMPS(sites,state)
@show flux(psi0)
# Plan to do 40 DMRG sweeps:
sweeps = Sweeps(60)
# Set maximum MPS bond dimensions for each sweep
#setmaxdim!(sweeps, 10, 30, 50, 80, 100, 150, 200,300)
setmaxdim!(sweeps, 10, 10, 10, 20, 20, 30, 30, 50, 80, 80, 100, 100, 150, 180, 200, 250, 300)
# Set maximum truncation error allowed when adapting bond dimensions
setcutoff!(sweeps, 1e-5)
setnoise!(sweeps, 1e-5, 1e-5, 1e-7, 1e-7, 1e-9, 1e-9, 1e-10, 1e-11, 1e-12, 1e-12)
@show sweeps
# Run the DMRG algorithm, returning energy and optimized MPS
if state_reqd == "gr"
energy, psi = dmrg(H, psi0, sweeps)
elseif state_reqd == "1st_exc"
energy, psi = dmrg(H, [gr_state], init_state, sweeps; weights)
end
return energy, psi
end
end
function main(spin_type, no_spins, state_reqd, bnd_str)
# Create N spin-one degrees of freedom
#sites = siteinds("S=$spin_type", no_spins, conserve_qns=true)
sites = siteinds("S=$spin_type", no_spins)
# Alternatively can make spin-half sites instead
#sites = siteinds("S=1/2",no_spins)
D_tuple=IterTools.product([-8,0,8], [-8,0,8])
energy_arr = zeros(length(D_tuple))
var_arr = zeros(length(D_tuple))
ovlp_gr = zeros(length(D_tuple))
for (index, D_tup) in enumerate(D_tuple)
@show D_tup
#print("one begins", index)
D1 = D_tup[1]
D2 = D_tup[2]
D_vec = [D1 + ((-1)^i)*D2 for i in 0:length(sites)-1] #local anisotropy as parameter for Hamiltonian
#println("D array", D_vec)
H = Ham_constructor(sites,[-5.58, 9.53, -8.97, 1.27, 6.59, -3.18, 5.04], bnd_str, true, D_vec)
if state_reqd=="gr"
energy, psi = DMRG_driver(sites, H, state_reqd)
elseif state_reqd == "1st_exc"
energy_gr, psi_gr = DMRG_driver(sites, H, "gr")
psi_ex_init = psi_gr
weights = 300
energy, psi = DMRG_driver(sites, H, state_reqd, psi_ex_init, psi_gr, weights)
@show inner(psi, psi_gr)
end
# Compute energy variance
#************************
H_sq = inner(H,psi,H,psi)
H_avg = inner(psi', H , psi)
var = H_sq-H_avg^2
@printf("Var of final state = %.12f\n", real(var))
@printf("Avg energy in final state = %.12f\n", real(H_avg))
#@printf("Per exciton energy in final state computed from <H> = %.12f\n", real(H_avg)/length(sites))
@printf("Per exciton energy in final state computed directly = %.12f\n", real(energy)/length(sites))
energy_arr[index] = real(H_avg)/length(sites)
var_arr[index] = real(var)
ovlp_gr[index] = real(inner(psi, psi_gr))
**** REST OF THE CODE FOLLOWS ****
end
end
```

The outputs are as follows:

For no_spins = 10, S=1 spin chains

#
Parameter set (D1 = 0, D2 =-8)

#**************************************

#
EXCITED STATE NOT CONVERGED; RETURNS GROUND STATE

.

.

.

After sweep 59 energy=-36.72260723778795 maxlinkdim=4 maxerr=3.20E-07 time=0.029

After sweep 60 energy=-36.722607237787955 maxlinkdim=4 maxerr=3.20E-07 time=0.029

inner(psi, psi_gr) = 0.9999999999906718 + 6.730727086790012e-16im #NOT CONVERGED EXCITED STATE AS OVERLAP WITH GR STATE ABSOLUTE

Var of final state = 0.000825051831

Avg energy in final state = -37.722607237769

Per exciton energy in final state computed directly = -3.672260723779

#
Parameter set (D1 = -8, D2 =-8)

#**************************************

#
EXCITED STATE CONVERGED

.

.

.

After sweep 59 energy=-56.3847361997004 maxlinkdim=128 maxerr=9.86E-12 time=2.008

After sweep 60 energy=-56.38473619970033 maxlinkdim=128 maxerr=9.86E-12 time=2.021

inner(psi, psi_gr) = 4.98411099630891e-6 - 1.21952253850921e-10im # WELL CONVERGED EXCITED STATE AS INNER WITH GROUND STATE IS NEARLY ZERO

Var of final state = 0.000000374540

Avg energy in final state = -56.384736199725

Per exciton energy in final state computed directly = -5.638473619970