Untitled

avatar
unknown
plain_text
a year ago
7.4 kB
7
Indexable
using ITensors
using ITensorCorrelators
using Printf
using METTS
using LinearAlgebra
using ITensors.HDF5
using Random
using Serialization
using HDF5
using Dumper

L = parse(Int64, ARGS[1])
W = parse(Int64, ARGS[2])
U = parse(Float64, ARGS[3])
filling = parse(Float64, ARGS[4])
t = parse(Float64, ARGS[5])
T = parse(Float64, ARGS[6])
tau = parse(Float64, ARGS[7])
maxD = parse(Int64, ARGS[8])
cutoff = parse(Float64, ARGS[9])
seed = parse(Int64, ARGS[10])
nmetts = parse(Int64, ARGS[11])

function ITensors.op(::OpName"hole", ::SiteType"Electron", s::Index)
    nup = op("Nup", s)
    ndn = op("Ndn", s)
    nupdn = op("Nupdn", s)
    id = op("Id", s)
    hole_op = id - ndn - nup
    return hole_op
end

# Function to calculate average and error
function avg_err(v::Vector)
    N = length(v)
    avg = sum(v) / N
    avg2 = sum(v .^ 2) / N
    return avg, √((avg2 - avg^2) / N)
end

function load_samples(filename)
    if isfile(filename)
        samples = Dict()
        open(filename, "r") do file
            for line in eachline(file)
                step, sample_str = split(line, ":")
                sample_vec_str = strip(sample_str) 
                sample_vec = eval(Meta.parse(sample_vec_str))  
                samples[parse(Int, step)] = sample_vec
            end
        end
        return samples, maximum(keys(samples))
    else
        return Dict(), 0
    end
end

# Function to save samples
function save_samples(filename, samples)
    open(filename, "w") do file
        for (step, sample) in samples
            write(file, "$step: $sample\n")
        end
    end
end

# Define Fermi-Hubbard Hamiltonian
function fermi_hubbard_cylinder_mpo(L, W, t, U, sites)
    N = L * W
    lattice = square_lattice(L, W; yperiodic=true)
    os = OpSum()
    for b in lattice
        os += -t, "Cdagup", b.s1, "Cup", b.s2
        os += -t, "Cdagdn", b.s1, "Cdn", b.s2
        os += -t, "Cdagup", b.s2, "Cup", b.s1
        os += -t, "Cdagdn", b.s2, "Cdn", b.s1
    end
    for n = 1:N
        os += U, "Nupdn", n
    end
    H = MPO(os, sites)
    return H
end

