From 177cc1e0605395294072ec609c6cb5486000b0ca Mon Sep 17 00:00:00 2001 From: Frederic Freyer Date: Sun, 20 Feb 2022 15:34:31 +0100 Subject: [PATCH] Update docs + cleanup (#161) * add field docs * add references * update introduction * simplify nflavors and hopping_matrix * update model docs * add choose_field * better type defaults * test pad_to_nflavors * update DQMC docs * update crosschecks * update HBC crosscheck * fix missing padding --- docs/make.jl | 1 + docs/src/DQMC/Introduction.md | 46 ++- docs/src/DQMC/dqmc.md | 7 +- docs/src/DQMC/fields.md | 62 ++++ docs/src/DQMC/model.md | 65 ++-- docs/src/examples/ALF1.md | 2 +- docs/src/examples/HBC_model.jl | 309 +++--------------- docs/src/examples/HofmannBergChowdhury.md | 245 +------------- docs/src/examples/triangular_Hubbard.jl | 2 +- docs/src/examples/triangular_Hubbard.md | 4 +- src/flavors/DQMC/DQMC_interface.jl | 67 +++- src/flavors/DQMC/fields.jl | 2 +- src/flavors/DQMC/main.jl | 7 +- src/flavors/DQMC/measurements/extensions.jl | 7 +- src/flavors/DQMC/measurements/measurements.jl | 4 +- src/flavors/DQMC/stack.jl | 12 +- src/models/DummyModel.jl | 3 +- src/models/HubbardModel.jl | 27 +- test/ED/ED_tests.jl | 2 +- test/fields.jl | 108 +++--- test/flavortests_DQMC.jl | 13 +- 21 files changed, 341 insertions(+), 654 deletions(-) create mode 100644 docs/src/DQMC/fields.md diff --git a/docs/make.jl b/docs/make.jl index 58090b7d..4186f7a1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,6 +10,7 @@ makedocs( "Introduction" => "DQMC/Introduction.md", "Model" => "DQMC/model.md", "Lattices" => "DQMC/lattice.md", + "Fields" => "DQMC/fields.md", "DQMC" => "DQMC/dqmc.md", "Configuration Recorder" => "DQMC/recorder.md", "Update Scheduler" => "DQMC/scheduler.md", diff --git a/docs/src/DQMC/Introduction.md b/docs/src/DQMC/Introduction.md index e61c0719..eee63582 100644 --- a/docs/src/DQMC/Introduction.md +++ b/docs/src/DQMC/Introduction.md @@ -1,20 +1,56 @@ # Determinant Quantum Monte Carlo -Determinant Quantum Monte Carlo is a Quantum Monte Carlo algorithm for fermionic Hamiltonians. The general idea is to use the Hubbard-Stranovich transformation to simplify the Hamiltonian to one with only quadratic fermionic degrees of freeedom. This introduces a bosonic fields for quartic terms, which can be sampled by a Monte Carlo procedure. +Determinant Quantum Monte Carlo is a Quantum Monte Carlo algorithm for fermionic Hamiltonians. The general idea is to use the Hubbard-Stranovich transformation to simplify any term with four fermionic operators to two. This comes at the cost of introducing a bosonic field which is sampled by a Monte Carlo procedure. The minimal working example for a DQMC simulation is the following. ```julia using MonteCarlo -model = HubbardModelAttractive(4, 2) -dqmc = DQMC(model, beta=1.0) +dqmc = DQMC(HubbardModel(), beta = 1.0) run(dqmc) ``` -This will set up and run a DQMC simulation at inverse temperature $\beta = 1.0$ using an attractive Hubbard model with a four by four square lattice. Note that by default no measurements are taken. +This will set up and run a DQMC simulation at inverse temperature $\beta = 1.0$ using an attractive Hubbard model with a two by two square lattice. Of course this example is not very useful. In practice you will want to modify lots of parameters and components of the simulation. We give a brief overview of the different components here. You can also check the examples section for some more involved and realistic examples. + +#### Model + +The model is a description of the Hamiltonian. Currently MonteCarlo.jl only implements the `HubbardModel` which includes nearest neighbor hoppings, a chemical potential and the Hubbard interaction. One aim of this package is to make it simple to implement different models, though currently only Hubbard interactions are supported. + +#### Lattice + +The lattice is interchangeable part of the model. We currently support a few lattices natively (Chain, Square, Cubic, Triangular) and provide a Wrapper for [LatticePhysics.jl](https://github.com/janattig/LatticePhysics.jl) + +#### Field + +The field defines how the interaction is handled, i.e. what Monte Carlo configurations look like, how the interactive term is interpreted and how updates work. We currently implement Hirsch fields and Gauß-Hermite Quadrature fields, both with a density- and magnetic-channel version. Any of these choices works with any choice of the Hubbard model, though performance and accuracy can vary. + +#### Updates + +Beyond local updates we provide a selection of global and experimental parallel updates. Like models, this part is designed to be relatively easy to extend. Most global updates just generate a new configuration and call another function to accept or reject it. + +The updates themselves are collected in a scheduler. We currently provide two types of them - one which simply iterates through updates and one which adaptively favors updates with high acceptance rates. + +#### Configuration Recorder + +We provide an interface to record configurations during a MonteCarlo run so that they can later be replayed. This interface also implements functionality for compressing configurations. + +#### Measurements + +We currently include measurements of the following observables: +- Equal time Greens function (and to a less direct extend time displaced Greens functions) +- occupations +- charge density correlations and susceptibilities +- magnetization in x, y and z direction +- spin density correlations and susceptibilities in x, y and z direction +- pairing correlation and susceptibility for generic symmetries +- current current susceptibilities which are required for the superfluid stiffness (this can also be calculated more directly) +- total and partial energies (quadratic terms and quartic terms separated) + +#### DQMC + +The `DQMC` struct represents the full simulation. All the components above are collected here. It also contains a bunch paramters, such as the inverse temperature `beta`, the time step discretization `delta_tau`, the number of thermazition and measurements sweeps and more. -In the following pages we will discuss the various components that go into a DQMC simulation. ## Derivation diff --git a/docs/src/DQMC/dqmc.md b/docs/src/DQMC/dqmc.md index 0d54d8d0..35e7d409 100644 --- a/docs/src/DQMC/dqmc.md +++ b/docs/src/DQMC/dqmc.md @@ -23,6 +23,7 @@ Beyond this there are a couple of keyword arguments which are more involved and * `thermalization_measurements = Dict{Symbol, AbstractMeasurement}()`: Same as the above, but the measurements run during the thermalization stage. Might be useful to judge convergence or the number of necessary sweeps. * `recorder = ConfigRecorder`: A recorder for the configurations generated during the simulation. * `recording_rate = measure_rate`: The rate at which configurations are recorded. +* `field::Type = choose_field(model)`: The type of the field used by the simulation Running a simulation also comes with a bunch of options via keyword arguments - most dealing with saving the simulation. The options for `run(dqmc)` include: @@ -31,5 +32,7 @@ Running a simulation also comes with a bunch of options via keyword arguments - * `safe_before::TimeType = now() + Year(100)`: Sets a time stamp before which the simulation will cancel itself and save. * `safe_every::TimePeriod = Hour(10000)`: Sets a time period interval for regular saves. * `grace_period::TimePeriod = Minute(5)`: Sets a buffer time period for saving. By default saves at least 5 minutes before the requested time. -* `resumable_filename = "resumable_$(datestring).jld2`: Sets the name the savefile generated from safe_before and safe_every. -* `overwrite = true`: Enables or disables overwriting of existing files. \ No newline at end of file +* `resumable_filename = "resumable_$datestring.jld2`: Sets the name the savefile generated from safe_before and safe_every. +* `overwrite = true`: Enables or disables overwriting of existing files. +* `min_update_rate = 0.001`: Sets a minimum update acceptance rate under which a simulation is considered "failed". This will only trigger after `1 / min_update_rate` sweeps and only if all updates fall below this rate. +* `fail_filename = "failed_$datestring.jld2"`: Sets the filename of a failed simulation. \ No newline at end of file diff --git a/docs/src/DQMC/fields.md b/docs/src/DQMC/fields.md new file mode 100644 index 00000000..94d75a40 --- /dev/null +++ b/docs/src/DQMC/fields.md @@ -0,0 +1,62 @@ +# Fields + +In MonteCarlo.jl fields describe how to handle the interactive term(s) in a model/Hamiltonian. Currently we implement two transformations, the Hirsch transformation and a Gauß-Hermite Quadrature based approximation of the Hubbard Stratonovich transformation, each with a density- and spin-channel variant. + +## Hirsch Transformation + +The Hirsch transformation applies specifically to a Hubbard interaction $U (n_{i\uparrow} - 0.5) (n_{i\downarrow} - 0.5)`. Each number operator can take the values 0 or 1, which allows us to rewrite the term as either $0.5 U (n_{i\uparrow} + n_{i\downarrow} - 1)^2 - 0.25$ or $0.5 U (n_{i\uparrow} - n_{i\downarrow})^2 + 0.25$. These two cases are the density- and spin-channel respectively. The next step is to is to introduce a bosonic field $x$ such that + +```math +e^{|U|\Delta \tau \left( n_{i\uparrow} - \frac{1}{2} \right) \left(n_{i\downarrow} - \frac{1}{2} \right)} += \frac{1}{2} e^{-\frac{1}{4} \Delta \tau U} \sum_{x=\pm 1} e^{\alpha x (n_{i\uparrow} + n_{i\downarrow} - \frac{1}{2})} +``` + +or + +```math +e^{|U|\Delta \tau \left( n_{i\uparrow} - \frac{1}{2} \right) \left(n_{i\downarrow} - \frac{1}{2} \right)} += \frac{1}{2} e^{-\frac{1}{4} \Delta \tau U} \sum_{x=\pm 1} e^{\alpha x (n_{i\uparrow} - n_{i\downarrow})}. +``` + +For both of these cases we define a positive $U$ as the attractive case. The constant $\alpha = acosh(exp(0.5 \Delta\tau U))$ in the attractive case and $\alpha = acosh(exp(-0.5 \Delta\tau U))$ in the repulsive case. Both of these can be verified by expanding the sum over $x$ and going through the different possible values of the number operator $n$. + +In MonteCarlo.jl the first case is implemented as `DensityHirschField` and the second as `MagneticHirschField`. Both can be used with positive and negative $U$, though $\alpha$ may become complex. This will of course influence performance. Furthermore the density channel is symmetric between spin up and down, which is used to reduced the size of matrices involved. + +See [Quantum Monte Carlo Methods](https://doi.org/10.1017/CBO9780511902581) + +## Gauß-Hermite Quadrature + +In this case we use Gaussian quadrature to approximate the integral generated by the Hubbard Stratonovich transformation. This is a more general approach than the Hirsch transformation, though we currently only have an implementation for Hubbard interactions. + +We start of the same way as with the Hirsch transformation - by rewriting the our interaction as a squared term. Generally speaking we need to arrive at $exp(-0.5\Delta\tau (c_i^\dagger V_{ij} c_j + C)^2 + C^\prime)$. Using the Hubbard Stratonovich transformation this becomes + +```math +e^{-\frac{1}{2}\Delta\tau (c_i^\dagger V_{ij} c_j + C)^2 + C^\prime} = e^{-\Delta\tau C^\prime} \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{1}{2}x^2 \pm x \sqrt{\Delta\tau} (c_i^\dagger V_{ij} c_j + C)} dx +``` + +One method of solving this integral is Gauß-Hermite quadrature. With four weights and notes we get + +```math +\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{1}{2}x^2 \pm x \sqrt{\Delta\tau} (c_i^\dagger V_{ij} c_j + C)} dx = \đrac{1}{4} \sum_{x = \pm 1, \pm 2} \gamma(x) e^{\sqrt{\Delta\tau} \eta(x) (c_i^\dagger V_{ij} c_j^\dagger + C)} +``` + +where $\gamma$ and $\eta$ give the weights and nodes for a given $x$. For the density-channel version $V_\uparrow = V_\downarrow = \sqrt{U}$ and $C = 0.5 \sqrt{U}$, for the spin channel $V_\uparrow = - V\downarrow = \sqrt{U}$ and $C = 0$. These two cases are implemented as `DensityGHQField` and `MagneticGHQField`. Just like with the Hirsch field these two methods can create real or complex interaction matrices depending on the sign $U$ and the Density case makes use of symmetry. + +See [ALF](https://git.physik.uni-wuerzburg.de/ALF/ALF) Documentation (Auxiliary Field Quantum Monte Carlo section) and https://arxiv.org/pdf/2009.04491.pdf + +## Performance & Accuracy + +As a rough estimate of relative performance we can look at the runtime of a small Hubbard model (6x6 square lattice, $\beta = 1$, 10 000 sweeps total): + +| Field Type | U = +1 | U = -1 | +| ------------------- | ------ | ------ | +| DensityHirschField | 6.8s | 23s | +| MagneticHirschField | 51s | 17s | +| DensityGHQField | 6.2s | 22s | +| MagneticGHQField | 46s | 15s | + +As you can see the runtimes are relatively stable for the repulsive Hubbard model (negative U). This is because we either have to include spin up and down channels or work with complex matrices. In the attractive case both of these fall together. With a density-channel field we make use of spin symmetry and work with real matrices, but with a magnetic-channel field both of these slow us down dramatically. + +In terms of accuracy none of these should technically have a big advantage over the other other in a Hubbard model. However this doesn't generally seem to be true in practice. In one model we had issues resolving occupations accurately with `DensityHirschField` and found that spin-channel fields do a better job of resolving most observables. The exception here are spin densities, which are resolved more cleanly by density-channel fields. Overall the `MagneticGHQField` gave the best results. + +In general it is probably a good idea to simulate a smaller system with all the different field options to see which yields the best combination of performance and accuracy. \ No newline at end of file diff --git a/docs/src/DQMC/model.md b/docs/src/DQMC/model.md index 0c0d578d..f8b832f4 100644 --- a/docs/src/DQMC/model.md +++ b/docs/src/DQMC/model.md @@ -1,53 +1,54 @@ # Models -The `Model` is an abstract type whose concrete versions should implement lightweight description of a Hamiltonian. There are two Models implemented for DQMC - the attractive and repulsive Hubbard model. +In MonteCarlo.jl a `Model` describes a Hamiltonian. It's primary purpose is to collect paramaters from different terms in the Hamiltonian as well as the lattice and generate a hopping matrix for the simulation. We currently provide one model - the Hubbard model. +## Hubbard Model -## Attractive Hubbard Model - -The attractive (negative interaction term) Hubbard model is given by +The Hubbard model is given by ```math -\mathcal{H} = -t \sum_{\langle i,j \rangle, \sigma} \left( c^\dagger_{i\sigma} c_{j\sigma} + \text{h.c.} \right) - |U| \sum_j \left( n_{j\uparrow} - \frac{1}{2} \right) \left( n_{j\downarrow} - \frac{1}{2} \right) - \mu\sum_j n_{j}, +\mathcal{H} = -t \sum_{\langle i,j \rangle, \sigma} \left( c^\dagger_{i\sigma} c_{j\sigma} + \text{h.c.} \right) - U \sum_j \left( n_{j\uparrow} - \frac{1}{2} \right) \left( n_{j\downarrow} - \frac{1}{2} \right) - \mu\sum_j n_{j}, ``` -where $\sigma$ denotes spin, $t$ is the hopping amplitude, $U$ the on-site attractive interaction strength, $\mu$ the chemical potential and $\langle i, j \rangle$ indicates that the sum has to be taken over nearest neighbors. Note that (1) is written in particle-hole symmetric form such that $\mu = 0$ corresponds to half-filling. +where $\sigma$ denotes spin, $t$ is the hopping amplitude, $U$ the on-site Hubbard interaction strength, $\mu$ the chemical potential and $\langle i, j \rangle$ indicates that the sum has to be taken over nearest neighbors. Note that the Hamiltonian is written in particle-hole symmetric form such that $\mu = 0$ corresponds to half-filling. -The parameters $U = 1.0$, $t = 1.0$ and $\mu = 0$ as well as the lattice are saved in the model. They can be specified via keyword arguments when calling `HubbardModelAttractive()`. +Our implementation allows for both attractive (positive) and repulsive (negative) $U$. Note that for the repulsive case there is a sign problem for $\mu \ne 0$. The model also works with any lattice, assuming that lattice provides the required functionality. -## Repulsive Hubbard Model +You can create a Hubbard model with `HubbardModel()`. Optional keyword arguments include: +- `l::AbstractLattice = choose_lattice(HubbardModel, dims, L)` is the lattice used by the model. +- `dims::Integer = 2` is the dimensionality of the default lattice (Chain/Square/Cubic) +- `L::Integer = 2` is the linear system size of the default lattice (Chain/Square/Cubic) +- `t::Float64 = 1.0` is the hopping strength. +- `mu::Float64 = 0.0` is the chemical potential. (Must be 0 if U is negative.) +- `U::Float64 = 1.0` is the interaction strength. -The repulsive (positive interaction term) Hubbard model is given by +## Creating your own Model -```math -\mathcal{H} = -t \sum_{\langle i,j \rangle, \sigma} \left( c^\dagger_{i\sigma} c_{j\sigma} + \text{h.c.} \right) + |U| \sum_j \left( n_{j\uparrow} - \frac{1}{2} \right) \left( n_{j\downarrow} - \frac{1}{2} \right), -``` +To create your own model you will need to inherit from the abstract type `Model`. There is a set of mandatory and optional methods you must/can implement: + +#### Mandatory Methods -where $\sigma$ denotes spin, $t$ is the hopping amplitude, $U$ the on-site repulsive interaction strength and $\langle i, j \rangle$ indicates that the sum has to be taken over nearest neighbors. This model does not include a chemical potential due to the sign problem. +A custom model needs to implement these methods to function -The parameters $U = 1.0$ and $t = 1.0$ as well as the lattice are saved in the model. They can be specified via keyword arguments when calling `HubbardModelRepulsive()`. +- `lattice(model)` needs to return a MonteCarlo compatible lattice +- `nflavors(model)` needs to return the number of unique fermion flavors of the hopping matrix. For example, in a two spin model this would return 2 if the hopping matrix is different between spin up and down, or 1 if one sector is a copy of the other. Internally this is used together with `nflavors(field)` to optimize spin/flavor symmetric systems. +- `hopping_matrix(model)` needs to generate the hopping matrix, which includes all quadratic terms. (I.e. also the chemical potential.) The hopping matrix should only include as many flavors as necessary. If the hopping matrix contains two copies of the same matrix, one for spin up and one for spin down for example, then it should only return one of these. Expanding it to an appropriate size is handled internally. -## Creating your own Model +#### (Semi-) Optional Methods -To create your own model you will need to inherit from `Model`. There is a set of mandatory and optional methods you must/can implement: +These methods aren't strictly necessary to implement, but may boost performance when implemented. It is recommended to provide these if the defaults do not apply. -#### Mandatory Methods -* `rand(::Type{DQMC}, m::Model, nslices)` should return a full new random configuration, where `nslices` is the number of imaginary time slices. -* `nflavors(m::Model)` should return the number of fermion flavors used by the model. This controls the size of most simulation matrices, e.g. the greens matrix, hopping and interaction matrices, etc. -* `hopping_matrix(dqmc::DQMC, m::Model)` should return the hopping matrix of the model. This includes all terms that are not connected to a bosonic field. -* `interaction_matrix_exp!(dqmc::DQMC, m::Model, result::AbstractArray, conf, slice, power=1.0)` should calculate $exp(- power \cdot delta_tau \cdot V_{slice})$ for the current `conf` at time slice `slice` and save it to `result`. Note that this method is performance critical. -* `propose_local(dqmc::DQMC; m::Model, i, slice, conf)` should propose a local update at site `i` and the current tiem slice. It should calculate the determinant ratio and bosonic energy difference (which maybe 0 for some models) and return `determinant_ratio, bosonic_energy, passthrough` where `passthrough` are a tuple of variables you may want to use in `accept_local!`. Note this is also a performance critical method. -* `accept_local!(mc::DQMC, m::Model, i, slice, conf, detratio, ΔE_boson, passthrough)` should update the greens function when the proposed change is accepted. Specifically this means updating `dqmc.stack.greens`. This is a performance critical method. +- `save_model(file::JLDFile, model, entryname)` should write model information to the given file. It should also save the lattice via `save_lattice` and save a unqiue `tag`. If this is not implemented JLD2 will be asked to save the type as is, which makes it hard to load data when the model type is edited. +- `_load_model(data, ::Val{Symbol(tag)})` loads a model from `data`, which typically is a JLDFile. Note that saved tag is used to dispatch to the correct method. +- `intE_kernel(mc, model, G, ::Val{flv})` should be implemented to enable measurements of the energy from the interactive term as well as the total energy. +- `choose_field(model)` returns the default field type for a given model. If this is not implemented a field must be passed to `DQMC`. -#### (Semi-) Optional Methods +#### Optional Methods -* `greenseltype(::Type{DQMC}, m::Model) = ComplexF64` sets the element type of the greens function. If your greens function is real this function should be implemented for better performance. -* `hoppingeltype(::Type{DQMC}, m::Modle) = Float64` sets the expected element type of the hopping matrix. Adjust this if it is wrong. -* `interaction_matrix_type(::Type{DQMC}, m::Model) = Matrix{greenseltype(DQMC, m)}` Sets the matrix type of the interaction matrix. For the Hubbard model this is `Diagonal{Float64}`. -* `hopping_matrix_type(::Type{DQMC}, m::Model) = Matrix{hoppingeltype(DQMC, m)}` Sets the matrix type of the hopping matrix. If there is a more efficient way to represent your matrix you may use this to change to that representation. The repulsive Hubbard model uses `BlockDiagonal` matrices, for example, which are implemented by MonteCarlo.jl -* `greens_matrix_type(::Type{DQMC}, m::Model) = Matrix{greenseltype(DQMC, m)}` Sets the matrix type of the greens matrix. The repulsive Hubbard model also uses `BlockDiagonal` here. -* `init_interaction_matrix(m::Model)` initializes the interaction matrix. If you are using a custom matrix type you must implement this. Note that the values are irrelevant. -* `energy_boson(dqmc::DQMC, m::Model, conf)` should be implemented if you want to measure energies in your simulation. -* `parameters(m::Model)` should collect the parameters from the model and lattice in a NamedTuple. This is purely utility. +- `hopping_eltype(model) = eltype(hopping_matrix(model))` returns the element type of the hopping matrix. +- `hopping_matrix_type(field, model) = typeof(pad_to_nflavors(field, model, hopping_matrix(model)))` return the matrix type of the hopping matrix. +- `greens_eltype(field, model) = generalized_eltype(interaction_eltype(field), hopping_eltype(model))` returns the element type of the greens function. This must be compatible with element types of the hopping matrix and the interaction matrix. +- `greens_matrix_type(field, model) = Matrix{greens_eltype(field, model)}` returns the full type of the greens function. This must be compatible with both the type of the hopping matrix and the interaction matrix. +* `parameters(m::Model)` should collect the parameters from the model and lattice in a NamedTuple. Also note that you may need to update the measurement kernels. More information about that on the measurement page. \ No newline at end of file diff --git a/docs/src/examples/ALF1.md b/docs/src/examples/ALF1.md index 4b61b69b..a4a285a2 100644 --- a/docs/src/examples/ALF1.md +++ b/docs/src/examples/ALF1.md @@ -99,7 +99,7 @@ using Revise, MonteCarlo, Printf, LinearAlgebra mcs = [] @time for beta in [1.0, 6.0, 12.0] - m = HubbardModelRepulsive(4, 2, U = 4) + m = HubbardModel(4, 2, U = -4) mc = DQMC( m, beta=beta, thermalization=5_000, sweeps=15_000, print_rate=5_000, delta_tau = 0.05#, measure_rate=5 diff --git a/docs/src/examples/HBC_model.jl b/docs/src/examples/HBC_model.jl index b2e8ed0c..50c33db9 100644 --- a/docs/src/examples/HBC_model.jl +++ b/docs/src/examples/HBC_model.jl @@ -74,38 +74,26 @@ end -using MonteCarlo -using MonteCarlo: StructArray, CMat64, CVec64 -using MonteCarlo: AbstractLattice, @bm, HubbardConf, @turbo, BlockDiagonal, vmul!, conf - +using MonteCarlo, LinearAlgebra +using MonteCarlo: StructArray, CMat64, CVec64, AbstractField +using MonteCarlo: AbstractLattice, BlockDiagonal, vmul!, conf """ t5 = 0 <=> F = 0.2 t5 = (1 - sqrt(2)) / 4 <=> F = 0.009 """ -MonteCarlo.@with_kw_noshow struct HBCModel{LT<:AbstractLattice} <: HubbardModel - # parameters with defaults based on paper +MonteCarlo.@with_kw_noshow struct HBCModel{LT<:AbstractLattice} <: MonteCarlo.Model + # user optional mu::Float64 = 0.0 U::Float64 = 1.0 @assert U >= 0. "U must be positive." t1::Float64 = 1.0 t2::Float64 = 1.0 / sqrt(2.0) t5::Float64 = (1 - sqrt(2)) / 4 - - # lattice l::LT - - # two fermion flavors (up, down) - flv::Int = 2 - - # temp storage to avoid allocations in propose_local and accept_local - IG::CMat64 = StructArray(Matrix{ComplexF64}(undef, length(l), 2)) - IGR::CMat64 = StructArray(Matrix{ComplexF64}(undef, length(l), 2)) - R::Diagonal{ComplexF64, CVec64} = Diagonal(StructArray(Vector{ComplexF64}(undef, 2))) end -# Constructors HBCModel(params::Dict{Symbol}) = HBCModel(; params...) HBCModel(params::NamedTuple) = HBCModel(; params...) function HBCModel(lattice::AbstractLattice; kwargs...) @@ -116,13 +104,14 @@ function HBCModel(L, dims; kwargs...) HBCModel(l = l; kwargs...) end +MonteCarlo.nflavors(::HBCModel) = 2 +MonteCarlo.lattice(m::HBCModel) = m.l +MonteCarlo.choose_field(m::HBCModel) = m.U < 0.0 ? MagneticHirschField : DensityHirschField -# type information for DQMC stack -MonteCarlo.hoppingeltype(::Type{DQMC}, ::HBCModel) = ComplexF64 -MonteCarlo.hopping_matrix_type(::Type{DQMC}, ::HBCModel) = BlockDiagonal{ComplexF64, 2, CMat64} -MonteCarlo.greenseltype(::Type{DQMC}, ::HBCModel) = ComplexF64 -MonteCarlo.greens_matrix_type( ::Type{DQMC}, ::HBCModel) = BlockDiagonal{ComplexF64, 2, CMat64} - +MonteCarlo.hopping_eltype(::HBCModel) = ComplexF64 +MonteCarlo.hopping_matrix_type(::AbstractField, ::HBCModel) = BlockDiagonal{ComplexF64, 2, CMat64} +MonteCarlo.greens_eltype(::AbstractField, ::HBCModel) = ComplexF64 +MonteCarlo.greens_matrix_type( ::AbstractField, ::HBCModel) = BlockDiagonal{ComplexF64, 2, CMat64} # cosmetics import Base.summary @@ -138,256 +127,43 @@ Base.show(io::IO, m::MIME"text/plain", model::HBCModel) = print(io, model) @inline MonteCarlo.parameters(m::HBCModel) = (N = length(m.l), t1 = m.t1, t2 = m.t2, t5 = m.t5, U = -m.U, mu = m.mu) -# generate hopping matrix -function MonteCarlo.hopping_matrix(mc::DQMC, m::HBCModel{<: LatPhysLattice}) - # number of sites +function MonteCarlo.hopping_matrix(m::HBCModel{<: LatPhysLattice}) N = length(m.l) + t1 = diagm(0 => fill(-ComplexF64(m.mu), N)) + t2 = diagm(0 => fill(-ComplexF64(m.mu), N)) - # spin up and spin down blocks of T - tup = diagm(0 => fill(-ComplexF64(m.mu), N)) - tdown = diagm(0 => fill(-ComplexF64(m.mu), N)) - - # positive and negative prefactors for t1, t2 - t1p = m.t1 * cis(+pi/4) # ϕ_ij^↑ = + π/4 - t1m = m.t1 * cis(-pi/4) # ϕ_ij^↓ = - π/4 + t1p = m.t1 * cis(+pi/4) + t1m = m.t1 * cis(-pi/4) t2p = + m.t2 t2m = - m.t2 + # Nearest neighbor hoppings for b in bonds(m.l.lattice) - # NN paper direction if b.label == 1 - tup[b.from, b.to] = - t1p - tdown[b.from, b.to] = - t1m - - # NN reverse direction + # ϕ_ij^↑ = + π/4 + t1[b.from, b.to] = - t1p + # ϕ_ij^↓ = - π/4 + t2[b.from, b.to] = - t1m elseif b.label == 2 - tup[b.from, b.to] = - t1m - tdown[b.from, b.to] = - t1p + # TODO do we use reverse NN? - doesn't look like it (sign problem) + t1[b.from, b.to] = - t1m + t2[b.from, b.to] = - t1p - # NNN solid bonds elseif b.label == 3 - tup[b.from, b.to] = - t2p - tdown[b.from, b.to] = - t2p - - # NNN dashed bonds + t1[b.from, b.to] = - t2p + t2[b.from, b.to] = - t2p elseif b.label == 4 - tup[b.from, b.to] = - t2m - tdown[b.from, b.to] = - t2m - - # Fifth nearest neighbors + t1[b.from, b.to] = - t2m + t2[b.from, b.to] = - t2m else - tup[b.from, b.to] = - m.t5 - tdown[b.from, b.to] = - m.t5 - end - end - - return BlockDiagonal(StructArray(tup), StructArray(tdown)) -end - - -""" -Calculate the interaction matrix exponential `expV = exp(- power * delta_tau * V(slice))` -and store it in `result::Matrix`. - -This is a performance critical method. -""" -@inline @bm function MonteCarlo.interaction_matrix_exp!(mc::DQMC, m::HBCModel, - result::Diagonal, conf::HubbardConf, slice::Int, power::Float64=1.) - dtau = mc.parameters.delta_tau - lambda = acosh(exp(0.5 * m.U * dtau)) - N = length(lattice(m)) - - # spin up block - @inbounds for i in 1:N - result.diag[i] = exp(sign(power) * lambda * conf[i, slice]) - end - - # spin down block - @inbounds for i in 1:N - result.diag[N+i] = exp(sign(power) * lambda * conf[i, slice]) - end - nothing -end - - - -@inline @bm function MonteCarlo.propose_local( - mc::DQMC, model::HBModel, i::Int, slice::Int, conf::HubbardConf - ) - N = length(model.l) - G = mc.stack.greens - Δτ = mc.parameters.delta_tau - R = model.R - - α = acosh(exp(0.5Δτ * model.U)) - ΔE_boson = -2.0α * conf[i, slice] - Δ = exp(ΔE_boson) - 1.0 - - # Unrolled R = I + Δ * (I - G) - # up-up term - R.diag.re[1] = 1.0 + Δ * (1.0 - G.blocks[1].re[i, i]) - R.diag.im[1] = - Δ * G.blocks[1].im[i, i] - # down-down term - R.diag.re[2] = 1.0 + Δ * (1.0 - G.blocks[2].re[i, i]) - R.diag.im[2] = - Δ * G.blocks[2].im[i, i] - - # Calculate "determinant" - detratio = ComplexF64( - R.diag.re[1] * R.diag.re[2] - R.diag.im[1] * R.diag.im[2], - R.diag.re[1] * R.diag.im[2] + R.diag.im[1] * R.diag.re[2] - ) - - return detratio, ΔE_boson, Δ -end - -@inline @bm function MonteCarlo.accept_local!( - mc::DQMC, model::HBModel, i::Int, slice::Int, conf::HubbardConf, - detratio, ΔE_boson, Δ) - - @bm "accept_local (init)" begin - N = length(model.l) - G = mc.stack.greens - IG = model.IG - IGR = model.IGR - R = model.R - end - - # compute R⁻¹ Δ, using that R is Diagonal, Δ is Number - # using Δ / (a + ib) = Δ / (a^2 + b^2) * (a - ib) - @bm "accept_local (inversion)" begin - f = Δ / (R.diag.re[1]^2 + R.diag.im[1]^2) - R.diag.re[1] = +f * R.diag.re[1] - R.diag.im[1] = -f * R.diag.im[1] - f = Δ / (R.diag.re[2]^2 + R.diag.im[2]^2) - R.diag.re[2] = +f * R.diag.re[2] - R.diag.im[2] = -f * R.diag.im[2] - end - - # Compute (I - G) R^-1 Δ - # Note IG is reduced to non-zero entries. Full IG would be - # (I-G)[:, i] 0 - # 0 (I-G)[:, i+N] - # our IG is [(I-G)[:, i] (I-G)[:, i+N]] - @bm "accept_local (IG, R)" begin - # Calculate IG = I - G (relevant entries only) - @turbo for m in axes(IG, 1) - IG.re[m, 1] = -G.blocks[1].re[m, i] - end - @turbo for m in axes(IG, 1) - IG.re[m, 2] = -G.blocks[2].re[m, i] - end - @turbo for m in axes(IG, 1) - IG.im[m, 1] = -G.blocks[1].im[m, i] - end - @turbo for m in axes(IG, 1) - IG.im[m, 2] = -G.blocks[2].im[m, i] - end - IG.re[i, 1] += 1.0 - IG.re[i, 2] += 1.0 - - # Calculate IGR = IG * R where R = R⁻¹ Δ from the - # previous calculation (relevant entries only) - # spin up-up block - @turbo for m in axes(IG, 1) - IGR.re[m, 1] = IG.re[m, 1] * R.diag.re[1] - end - @turbo for m in axes(IG, 1) - IGR.re[m, 1] -= IG.im[m, 1] * R.diag.im[1] - end - @turbo for m in axes(IG, 1) - IGR.im[m, 1] = IG.re[m, 1] * R.diag.im[1] - end - @turbo for m in axes(IG, 1) - IGR.im[m, 1] += IG.im[m, 1] * R.diag.re[1] - end - - # spin down-down block - @turbo for m in axes(IG, 1) - IGR.re[m, 2] = IG.re[m, 2] * R.diag.re[2] - end - @turbo for m in axes(IG, 1) - IGR.re[m, 2] -= IG.im[m, 2] * R.diag.im[2] - end - @turbo for m in axes(IG, 1) - IGR.im[m, 2] = IG.re[m, 2] * R.diag.im[2] - end - @turbo for m in axes(IG, 1) - IGR.im[m, 2] += IG.im[m, 2] * R.diag.re[2] + t1[b.from, b.to] = - m.t5 + t2[b.from, b.to] = - m.t5 end end - # Update G according to G = G - (I - G)[:, i:N:end] * R⁻¹ * Δ * G[i:N:end, :] - # We already have IG = (I - G)[:, i:N:end] * R⁻¹ * Δ - @bm "accept_local (finalize computation)" begin - # get blocks to write less - G1 = G.blocks[1] - G2 = G.blocks[2] - temp1 = mc.stack.greens_temp.blocks[1] - temp2 = mc.stack.greens_temp.blocks[2] - - # compute temp = IG[:, i:N:end] * G[i:N:end, :] - # spin up-up block - @turbo for m in axes(G1, 1), n in axes(G1, 2) - temp1.re[m, n] = IGR.re[m, 1] * G1.re[i, n] - end - @turbo for m in axes(G1, 1), n in axes(G1, 2) - temp1.re[m, n] -= IGR.im[m, 1] * G1.im[i, n] - end - @turbo for m in axes(G1, 1), n in axes(G1, 2) - temp1.im[m, n] = IGR.im[m, 1] * G1.re[i, n] - end - @turbo for m in axes(G1, 1), n in axes(G1, 2) - temp1.im[m, n] += IGR.re[m, 1] * G1.im[i, n] - end - - # spin down-down block - @turbo for m in axes(G2, 1), n in axes(G2, 2) - temp2.re[m, n] = IGR.re[m, 2] * G2.re[i, n] - end - @turbo for m in axes(G2, 1), n in axes(G2, 2) - temp2.re[m, n] -= IGR.im[m, 2] * G2.im[i, n] - end - @turbo for m in axes(G2, 1), n in axes(G2, 2) - temp2.im[m, n] = IGR.im[m, 2] * G2.re[i, n] - end - @turbo for m in axes(G2, 1), n in axes(G2, 2) - temp2.im[m, n] += IGR.re[m, 2] * G2.im[i, n] - end - - # Calculate G = G - temp - # spin up-up block - @turbo for m in axes(G1, 1), n in axes(G1, 2) - G1.re[m, n] = G1.re[m, n] - temp1.re[m, n] - end - @turbo for m in axes(G1, 1), n in axes(G1, 2) - G1.im[m, n] = G1.im[m, n] - temp1.im[m, n] - end - - # spin down-down block - @turbo for m in axes(G2, 1), n in axes(G2, 2) - G2.re[m, n] = G2.re[m, n] - temp2.re[m, n] - end - @turbo for m in axes(G2, 1), n in axes(G2, 2) - G2.im[m, n] = G2.im[m, n] - temp2.im[m, n] - end - - # Update configuration - conf[i, slice] *= -1 - end - - nothing + return BlockDiagonal(StructArray(t1), StructArray(t2)) end - -# this enables boson energy measurements and global updates -@inline function MonteCarlo.energy_boson(mc::DQMC, m::HBCModel, hsfield = conf(mc)) - dtau = mc.parameters.delta_tau - lambda = acosh(exp(m.U * dtau/2)) - return lambda * sum(hsfield) -end - - -# to save only what is necessary and not rely on the model type remaining the same function MonteCarlo.save_model( file::MonteCarlo.JLDFile, m::HBCModel, entryname::String="Model" @@ -401,15 +177,15 @@ function MonteCarlo.save_model( write(file, entryname * "/t2", m.t2) write(file, entryname * "/t5", m.t5) MonteCarlo.save_lattice(file, m.l, entryname * "/l") - write(file, entryname * "/flv", m.flv) nothing end -# I used HBModel as a type name and tag before. Adding this will alow the old -# tag to still work -MonteCarlo._load(data, ::Val{:HBModel}) = MonteCarlo._load(data, Val(:HBCModel)) -function MonteCarlo._load(data, ::Val{:HBCModel}) +# load_parameters(data, ::Type{<: DQMCParameters}) +# +# Loads a DQMCParameters object from a given `data` dictionary produced by +# `JLD.load(filename)`. +function MonteCarlo._load_model(data, ::Val{:HBCModel}) if !(data["VERSION"] == 1) throw(ErrorException("Failed to load HBCModel version $(data["VERSION"])")) end @@ -422,24 +198,21 @@ function MonteCarlo._load(data, ::Val{:HBCModel}) t2 = data["t2"], t5 = data["t5"], l = l, - flv = data["flv"] ) end - ################################################################################ ### Measurement kernels ################################################################################ - -MonteCarlo.checkflavors(model::HBCModel) = MonteCarlo.checkflavors(model, 2) - -function MonteCarlo.intE_kernel(mc, model::HBCModel, G::GreensMatrix) +function MonteCarlo.intE_kernel(mc, model::HBCModel, G::GreensMatrix, ::Val{2}) # ⟨U (n↑ - 1/2)(n↓ - 1/2)⟩ = ... # = U [G↑↑ G↓↓ - G↓↑ G↑↓ - 0.5 G↑↑ - 0.5 G↓↓ + G↑↓ + 0.25] # = U [(G↑↑ - 1/2)(G↓↓ - 1/2) + G↑↓(1 + G↑↓)] # with up-up = down-down and up-down = 0 - model.U * sum((diag(G.val) .- 0.5).^2) -end \ No newline at end of file +end + +true \ No newline at end of file diff --git a/docs/src/examples/HofmannBergChowdhury.md b/docs/src/examples/HofmannBergChowdhury.md index 9790a277..9462a62a 100644 --- a/docs/src/examples/HofmannBergChowdhury.md +++ b/docs/src/examples/HofmannBergChowdhury.md @@ -167,7 +167,7 @@ Now that we have the lattice we can generate a fitting hopping matrix. But befor [LoopVectorization.jl](https://github.com/JuliaSIMD/LoopVectorization.jl) is a great tool when pushing for peak single threaded/single core linear algebra performance. The linear algebra needed for DQMC is reimplemented in MonteCarlo.jl using it for both `Float64` and `ComplexF64`. The latter uses `MonteCarlo.CMat64` and `MonteCarlo.CVec64` as concrete array types which are based on [StructArrays.jl](https://github.com/JuliaArrays/StructArrays.jl) under the hood. They should be used in this model. Furthermore we can make use of `MonteCarlo.BlockDiagonal` as we have no terms with differing spin indices. Thus we set ```julia -MonteCarlo.@with_kw_noshow struct HBCModel{LT<:AbstractLattice} <: HubbardModel +MonteCarlo.@with_kw_noshow struct HBCModel{LT<:AbstractLattice} <: Model # parameters with defaults based on paper mu::Float64 = 0.0 U::Float64 = 1.0 @@ -178,14 +178,6 @@ MonteCarlo.@with_kw_noshow struct HBCModel{LT<:AbstractLattice} <: HubbardModel # lattice l::LT - - # two fermion flavors (up, down) - flv::Int = 2 - - # temp storage to avoid allocations in propose_local and accept_local - IG::CMat64 = StructArray(Matrix{ComplexF64}(undef, length(l), 2)) - IGR::CMat64 = StructArray(Matrix{ComplexF64}(undef, length(l), 2)) - R::Diagonal{ComplexF64, CVec64} = Diagonal(StructArray(Vector{ComplexF64}(undef, 2))) end MonteCarlo.hoppingeltype(::Type{DQMC}, ::HBCModel) = ComplexF64 @@ -245,243 +237,18 @@ end We note that the hermitian conjugates of a hopping $$c_j^\dagger c_i$$ can also be understood as reversing the bond direction. Since we include both directions in our lattice definitions, second and fifth nearest neighbor hermitian conjugates are taken care of. First nearest neighbors get a phase shift from complex conjugation, which is included by swapping `t1p` and `t1m` between group one and two. -The interaction matrix can almost be copied from the repulsive Hubbard model. The only difference is that the spin up and spin down blocks get the same sign. - -```julia -@inline @bm function MonteCarlo.interaction_matrix_exp!(mc::DQMC, m::HBCModel, - result::Diagonal, conf::HubbardConf, slice::Int, power::Float64=1.) - dtau = mc.parameters.delta_tau - lambda = acosh(exp(0.5 * m.U * dtau)) - N = length(lattice(m)) - - # spin up block - @inbounds for i in 1:N - result.diag[i] = exp(sign(power) * lambda * conf[i, slice]) - end - - # spin down block - @inbounds for i in 1:N - result.diag[N+i] = exp(sign(power) * lambda * conf[i, slice]) - end - nothing -end -``` - -In this case we do not need to set the type for the interaction matrix explicitly like we did for the hopping and greens matrices, because the (abstract) `HubbardModel` already uses `Diagonal` interaction matrices. - - -## Local Updates - - -Our next task is to implement `propose_local!` and `accept_local!`. Since those only rely on specific indices, columns or rows for a large part of their calculation we have to dig into the optimized matrix types a bit. `propose_local` aims to calculate the determinant ratio $$R$$ and bosonic energy difference $$\Delta E_{Boson} = V(C_{new}) - V(c_{old})$$ where $$C$$ is the auxiliary field configuration. The determinant ratio is defined as - -```math -R = \prod_\sigma \left[ - 1 + \left( \exp(\Delta E_{Boson}) - 1 \right) - \left( 1 - G_{ii}^{\sigma, \sigma}(\tau, \tau) \right) - \right] -``` - -where $$i$$ and $$\tau$$ are the lattice index and time slice index of the proposed change in the auxiliary field. This formula already assumes that the greens matrix $$G$$ is zero for all differing spin indices (i.e. spin up-down or down-up). Therefore it is just a product of two terms. With this `propose_local` is implemented as +To finish off the mandatory model interface we need to provide two more methods. The first is `lattice(model)` which simply return the lattice of the model. The other is `nflavors(model)` which returns the number of "active" flavors of the hopping matrix. This model has two flavors total, spin up and spin down, and both of these have an active effect on the hopping matrix, i.e. the values for spin up and spin down are different. Thus this method should return 2. ```julia -@inline @bm function MonteCarlo.propose_local( - mc::DQMC, model::HBModel, i::Int, slice::Int, conf::HubbardConf - ) - N = length(model.l) - G = mc.stack.greens - Δτ = mc.parameters.delta_tau - R = model.R - - α = acosh(exp(0.5Δτ * model.U)) - ΔE_boson = -2.0α * conf[i, slice] - Δ = exp(ΔE_boson) - 1.0 - - # Unrolled R = I + Δ * (I - G) - # up-up term - R.diag.re[1] = 1.0 + Δ * (1.0 - G.blocks[1].re[i, i]) - R.diag.im[1] = - Δ * G.blocks[1].im[i, i] - # down-down term - R.diag.re[2] = 1.0 + Δ * (1.0 - G.blocks[2].re[i, i]) - R.diag.im[2] = - Δ * G.blocks[2].im[i, i] - - # Calculate "determinant" - detratio = ComplexF64( - R.diag.re[1] * R.diag.re[2] - R.diag.im[1] * R.diag.im[2], - R.diag.re[1] * R.diag.im[2] + R.diag.im[1] * R.diag.re[2] - ) - - return detratio, ΔE_boson, Δ -end -``` - -Note that the fields of our special matrix types are directly indexed here. A `BlockDiagonal` matrix contains all of its data in `B.blocks`. We define the first (upper left) block as spin up and the second (lower right) as spin down. Sitting at each block is a complex matrix represented by `CMat64`. It contains two real valued matrices at `x.re` and `x.im` representing the real and imaginary parts respectively. - -For `accept_local` we need to update the auxiliary field and the currently active greens function. To avoid recalculating $$\Delta$$ it is returned in `propose_local` and will be passed to `accept_local`. The updated greens function is given by - -```math -G_{jk}^{\sigma \sigma^\prime} = - G_{jk}^{\sigma \sigma^\prime} - - \left(I - G^{\sigma \sigma^\prime}(\tau, \tau) \right)_{ji} - R_{\sigma, \sigma^\prime}^{-1} - \Delta_{ii}^{\sigma \sigma^\prime}(\tau) - G_{ik}^{\sigma \sigma^\prime}(\tau, \tau) -``` - -where $$i$$ is again the site index of the proposed flip. Let's go through some observations/simplifications. First we note that for $$\sigma \ne \sigma^\prime$$ the greens function is and remains zero. The inversion of $$R$$ is an inversion of a diagonal matrix and thus simplifies to calculating the inverse of each element. Finally, $$\Delta$$ has the same value for spin up and spin down so it is simply a number. - -Using these observations and applying optimizations relevant to our matrix types `accept_local` can be implemented as - -```julia -@inline @bm function MonteCarlo.accept_local!( - mc::DQMC, model::HBModel, i::Int, slice::Int, conf::HubbardConf, - detratio, ΔE_boson, Δ) - - @bm "accept_local (init)" begin - N = length(model.l) - G = mc.stack.greens - IG = model.IG - IGR = model.IGR - R = model.R - end - - # compute R⁻¹ Δ, using that R is Diagonal, Δ is Number - # using Δ / (a + ib) = Δ / (a^2 + b^2) * (a - ib) - @bm "accept_local (inversion)" begin - f = Δ / (R.diag.re[1]^2 + R.diag.im[1]^2) - R.diag.re[1] = +f * R.diag.re[1] - R.diag.im[1] = -f * R.diag.im[1] - f = Δ / (R.diag.re[2]^2 + R.diag.im[2]^2) - R.diag.re[2] = +f * R.diag.re[2] - R.diag.im[2] = -f * R.diag.im[2] - end - - # Compute (I - G) R^-1 Δ - # Note IG is reduced to non-zero entries. Full IG would be - # (I-G)[:, i] 0 - # 0 (I-G)[:, i+N] - # our IG is [(I-G)[:, i] (I-G)[:, i+N]] - @bm "accept_local (IG, R)" begin - # Calculate IG = I - G (relevant entries only) - @turbo for m in axes(IG, 1) - IG.re[m, 1] = -G.blocks[1].re[m, i] - end - @turbo for m in axes(IG, 1) - IG.re[m, 2] = -G.blocks[2].re[m, i] - end - @turbo for m in axes(IG, 1) - IG.im[m, 1] = -G.blocks[1].im[m, i] - end - @turbo for m in axes(IG, 1) - IG.im[m, 2] = -G.blocks[2].im[m, i] - end - IG.re[i, 1] += 1.0 - IG.re[i, 2] += 1.0 - - # Calculate IGR = IG * R where R = R⁻¹ Δ from the - # previous calculation (relevant entries only) - # spin up-up block - @turbo for m in axes(IG, 1) - IGR.re[m, 1] = IG.re[m, 1] * R.diag.re[1] - end - @turbo for m in axes(IG, 1) - IGR.re[m, 1] -= IG.im[m, 1] * R.diag.im[1] - end - @turbo for m in axes(IG, 1) - IGR.im[m, 1] = IG.re[m, 1] * R.diag.im[1] - end - @turbo for m in axes(IG, 1) - IGR.im[m, 1] += IG.im[m, 1] * R.diag.re[1] - end - - # spin down-down block - @turbo for m in axes(IG, 1) - IGR.re[m, 2] = IG.re[m, 2] * R.diag.re[2] - end - @turbo for m in axes(IG, 1) - IGR.re[m, 2] -= IG.im[m, 2] * R.diag.im[2] - end - @turbo for m in axes(IG, 1) - IGR.im[m, 2] = IG.re[m, 2] * R.diag.im[2] - end - @turbo for m in axes(IG, 1) - IGR.im[m, 2] += IG.im[m, 2] * R.diag.re[2] - end - end - - # Update G according to G = G - (I - G)[:, i:N:end] * R⁻¹ * Δ * G[i:N:end, :] - # We already have IG = (I - G)[:, i:N:end] * R⁻¹ * Δ - @bm "accept_local (finalize computation)" begin - # get blocks to write less - G1 = G.blocks[1] - G2 = G.blocks[2] - temp1 = mc.stack.greens_temp.blocks[1] - temp2 = mc.stack.greens_temp.blocks[2] - - # compute temp = IG[:, i:N:end] * G[i:N:end, :] - # spin up-up block - @turbo for m in axes(G1, 1), n in axes(G1, 2) - temp1.re[m, n] = IGR.re[m, 1] * G1.re[i, n] - end - @turbo for m in axes(G1, 1), n in axes(G1, 2) - temp1.re[m, n] -= IGR.im[m, 1] * G1.im[i, n] - end - @turbo for m in axes(G1, 1), n in axes(G1, 2) - temp1.im[m, n] = IGR.im[m, 1] * G1.re[i, n] - end - @turbo for m in axes(G1, 1), n in axes(G1, 2) - temp1.im[m, n] += IGR.re[m, 1] * G1.im[i, n] - end - - # spin down-down block - @turbo for m in axes(G2, 1), n in axes(G2, 2) - temp2.re[m, n] = IGR.re[m, 2] * G2.re[i, n] - end - @turbo for m in axes(G2, 1), n in axes(G2, 2) - temp2.re[m, n] -= IGR.im[m, 2] * G2.im[i, n] - end - @turbo for m in axes(G2, 1), n in axes(G2, 2) - temp2.im[m, n] = IGR.im[m, 2] * G2.re[i, n] - end - @turbo for m in axes(G2, 1), n in axes(G2, 2) - temp2.im[m, n] += IGR.re[m, 2] * G2.im[i, n] - end - - # Calculate G = G - temp - # spin up-up block - @turbo for m in axes(G1, 1), n in axes(G1, 2) - G1.re[m, n] = G1.re[m, n] - temp1.re[m, n] - end - @turbo for m in axes(G1, 1), n in axes(G1, 2) - G1.im[m, n] = G1.im[m, n] - temp1.im[m, n] - end - - # spin down-down block - @turbo for m in axes(G2, 1), n in axes(G2, 2) - G2.re[m, n] = G2.re[m, n] - temp2.re[m, n] - end - @turbo for m in axes(G2, 1), n in axes(G2, 2) - G2.im[m, n] = G2.im[m, n] - temp2.im[m, n] - end - - # Update configuration - conf[i, slice] *= -1 - end - - nothing -end +MonteCarlo.lattice(m::HBCModel) = m.l +MonteCarlo.nflavors(::HBCModel) = 2 ``` - -## Utilities and other functionality - - -Now that we have the lattice, the hopping and interaction matrix as well as `propose_local` and `accept_local!` we're done implementing the model. There are a couple of things one might want to add. For example, adding `energy_boson()` would enable global updates and boson energy measurements. Adding `save_model` and `_load` should help with reducing file size and help future proof things, but isn't strictly necessary. And adding `intE_kernel` would allow the interactive and total energy to be measured. Beyond that one might add some constructors and convenience function like `parameters`. +There are a few more methods we can implement for convenience. The most important of these is `choose_field(model)`, which sets a default field for our model. The best choice here should be `DensityHirschField` or `DensityGHQField` as the model uses an attractive interaction. Beyond this we could implement `intE_kernel` to enable energy measurements, `parameters(model)`, `save_model`, `_load_model` and printing. The full code including these convenience functions can be found [here](HBC_model.jl) - # Simulation Setup @@ -492,7 +259,7 @@ To keep the runtime of this crosscheck reasonable we used the smallest linear sy uc = LatticePhysics.getUnitcellSquare(17) lpl = getLatticePeriodic(uc, 8) l = LatPhysLattice(lpl) -m = HBModel(l, t5 = 0.0, mu = -2.206) # other defaults match F = 0.2 setup +m = HBCModel(l, t5 = 0.0, mu = -2.206) # other defaults match F = 0.2 setup mc = DQMC( m, beta = beta, thermalization = 1000, sweeps = 5000, measure_rate = 5, print_rate = 100, recorder = Discarder() diff --git a/docs/src/examples/triangular_Hubbard.jl b/docs/src/examples/triangular_Hubbard.jl index 1fa4cc63..64cd192d 100644 --- a/docs/src/examples/triangular_Hubbard.jl +++ b/docs/src/examples/triangular_Hubbard.jl @@ -31,7 +31,7 @@ begin @time for beta in betas, mu in mus counter += 1 print("\r[", lpad("$counter", 2), "/$N]") - m = HubbardModelAttractive(l = lattice, t = 1.0, U = 4.0, mu = mu) + m = HubbardModel(l = lattice, t = 1.0, U = 4.0, mu = mu) dqmc = DQMC( m, beta = beta, delta_tau = 0.125, safe_mult = 8, thermalization = 1000, sweeps = 1000, measure_rate = 1, diff --git a/docs/src/examples/triangular_Hubbard.md b/docs/src/examples/triangular_Hubbard.md index d85d70f5..fc01a4ba 100644 --- a/docs/src/examples/triangular_Hubbard.md +++ b/docs/src/examples/triangular_Hubbard.md @@ -21,11 +21,11 @@ N = length(mus) * length(betas) @time for beta in betas, mu in mus counter += 1 print("\r[", lpad("$counter", 2), "/$N]") - m = HubbardModelAttractive(l = lattice, t = 1.0, U = 4.0, mu = mu) + m = HubbardModel(l = lattice, t = 1.0, U = 4.0, mu = mu) dqmc = DQMC( m, beta = beta, delta_tau = 0.125, safe_mult = 8, thermalization = 1000, sweeps = 1000, measure_rate = 1, - recorder = Discarder + recorder = Discarder() ) dqmc[:occ] = occupation(dqmc, m) dqmc[:PC] = pairing_correlation(dqmc, m, kernel = MonteCarlo.pc_kernel) diff --git a/src/flavors/DQMC/DQMC_interface.jl b/src/flavors/DQMC/DQMC_interface.jl index 75f99129..f45188e8 100644 --- a/src/flavors/DQMC/DQMC_interface.jl +++ b/src/flavors/DQMC/DQMC_interface.jl @@ -19,7 +19,7 @@ A matrix element is the hopping amplitude for a hopping process: \$j,\\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))) +hopping_matrix(m::Model) = throw(MethodError(hopping_matrix, (m, ))) nflavors(m::Model) = throw(MethodError(nflavors, (m,))) @@ -128,10 +128,27 @@ end """ greens_matrix_type(field, model) -Returns the (matrix) type of the greens and most work matrices. Defaults to -`Matrix{greens_eltype(T, m)}`. +Returns the (matrix) type of the greens and most work matrices. By default +attempts to derive the type from the hopping and interaction matrix types. """ -greens_matrix_type(f::AbstractField, m::Model) = Matrix{greens_eltype(f, m)} +function greens_matrix_type(f::AbstractField, m::Model) + IT = interaction_matrix_type(f, m) + HT = hopping_matrix_type(f, m) + T = greens_eltype(f, m) + + if IT <: BlockDiagonal && HT <: BlockDiagonal + if IT.parameters[2] == HT.parameters[2] # same number of blocks + return BlockDiagonal{T, HT.parameters[2], matrix_type(T)} + else + return matrix_type(T) + end + elseif IT <: Diagonal && HT <: BlockDiagonal + return BlockDiagonal{T, HT.parameters[2], matrix_type(T)} + else + return matrix_type(T) + end +end + ### Model @@ -140,17 +157,20 @@ greens_matrix_type(f::AbstractField, m::Model) = Matrix{greens_eltype(f, m)} """ hopping_eltype(model) -Returns the type of the elements of the hopping matrix. Defaults to `Float64`. +Returns the type of the elements of the hopping matrix. Defaults to calling +`eltype(hopping_matrix(model))`. """ -hopping_eltype(::Model) = Float64 +hopping_eltype(m::Model) = eltype(hopping_matrix(m)) """ hopping_matrix_type(field, model) -Returns the (matrix) type of the hopping matrix. Defaults to -`Matrix{hopping_eltype(model)}`. +Returns the (matrix) type of the hopping matrix. Defaults to building a full +hopping matrix and getting its type. """ -hopping_matrix_type(::AbstractField, m::Model) = Matrix{hopping_eltype(m)} +function hopping_matrix_type(f::AbstractField, m::Model) + typeof(pad_to_nflavors(f, m, hopping_matrix(m))) +end ### Field @@ -183,7 +203,7 @@ 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)) + flv = nflavors(f, m) N = length(lattice(m)) FullT = interaction_matrix_type(f, m) @@ -214,4 +234,29 @@ energy_boson(f::AbstractField, c = nothing) = throw(MethodError(energy_boson, (f conf(f::AbstractField) = f.conf conf!(f::AbstractField, c) = conf(f) .= c -temp_conf(f::AbstractField) = f.temp_conf \ No newline at end of file +temp_conf(f::AbstractField) = f.temp_conf + + +################################################################################ +### Convenience +################################################################################ + + +# These methods are just to simplify things. NOne of these should be implemented +# when extending DQMC. + +nflavors(mc::DQMC) = nflavors(field(mc), model(mc)) +nflavors(f::AbstractField, m::Model) = max(nflavors(f), nflavors(m)) + +@inline pad_to_nflavors(mc::DQMC, mat) = pad_to_nflavors(field(mc), model(mc), mat) +function pad_to_nflavors(f::AbstractField, m::Model, mat) + N = length(lattice(m)) + flv = nflavors(f, m) + if size(mat, 1) == N * flv + return mat + elseif size(mat, 1) == N + return BlockDiagonal(ntuple(_ -> mat, flv)) + else + error("Failed to expand size $(size(mat)) matrix to size ($N * $flv, $N * $flv) ") + end +end \ No newline at end of file diff --git a/src/flavors/DQMC/fields.jl b/src/flavors/DQMC/fields.jl index aec10c6b..e4431c4d 100644 --- a/src/flavors/DQMC/fields.jl +++ b/src/flavors/DQMC/fields.jl @@ -9,7 +9,7 @@ end function FieldCache(field, model) N = size(field.conf, 1) - flv = max(nflavors(field), nflavors(model)) + flv = nflavors(field, model) T = interaction_eltype(field) GET = greens_eltype(field, model) VT = vector_type(GET) diff --git a/src/flavors/DQMC/main.jl b/src/flavors/DQMC/main.jl index 6a2b0214..503b72ff 100644 --- a/src/flavors/DQMC/main.jl +++ b/src/flavors/DQMC/main.jl @@ -92,6 +92,9 @@ function DQMC( ) end +# Contains mandatory and optional method definitions expected to exist for a +# given model or field. Also sets up some convenience methods. +include("DQMC_interface.jl") # Contains different auxiliary field types include("fields.jl") @@ -123,10 +126,6 @@ include("updates/parallel_updates.jl") # and local updates. include("DQMC.jl") -# Contains mandatory and optional method definitions expected to exist for a -# given model. -include("DQMC_interface.jl") - # saving and loading of DQMC include("FileIO.jl") diff --git a/src/flavors/DQMC/measurements/extensions.jl b/src/flavors/DQMC/measurements/extensions.jl index 119ac226..94a702ee 100644 --- a/src/flavors/DQMC/measurements/extensions.jl +++ b/src/flavors/DQMC/measurements/extensions.jl @@ -158,14 +158,17 @@ end function dia_K_x(mc, G, idxs) - T = Matrix(MonteCarlo.hopping_matrix(mc, mc.model)) + if !isdefined(mc.stack, :hopping_matrix) + MonteCarlo.init_hopping_matrices(mc, mc.model) + end + T = Matrix(mc.stack.hopping_matrix) push!(mc.lattice_iterator_cache, MonteCarlo.Dir2SrcTrg(), lattice(mc)) dir2srctrg = mc.lattice_iterator_cache[MonteCarlo.Dir2SrcTrg()] N = length(lattice(mc)) Kx = ComplexF64(0) - flv = max(nflavors(field(mc)), nflavors(model(mc))) + flv = nflavors(mc) if flv == 1; f = 2.0 elseif flv == 2; f = 1.0 else error("The diamagnetic contribution to the superfluid density has no implementation for $flv flavors") diff --git a/src/flavors/DQMC/measurements/measurements.jl b/src/flavors/DQMC/measurements/measurements.jl index ac096428..9feac5e3 100644 --- a/src/flavors/DQMC/measurements/measurements.jl +++ b/src/flavors/DQMC/measurements/measurements.jl @@ -26,7 +26,7 @@ is derived from the hopping matrix and does not include on-site "hoppings". """ function hopping_directions(dqmc::DQMC, model) dir2srctrg = dqmc[Dir2SrcTrg()] - T = hopping_matrix(dqmc, model) + T = hopping_matrix(model) valid_directions = Int64[] for i in 2:length(dir2srctrg) @@ -175,7 +175,7 @@ function current_current_susceptibility( greens_iterator = CombinedGreensIterator(dqmc), wrapper = nothing, lattice_iterator = EachLocalQuadByDistance(directions), kwargs... ) - @assert is_approximately_hermitian(hopping_matrix(dqmc, model)) "CCS assumes Hermitian matrix" + @assert is_approximately_hermitian(hopping_matrix(model)) "CCS assumes Hermitian matrix" li = wrapper === nothing ? lattice_iterator : wrapper(lattice_iterator) Measurement(dqmc, model, greens_iterator, li, cc_kernel; kwargs...) end diff --git a/src/flavors/DQMC/stack.jl b/src/flavors/DQMC/stack.jl index 992f2670..dbf94a18 100644 --- a/src/flavors/DQMC/stack.jl +++ b/src/flavors/DQMC/stack.jl @@ -201,14 +201,8 @@ function init_hopping_matrices(mc::DQMC{M,CB}, m::Model) where {M, CB<:Checkerbo nothing end function init_hopping_matrix_exp(mc::DQMC, m::Model) - N = length(lattice(m)) - flv = nflavors(field(mc)) dtau = mc.parameters.delta_tau - - T = hopping_matrix(mc, m) - size(T) == (flv*N, flv*N) || error( - "Hopping matrix should have size $((flv*N, flv*N)) but has size $(size(T))." - ) + T = pad_to_nflavors(mc, hopping_matrix(m)) if !is_approximately_hermitian(T) @error( @@ -242,7 +236,7 @@ rem_eff_zeros!(X::AbstractArray) = map!(e -> abs.(e)<1e-15 ? zero(e) : e,X,X) function init_checkerboard_matrices(mc::DQMC, m::Model) s = mc.stack l = lattice(m) - flv = nflavors(field(mc)) + flv = nflavors(mc) H = heltype(mc) N = length(l) dtau = mc.parameters.delta_tau @@ -252,7 +246,7 @@ function init_checkerboard_matrices(mc::DQMC, m::Model) n_grps = s.n_groups cb = s.checkerboard - T = reshape(hopping_matrix(mc, m), (N, flv, N, flv)) + T = reshape(pad_to_nflavors(mc, hopping_matrix(m)), (N, flv, N, flv)) s.chkr_hop_half = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps) s.chkr_hop_half_inv = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps) diff --git a/src/models/DummyModel.jl b/src/models/DummyModel.jl index 9fa177b3..56f7a719 100644 --- a/src/models/DummyModel.jl +++ b/src/models/DummyModel.jl @@ -21,7 +21,8 @@ end choose_field(::DummyModel) = DensityHirschField nflavors(::DummyModel) = 1 -lattice(m::DummyModel) = m.l +lattice(m::DummyModel) = get(m.data, "l", Chain(1)) +hopping_matrix(m::DummyModel) = fill(1.0, length(lattice(m)), length(lattice(m))) function save_model(file::JLDFile, m::DummyModel, entryname::String="Model") diff --git a/src/models/HubbardModel.jl b/src/models/HubbardModel.jl index eb569106..c74a12bd 100644 --- a/src/models/HubbardModel.jl +++ b/src/models/HubbardModel.jl @@ -86,7 +86,7 @@ nflavors(::HubbardModel) = 1 hopping_eltype(model::HubbardModel) = typeof(model.t) function hopping_matrix_type(field::AbstractField, model::HubbardModel) - flv = max(nflavors(field), nflavors(model)) + flv = nflavors(field, model) T = hopping_eltype(model) MT = matrix_type(T) if flv == 1 @@ -96,24 +96,9 @@ function hopping_matrix_type(field::AbstractField, model::HubbardModel) end end -function greens_matrix_type(f::AbstractHirschField{T}, m::HubbardModel) where {T} - if max(nflavors(f), nflavors(m)) == 1 - return matrix_type(T) - else - return BlockDiagonal{T, 2, matrix_type(T)} - end -end -function greens_matrix_type(f::AbstractGHQField{T}, m::HubbardModel) where {T} - if max(nflavors(f), nflavors(m)) == 1 - return matrix_type(T) - else - return BlockDiagonal{T, 2, matrix_type(T)} - end -end - """ - hopping_matrix(mc, model) + hopping_matrix(model) Calculates the hopping matrix \$T_{i, j}\$ where \$i, j\$ are site indices. @@ -121,7 +106,7 @@ site indices. This isn't a performance critical method as it is only used once before the actual simulation. """ -function hopping_matrix(mc::DQMC, m::HubbardModel) +function hopping_matrix(m::HubbardModel) N = length(m.l) T = diagm(0 => fill(-m.mu, N)) @@ -133,11 +118,7 @@ function hopping_matrix(mc::DQMC, m::HubbardModel) end end - if max(nflavors(field(mc)), nflavors(m)) == 1 - return T - else - return BlockDiagonal(T, copy(T)) - end + return T end function save_model(file::JLDFile, m::HubbardModel, entryname::String = "Model") diff --git a/test/ED/ED_tests.jl b/test/ED/ED_tests.jl index a735a618..e6e2c234 100644 --- a/test/ED/ED_tests.jl +++ b/test/ED/ED_tests.jl @@ -159,7 +159,7 @@ end N = length(lattice(model)) # Direct calculation simialr to what DQMC should be doing - T = Matrix(MonteCarlo.hopping_matrix(dqmc, model)) + T = Matrix(dqmc.stack.hopping_matrix) # Doing an eigenvalue decomposition makes this pretty stable vals, U = eigen(exp(-T)) D = Diagonal(vals)^(dqmc.parameters.beta) diff --git a/test/fields.jl b/test/fields.jl index a210f6ed..2d06c074 100644 --- a/test/fields.jl +++ b/test/fields.jl @@ -1,50 +1,70 @@ # these are just shorthads -using MonteCarlo: FVec64, FMat64, CVec64, CMat64 +using MonteCarlo: FVec64, FMat64, CVec64, CMat64, BlockDiagonal -@testset "Field Cache & Interaction Matrix" begin +@testset "Typing (Field Cache, hopping, interaction & greens matrix" begin T2 = NTuple{2} - - mc = DQMC(HubbardModel(8, 1, U = 1.0), field = DensityHirschField, beta = 1.0) - MonteCarlo.init!(mc) - @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{Float64, Float64, FVec64, Float64} - @test MonteCarlo.interaction_eltype(mc.field) == Float64 - @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{Float64, FVec64} - @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV) - - mc = DQMC(HubbardModel(8, 1, U = -1.0), field = DensityHirschField, beta = 1.0) - MonteCarlo.init!(mc) - @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{ComplexF64, ComplexF64, CVec64, ComplexF64} - @test MonteCarlo.interaction_eltype(mc.field) == ComplexF64 - @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{ComplexF64, CVec64} - @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV) - - mc = DQMC(HubbardModel(8, 1, U = 1.0), field = MagneticHirschField, beta = 1.0) - MonteCarlo.init!(mc) - @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{CVec64, CVec64, T2{CVec64}, ComplexF64} - @test MonteCarlo.interaction_eltype(mc.field) == ComplexF64 - @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{ComplexF64, CVec64} - @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV) - - mc = DQMC(HubbardModel(8, 1, U = -1.0), field = MagneticHirschField, beta = 1.0) - MonteCarlo.init!(mc) - @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{FVec64, FVec64, T2{FVec64}, Float64} - @test MonteCarlo.interaction_eltype(mc.field) == Float64 - @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{Float64, FVec64} - @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV) - - mc = DQMC(HubbardModel(8, 1, U = 1.0), field = MagneticGHQField, beta = 1.0) - MonteCarlo.init!(mc) - @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{CVec64, CVec64, T2{CVec64}, ComplexF64} - @test MonteCarlo.interaction_eltype(mc.field) == ComplexF64 - @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{ComplexF64, CVec64} - @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV) - - mc = DQMC(HubbardModel(8, 1, U = -1.0), field = MagneticGHQField, beta = 1.0) - MonteCarlo.init!(mc) - @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{FVec64, FVec64, T2{FVec64}, Float64} - @test MonteCarlo.interaction_eltype(mc.field) == Float64 - @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{Float64, FVec64} - @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV) + fields = (DensityHirschField, MagneticHirschField, DensityGHQField, MagneticGHQField) + + for U in (+1.0, -1.0) + for (i, field) in enumerate(fields) + + @testset "U = $U field = $field" begin + mc = DQMC(HubbardModel(U = U), field = field, beta = 1.0) + MonteCarlo.init!(mc) + + # missmatch between U and Density/Magnetic results in complex interaction + if (i + (U == +1.0)) % 2 == 0 # (U == 1 && Density) || (U == -1 && Magnetic) + T = Float64; VT = FVec64; mT = FMat64 + else + T = ComplexF64; VT = CVec64; mT = CMat64 + end + @test MonteCarlo.interaction_eltype(mc.field) == T + @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{T, VT} + @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa Diagonal{T, VT} + @test mc.stack.eV isa Diagonal{T, VT} + + # hopping are always real and 1 flavor. They should get dublicated + # if interaction requires two flavors (magnetic) + @test MonteCarlo.nflavors(mc.model) == 1 + @test MonteCarlo.hopping_eltype(mc.model) == Float64 + + # Greens matrix takes eltype from interaction because that might be complex + @test MonteCarlo.greens_eltype(mc.field, mc.model) == T + + # 1 flavor for Density, 2 for Magnetic in interaction/total + @test MonteCarlo.nflavors(mc.field) == 2 - (i % 2) + @test MonteCarlo.nflavors(mc.field, mc.model) == 2 - (i % 2) + + x = rand(4, 4) + + if i % 2 == 1 + @test MonteCarlo.hopping_matrix_type(mc.field, mc.model) == FMat64 + @test mc.stack.hopping_matrix isa FMat64 + + @test MonteCarlo.pad_to_nflavors(mc.field, mc.model, x) == x + + @test MonteCarlo.greens_matrix_type(mc.field, mc.model) == mT + @test mc.stack.greens isa mT + + # single flavor -> single value of greens eltype + @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{T, T, VT, T} + else + @test MonteCarlo.hopping_matrix_type(mc.field, mc.model) == BlockDiagonal{Float64, 2, FMat64} + @test mc.stack.hopping_matrix isa BlockDiagonal{Float64, 2, FMat64} + + @test MonteCarlo.pad_to_nflavors(mc.field, mc.model, x) == BlockDiagonal(x, x) + + @test MonteCarlo.greens_matrix_type(mc.field, mc.model) == BlockDiagonal{T, 2, mT} + @test mc.stack.greens isa BlockDiagonal{T, 2, mT} + + # two flavor -> matrix types of greens eltypes + @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{VT, VT, T2{VT}, T} + end + end + + end + end + end @testset "Lookup tables" begin diff --git a/test/flavortests_DQMC.jl b/test/flavortests_DQMC.jl index a903715d..4adb9579 100644 --- a/test/flavortests_DQMC.jl +++ b/test/flavortests_DQMC.jl @@ -70,7 +70,7 @@ end # DQMC mandatory m = DummyModel(SquareLattice(2)) f = DummyField() - @test_throws MethodError MonteCarlo.hopping_matrix(dqmc, m) + @test_throws MethodError MonteCarlo.hopping_matrix(m) @test_throws MethodError rand(f) @test_throws MethodError rand!(f) @test_throws MethodError MonteCarlo.nflavors(m) @@ -92,13 +92,14 @@ end MonteCarlo.nflavors(::DummyField) = 2 # DQMC optional - @test MonteCarlo.hopping_eltype(m) == Float64 + # method errors come from hopping_matrix(model) + @test_throws MethodError MonteCarlo.hopping_eltype(m) @test MonteCarlo.interaction_eltype(f) == Float64 - @test MonteCarlo.greens_eltype(f, m) == Float64 + @test_throws MethodError MonteCarlo.greens_eltype(f, m) == Float64 - @test MonteCarlo.hopping_matrix_type(f, m) == Matrix{Float64} + @test_throws MethodError MonteCarlo.hopping_matrix_type(f, m) @test MonteCarlo.interaction_matrix_type(f, m) == Matrix{Float64} - @test MonteCarlo.greens_matrix_type(f, m) == Matrix{Float64} + @test_throws MethodError MonteCarlo.greens_matrix_type(f, m) == Matrix{Float64} M = MonteCarlo.init_interaction_matrix(f, m) @test M isa Matrix{Float64} @@ -565,7 +566,7 @@ end N = length(lattice(model)) # Direct calculation similar to what DQMC should be doing - T = Matrix(MonteCarlo.hopping_matrix(dqmc, model)) + T = Matrix(dqmc.stack.hopping_matrix) # Doing an eigenvalue decomposition makes this pretty stable vals, U = eigen(exp(-T)) D = Diagonal(vals)^(dqmc.parameters.beta)