Skip to content

Commit

Permalink
Add Gauss Hermite Quadrature Fields (#149)
Browse files Browse the repository at this point in the history
* add basic tests

* add attractive GHQ model

* add file IO

init! needs to happen before generate_groups because otherwise the stack size is based on an unassigned n_elements

* cleanup

* minor optimizations

* add fields as an indepent type

* cleanup

* add GHQ field

* use 1..4 as GHQ field values

* move and cleanup linalg

* add FieldCache & update matrix types

* add complex update code

* remove test model

* add vsubkron! and vsub!

* add vmul! with slices

* add vldiv22!

* generalize vmul! and vsub!

* remove redundant methods

* remove redundant accept_local!

* move linalg code

* cleanup

* update tests & fix global updates

* add linalg tests

* cleanup DQMC interface

* add DummyField

* update measurements

* merge Hubbard models

* fix typo

* cleanup

* pretend complex global updates work

Sorry future me

* add ED and update tests for odd fields
  • Loading branch information
ffreyer authored Jan 17, 2022
1 parent 9f748ee commit 592e824
Show file tree
Hide file tree
Showing 48 changed files with 2,175 additions and 1,423 deletions.
11 changes: 5 additions & 6 deletions src/MonteCarlo.jl
Original file line number Diff line number Diff line change
Expand Up @@ -59,10 +59,6 @@ const JLDFile = Union{FileWrapper{<: JLD.JldFile}, FileWrapper{<: JLD2.JLDFile},

include("helpers.jl")
export enable_benchmarks, disable_benchmarks, print_timer, reset_timer!
include("linalg/general.jl")
include("linalg/UDT.jl")
include("linalg/complex.jl") # TODO
include("linalg/blockdiagonal.jl")
include("flavors/abstract.jl")
include("models/abstract.jl")
include("lattices/abstract.jl")
Expand Down Expand Up @@ -106,7 +102,7 @@ export ReplicaExchange, ReplicaPull, connect, disconnect
# export mask, uniform_fourier, structure_factor, SymmetryWrapped, swave, eswave

include("models/Ising/IsingModel.jl")
include("models/HubbardModel/HubbardModel.jl")
include("models/HubbardModel.jl")
export IsingEnergyMeasurement, IsingMagnetizationMeasurement

include("FileIO.jl")
Expand All @@ -115,9 +111,12 @@ export save, load, resume!

export reset!
export run!, resume!, replay!
export Model, IsingModel, HubbardModel, HubbardModelAttractive, HubbardModelRepulsive
export Model, IsingModel
# export RepulsiveGHQHubbardModel, AttractiveGHQHubbardModel
export HubbardModel, HubbardModelAttractive, HubbardModelRepulsive
export MonteCarloFlavor, MC, DQMC
export greens, greens!, lattice, model, parameters
export DensityHirschField, MagneticHirschField, MagneticGHQField

# For extending
export AbstractMeasurement, Model
Expand Down
1 change: 1 addition & 0 deletions src/configurations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,7 @@ function Base.push!(cr::BufferedConfigRecorder, mc, model, sweep)
if (sweep % cr.rate == 0)
_push!(cr, compress(mc, model, conf(mc)))
end
nothing
end
function _push!(cr::BufferedConfigRecorder, data)
if cr.idx == -1
Expand Down
50 changes: 28 additions & 22 deletions src/flavors/DQMC/DQMC.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,33 +37,28 @@ function DQMC(model::M;
last_sweep = 0,
measure_rate = 10,
recording_rate = measure_rate,
recorder = ConfigRecorder(DQMC, M, recording_rate),
scheduler = SimpleScheduler(LocalSweep()),
field = choose_field(model),
recorder = ConfigRecorder(field, recording_rate),
kwargs...
) where M<:Model
# default params
# paramskwargs = filter(kw->kw[1] in fieldnames(DQMCParameters), kwargs)
parameters = DQMCParameters(measure_rate = measure_rate; kwargs...)

HET = hoppingeltype(DQMC, model)
GET = greenseltype(DQMC, model)
HMT = hopping_matrix_type(DQMC, model)
GMT = greens_matrix_type(DQMC, model)
IMT = interaction_matrix_type(DQMC, model)
stack = DQMCStack{GET, HET, GMT, HMT, IMT}()
ut_stack = UnequalTimeStack{GET, GMT}()

seed == -1 || Random.seed!(seed)
conf = rand(DQMC, model, parameters.slices)
field_data = field(parameters, model)
rand!(field_data)

stack = DQMCStack(field_data, model)
ut_stack = UnequalTimeStack{geltype(stack), gmattype(stack)}()

analysis = DQMCAnalysis()
CB = checkerboard ? CheckerboardTrue : CheckerboardFalse

mc = DQMC(
CB,
model, conf, deepcopy(conf), last_sweep,
stack, ut_stack, scheduler,
parameters, analysis,
recorder, thermalization_measurements, measurements
CB, model, field_data, last_sweep, stack, ut_stack, scheduler,
parameters, analysis, recorder, thermalization_measurements, measurements
)

init!(mc)
Expand All @@ -85,7 +80,8 @@ DQMC(m::Model, params::NamedTuple) = DQMC(m; params...)
@inline beta(mc::DQMC) = mc.parameters.beta
@inline nslices(mc::DQMC) = mc.parameters.slices
@inline model(mc::DQMC) = mc.model
@inline conf(mc::DQMC) = mc.conf
@inline field(mc::DQMC) = mc.field
@inline conf(mc::DQMC) = error("Replace this with field?")
@inline current_slice(mc::DQMC) = mc.stack.current_slice
@inline last_sweep(mc::DQMC) = mc.last_sweep
@inline configurations(mc::DQMC) = mc.recorder
Expand All @@ -101,7 +97,7 @@ import Base.show
Base.summary(mc::DQMC) = "DQMC simulation of $(summary(mc.model))"
function Base.show(io::IO, mc::DQMC)
print(io, "Determinant quantum Monte Carlo simulation\n")
print(io, "Model: ", mc.model, "\n")
print(io, "Model: ", summary(mc.model), "\n")
print(io, "Beta: ", mc.parameters.beta, " (T ≈ $(round(1/mc.parameters.beta, sigdigits=3)))\n")
N_th_meas = length(mc.thermalization_measurements)
N_me_meas = length(mc.measurements)
Expand Down Expand Up @@ -167,6 +163,7 @@ See also: [`resume!`](@ref)
total_sweeps = sweeps + thermalization

# Generate measurement groups
init!(mc)
th_groups = generate_groups(
mc, mc.model,
[mc.measurements[k] for k in keys(mc.thermalization_measurements) if !(k in ignore)]
Expand All @@ -183,10 +180,18 @@ See also: [`resume!`](@ref)

# fresh stack
verbose && println("Preparing Green's function stack")
init!(mc)
reverse_build_stack(mc, mc.stack)
propagate(mc)

# Check assumptions for global updates
copyto!(mc.stack.tmp2, mc.stack.greens)
udt_AVX_pivot!(mc.stack.tmp1, mc.stack.tempvf, mc.stack.tmp2, mc.stack.pivot, mc.stack.tempv)
ud = det(Matrix(mc.stack.tmp1))
td = det(Matrix(mc.stack.tmp2))
if !(0.9999999 <= abs(td) <= 1.0000001) || !(0.9999999 <= abs(ud) <= 1.0000001)
@error("Assumptions for global updates broken! ($td, $ud should be 1)")
end

min_sweeps = round(Int, 1 / min_update_rate)

_time = time() # for step estimations
Expand All @@ -212,7 +217,7 @@ See also: [`resume!`](@ref)
t0 = time()
end
else
push!(mc.recorder, mc, mc.model, mc.last_sweep)
push!(mc.recorder, field(mc), mc.last_sweep)
if iszero(mc.last_sweep % mc.parameters.measure_rate)
for (requirement, group) in groups
apply!(requirement, group, mc, mc.model, mc.last_sweep)
Expand Down Expand Up @@ -357,6 +362,7 @@ function replay!(
"There are no measurements set up for this simulation!"
)
# Generate measurement groups
init!(mc)
groups = generate_groups(
mc, mc.model,
[mc.measurements[k] for k in keys(mc.measurements) if !(k in ignore)]
Expand All @@ -368,17 +374,17 @@ function replay!(
end

verbose && println("Preparing Green's function stack")
init!(mc)
build_stack(mc, mc.stack)
propagate(mc)
mc.stack.current_slice = 1
mc.conf = rand(DQMC, mc.model, nslices(mc))
rand!(field(mc))

_time = time()
verbose && println("\n\nReplaying measurement stage - ", length(configurations))
prepare!(mc.measurements, mc, mc.model)
for i in mc.last_sweep+1:mc.parameters.measure_rate:length(configurations)
copyto!(mc.conf, decompress(mc, mc.model, configurations[i]))
# copyto!(mc.conf, decompress(mc, mc.model, configurations[i]))
decompress!(field(mc), configurations[i])
calculate_greens(mc, 0) # outputs to mc.stack.greens
for (requirement, group) in groups
apply!(requirement, group, mc, mc.model, i)
Expand Down
217 changes: 217 additions & 0 deletions src/flavors/DQMC/DQMC_interface.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,217 @@
################################################################################
### Mandatory Interface
################################################################################

### Model
########################################

"""
hopping_matrix(mc::DQMC, m::Model)
Calculates the hopping matrix \$T_{i\\sigma, j\\sigma '}\$ where \$i, j\$ are
site indices and \$\\sigma , \\sigma '\$ are flavor indices (e.g. spin indices).
The hopping matrix should also contain potential chemical potential terms on the
diagonal.
A matrix element is the hopping amplitude for a hopping process: \$j,\\sigma '
\\rightarrow i,\\sigma\$.
Regarding the order of indices, if `T[i, σ, j, σ']` is your desired 4D hopping
array, then `reshape(T, (n_sites * n_flavors, :))` is the hopping matrix.
"""
hopping_matrix(mc::DQMC, m::Model) = throw(MethodError(hopping_matrix, (mc, m)))

nflavors(m::Model) = throw(MethodError(nflavors, (m,)))

### Field
########################################

"""
interaction_matrix_exp!(
mc::DQMC, m::Model, field::AbstractField,
result::AbstractArray, slice::Int, power::Float64 = 1.0
)
Calculate the, exponentiated interaction matrix
`exp(- power * delta_tau * V(slice))` and store it in `result::AbstractArray`.
This only includes terms with 4 operators, i.e. not the chemical potential or
any hoppings. By default the calculation will be performed by the appropriate
field type (i.e. by `interaction_matrix_exp!(field, result, slice, power)`)
By default this function will call
`interaction_matrix_exp!(field, result, slice, power)`.
This is a performance critical method and one might consider efficient in-place
(in `result`) construction.
"""
@inline function interaction_matrix_exp!(
mc, model, field, result, slice, power = +1.0
)
interaction_matrix_exp!(field, result, slice, power)
end
function interaction_matrix_exp!(f, A, s, p)
throw(MethodError(interaction_matrix_exp!, (f, A, s, p)))
end

"""
propose_local(mc::DQMC, m::Model, field::AbstractField, i::Int, slice::Int)
Propose a local move for lattice site `i` at time slice `slice` for a `field`
holding the current configuration. Returns the Green's function determinant
ratio, the boson energy difference `ΔE_boson = E_boson_new - E_boson`,
and any extra information `passthrough` that might be useful in `accept_local`.
By default this function will call `propose_local(mc, field, i, slice)`.
See also [`accept_local!`](@ref).
"""
@inline propose_local(mc, m, field, i, slice) = propose_local(mc, field, i, slice)
propose_local(mc, f, i, s) = throw(MethodError(propose_local, (mc, f, i, s)))


"""
accept_local!(
mc::DQMC, m::Model, field::AbstractField, i::Int, slice::Int,
detratio, ΔE_boson, passthrough
)
Accept a local move for site `i` at imaginary time slice `slice` of current
configuration in `field`. Arguments `detratio`, `ΔE_boson` and `passthrough`
correspond to output of `propose_local` for that local move.
By default this function will call
`accept_local!(mc, field, i, slice, detration, ΔE_boson, passthrough)`
See also [`propose_local`](@ref).
"""
@inline function accept_local!(
mc, m, field, i, slice, detratio, ΔE_boson, passthrough
)
accept_local!(mc, field, i, slice, detratio, ΔE_boson, passthrough)
end
function accept_local!(mc, f, i, s, d, ΔE, pt)
throw(MethodError(accept_local!, (mc, f, i, s, d, ΔE, pt)))
end

Base.rand(f::AbstractField) = throw(MethodError(rand, (f,)))
Random.rand!(f::AbstractField) = throw(MethodError(rand!, (f,)))

nflavors(f::AbstractField) = throw(MethodError(nflavors, (f,)))

# For ConfigRecorder
compress(f::AbstractField) = throw(MethodError(compress, (f, )))
compressed_conf_type(f::AbstractField) = throw(MethodError(compressed_conf_type, (f, )))
decompress(f::AbstractField, c) = throw(MethodError(decompress, (f, c)))
decompress!(f::AbstractField, c) = throw(MethodError(decompress!, (f, c)))


################################################################################
### Optional Interface
################################################################################


### Mixed (Field + Model)
########################################

"""
greens_eltype(field, model)
Returns the type of the elements of the Green's function matrix. Defaults to
Float64 if both the hopping and interaction matrix contain floats and ComplexF64
otherwise.
"""
function greens_eltype(field::AbstractField, model::Model)
generalized_eltype(interaction_eltype(field), hopping_eltype(model))
end

"""
greens_matrix_type(field, model)
Returns the (matrix) type of the greens and most work matrices. Defaults to
`Matrix{greens_eltype(T, m)}`.
"""
greens_matrix_type(f::AbstractField, m::Model) = Matrix{greens_eltype(f, m)}


### Model
########################################

"""
hopping_eltype(model)
Returns the type of the elements of the hopping matrix. Defaults to `Float64`.
"""
hopping_eltype(::Model) = Float64

"""
hopping_matrix_type(field, model)
Returns the (matrix) type of the hopping matrix. Defaults to
`Matrix{hopping_eltype(model)}`.
"""
hopping_matrix_type(::AbstractField, m::Model) = Matrix{hopping_eltype(m)}


### Field
########################################

"""
interaction_eltype(model)
Returns the type of the elements of the interaction matrix. Defaults to `Float64`.
"""
interaction_eltype(::AbstractField) = Float64


"""
interaction_matrix_type(field, model)
Returns the (matrix) type of the interaction matrix. Defaults to
`Matrix{interaction_eltype(model)}`.
"""
interaction_matrix_type(f::AbstractField, m::Model) = Matrix{interaction_eltype(f)}


"""
init_interaction_matrix(field::AbstractField, model::Model)
Returns an initial interaction matrix. This only used to allocate a correctly
sized matrix.
By default this uses the matrix type from `interaction_matrix_type` and uses
`max(nflavors(field), nflavors(model)) * length(lattice(model))` as the size.
"""
function init_interaction_matrix(f::AbstractField, m::Model)
flv = max(nflavors(f), nflavors(m))
N = length(lattice(m))
FullT = interaction_matrix_type(f, m)

if FullT <: BlockDiagonal
MT = matrix_type(interaction_eltype(f))
return BlockDiagonal(MT(undef, N, N), MT(undef, N, N))
elseif FullT <: Diagonal
VT = vector_type(interaction_eltype(f))
Diagonal(VT(undef, flv * N))
else
FullT(undef, flv * N, flv * N)
end
end

"""
energy_boson(mc::DQMC, model::Model, conf)
Calculate bosonic part (non-Green's function determinant part) of energy for
configuration `conf` for Model `m`.
This is required for global and parallel updates as well as boson energy
measurements, but not for local updates. By default calls
`energy_boson(field(mc), conf)`
"""
energy_boson(mc::DQMC, m::Model, c = nothing) = energy_boson(field(mc), c)
energy_boson(f::AbstractField, c = nothing) = throw(MethodError(energy_boson, (f, c)))


conf(f::AbstractField) = f.conf
conf!(f::AbstractField, c) = conf(f) .= c
temp_conf(f::AbstractField) = f.temp_conf
Loading

0 comments on commit 592e824

Please sign in to comment.