# Main simulation function
function main(;L=4, W=4, U=10, filling=0.875, t=1.00, T=0.5, tau=0.1, maxD=250, cutoff=1e-6, seed=1, NMETTS=1000)
    ITensors.Strided.set_num_threads(1)
    BLAS.set_num_threads(1)
    ITensors.disable_threaded_blocksparse()
    filename = @sprintf("/data/condmat/asinha/Research/Projects/metts.cylinder.hubbard/ps.metts.hubbard/outfiles.metts/L.%d.W.%d/U.%.3f/filling.%.5f/T.%.5f/D.2000/tau.0.20.cutoff.%.1e.seed.%d/outfile.h5", L, W, U, filling, T, cutoff, seed)
    mkpath(dirname(filename))

    outfile = DumpFile(filename)
    Random.seed!(seed)
    beta = 1.0 / T
    N = L * W

    sites = siteinds("Electron", N; conserve_qns=true)
    H = fermi_hubbard_cylinder_mpo(L, W, t, U, sites)

    samples_filename = joinpath(dirname(filename), "samples.txt")
    samples, start_step = load_samples(samples_filename)
    @show start_step

    # Initialize state
    state = fill("Emp", N)
    if start_step > 0
        last_sample = samples[start_step]
        for (i, occupancy) in enumerate(last_sample)
            state[i] = occupancy == 1 ? "Emp" : occupancy == 2 ? "Up" : occupancy == 3 ? "Dn" : "UpDn"
        end
    else
        Ne = round(Int, filling * N)
        Ne % 2 == 0 || error("Ne must be even for equal number of up and down electrons")
        up_indices = Random.randperm(N)[1:Ne÷2]
        down_indices = setdiff(Random.randperm(N), up_indices)[1:Ne÷2]
        for idx in up_indices
            state[idx] = "Up"
        end
        for idx in down_indices
            state[idx] = "Dn"
        end
    end
    
    psi = randomMPS(sites, state)
    energies = Float64[]

    for step in (start_step + 1):NMETTS
        println("Making METTS number $step")
        psi = timeevo_tdvp_extend(H, psi, -(beta/ 2); tau=tau, cutoff=cutoff, normalize=true, solver_backend="applyexp", maxm=maxD, tol=1e-6, kkrylov=2, tau0=0.03, nsubdiv=3)

        GC.gc()

        println("measure energy")
        @time energy = ITensors.inner(psi', H, psi)
        println("measure Sz")
        @time szs = expect(psi, "Sz")
        println("measure Ntot")
        @time  Ntotal = expect(psi, "Ntot")
        println("measure Nupdn")
        @time  dob = expect(psi, "Nupdn")
        println("measure hole")
        @time  holes = expect(psi, "hole")
        println("measure SzSz")
        @time szszs = correlation_matrix(psi, "Sz", "Sz"; ishermitian=true)
        println("measure SpSm")
        @time spsms = correlation_matrix(psi, "S+", "S-"; ishermitian=false)
        println("measure SmSp")
        @time spsms = correlation_matrix(psi, "S-", "S+"; ishermitian=false)t
        println("measure NupNup")
        @time nupnup = correlation_matrix(psi, "Nup", "Nup"; ishermitian=true)
        println("measure NdnNdn")
        @time ndnndn = correlation_matrix(psi, "Ndn", "Ndn"; ishermitian=true)
        println("measure NupNdn")
        @time nupndn = correlation_matrix(psi, "Nup", "Ndn"; ishermitian=true)
        println("measure NdnNup")
        @time ndnnup = correlation_matrix(psi, "Ndn", "Nup"; ishermitian=true)
        println("measure NtotNtot")
        @time ntotntot = correlation_matrix(psi, "Ntot", "Ntot"; ishermitian=true)
        println("measure HoleHole")
        @time hole_corr = correlation_matrix(psi, "hole", "hole"; ishermitian=true)
        println("measure DoublonDoublon")
        @time doublon_doublon_corr = correlation_matrix(psi, "Nupdn", "Nupdn"; ishermitian=true)
        println("measure Hole-Doublon")
        @time hole_doublon_corr = correlation_matrix(psi, "hole", "Nupdn"; ishermitian=true)
        println("measure Total S.S correlations")
        @time SS = correlation_matrix(psi, "Sz", "Sz"; ishermitian=true) + 0.5 * (correlation_matrix(psi, "S+", "S-"; ishermitian=false) + correlation_matrix(psi, "S-", "S+"; ishermitian=false))

        push!(energies, energy)
        @printf("  Energy of METTS %d = %.4f\n", step, energy)
        a_E, err_E = avg_err(energies)
        @printf(
            "  Estimated Energy = %.4f +- %.4f  [%.4f,%.4f]\n",
            a_E,
            err_E,
            a_E - err_E,
            a_E + err_E
        )

        dump!(outfile, "energy", energy)
        dump!(outfile, "SZ", szs)
        dump!(outfile, "Ntotal", Ntotal)
        dump!(outfile, "DO", dob)
        dump!(outfile, "Hole", holes)
        dump!(outfile, "SZSZ", szszs)
        dump!(outfile, "SPSM", spsms)
        dump!(outfile, "NUPNUP", nupnup)
        dump!(outfile, "NDNNDN", ndnndn)
        dump!(outfile, "NUPNDN", nupndn)
        dump!(outfile, "NDNNUP", ndnnup)
        dump!(outfile, "NtotNtot", ntotntot)
        dump!(outfile, "SS", SS)
        dump!(outfile, "HoleHole", hole_corr)
        dump!(outfile, "HoleDoublon", hole_doublon_corr)
        dump!(outfile, "DoublonDoublon", doublon_doublon_corr)

        samp = collapse_with_qn!(psi, "X")
        @show samp

        # save sample
        samples[step] = samp

        save_samples(samples_filename, samples)

        new_state = [samp[j] == 1 ? "Emp" :
                     samp[j] == 2 ? "Up" :
                     samp[j] == 3 ? "Dn" : "UpDn" for j in 1:N]

        psi = randomMPS(sites, new_state)
        GC.gc()
    end

    return nothing
end

# Run the simulation with the parsed parameters
main(;L=L, W=W, U=U, filling=filling, t=t, T=T, tau=tau, maxD=maxD, cutoff=cutoff, seed=seed, NMETTS=nmetts)
Editor is loading...
Leave a Comment