Rimu

Random integrators for many-body quantum systems

Welcome to Rimu version 0.11.2. Read the documentation online.

Rimu.PACKAGE_VERSION

Constant that contains the current VersionNumber of Rimu.

AllOverlaps(num_replicas=2; operator=nothing, transform=nothing, vecnorm=true) <: ReplicaStrategy{num_replicas}

Run num_replicas replicas and report overlaps between all pairs of replica vectors. If operator is not nothing, the overlap dot(c1, operator, c2) is reported as well. If operator is a tuple of operators, the overlaps are computed for all operators.

Column names in the report are of the form c{i}_dot_c{j} for vector-vector overlaps, and c{i}_Op{k}_c{j} for operator overlaps.

See lomc!, ReplicaStrategy and AbstractHamiltonian (for an interface for implementing operators).

Transformed Hamiltonians

If a transformed Hamiltonian G has been passed to lomc! then overlaps can be calculated by passing the same transformed Hamiltonian to AllOverlaps by setting transform=G. A warning is given if these two Hamiltonians do not match.

Implemented transformations are:

• GutzwillerSampling
• GuidingVectorSampling

In the case of a transformed Hamiltonian the overlaps are defined as follows. For a similarity transformation G of the Hamiltonian (see e.g. GutzwillerSampling.)

\[ \hat{G} = f \hat{H} f^{-1}.\]

The expectation value of an operator $\hat{A}$ is

\[ \langle \hat{A} \rangle = \langle \psi | \hat{A} | \psi \rangle = \frac{\langle \phi | f^{-1} \hat{A} f^{-1} | \phi \rangle}{\langle \phi | f^{-2} | \phi \rangle}\]

where

\[ | \phi \rangle = f | \psi \rangle\]

is the (right) eigenvector of $\hat{G}$ and $| \psi \rangle$ is an eigenvector of $\hat{H}$.

For a K-tuple of input operators (\hat{A}_1, ..., \hat{A}_K), overlaps of $\langle \phi | f^{-1} \hat{A} f^{-1} | \phi \rangle$ are reported as c{i}_Op{k}_c{j}. The correct vector-vector overlap $\langle \phi | f^{-2} | \phi \rangle$ is reported last as c{i}_Op{K+1}_c{j}. This is in addition to the bare vector-vector overlap $\langle \phi | f^{-2} | \phi \rangle$ that is reported as c{i}_dot_c{j}.

In either case the c{i}_dot_c{j} overlap can be omitted with the flag vecnorm=false.

ConstantTimeStep <: TimeStepStrategy

Keep dτ constant.

DontUpdate(; targetwalkers = 1_000_000) <: ShiftStrategy

Don't update the shift. Return when targetwalkers is reached.

See ShiftStrategy, lomc!.

DoubleLogProjected(; target, projector, ζ = 0.08, ξ = ζ^2/4) <: ShiftStrategy

Strategy for updating the shift according to the log formula with damping parameter ζ and ξ after projecting onto projector.

\[S^{n+1} = S^n -\frac{ζ}{dτ}\ln\left(\frac{P⋅Ψ^{(n+1)}}{P⋅Ψ^{(n)}}\right)-\frac{ξ}{dτ}\ln\left(\frac{P⋅Ψ^{(n+1)}}{\text{target}}\right)\]

Note that adjusting the keyword maxlength in lomc! is advised as the default may not be appropriate.

See ShiftStrategy, lomc!.

DoubleLogSumUpdate(; targetwalkers = 1000, ζ = 0.08, ξ = ζ^2/4, α = 1/2) <: ShiftStrategy

Strategy for updating the shift according to the log formula with damping parameters ζ and ξ.

\[S^{n+1} = S^n -\frac{ζ}{dτ}\ln\left(\frac{N_\mathrm{w}^{n+1}}{N_\mathrm{w}^n}\right) - \frac{ξ}{dτ}\ln\left(\frac{N_\mathrm{w}^{n+1}}{N_\mathrm{w}^\text{target}}\right),\]

where $N_\mathrm{w} =$ (1-α)*walkernumber() + α*UniformProjector()⋅ψ computed with walkernumber() and UniformProjector(). When ξ = ζ^2/4 this corresponds to critical damping with a damping time scale T = 2/ζ.

See ShiftStrategy, lomc!.

DoubleLogUpdate(; targetwalkers = 1000, ζ = 0.08, ξ = ζ^2/4) <: ShiftStrategy

Strategy for updating the shift according to the log formula with damping parameter ζ and ξ.

\[S^{n+1} = S^n -\frac{ζ}{dτ}\ln\left(\frac{\|Ψ\|_1^{n+1}}{\|Ψ\|_1^n}\right)-\frac{ξ}{dτ}\ln\left(\frac{\|Ψ\|_1^{n+1}}{\|Ψ\|_1^\text{target}}\right)\]

When ξ = ζ^2/4 this corresponds to critical damping with a damping time scale T = 2/ζ.

See ShiftStrategy, lomc!.

DoubleLogUpdateAfterTargetWalkers(targetwalkers, ζ = 0.08, ξ = 0.0016) <: ShiftStrategy

Strategy for updating the shift: After targetwalkers is reached, update the shift according to the log formula with damping parameter ζ and ξ.

See DoubleLogUpdate, ShiftStrategy, lomc!.

 FciqmcRunStrategy{T}

Abstract type representing the strategy for running and terminating lomc!(). The type parameter T is relevant for reporting the shift and the norm.

Implemented strategies:

• RunTillLastStep

FirstOrderTransitionOperator(hamiltonian, shift, dτ) <: AbstractHamiltonian

First order transition operator

\[𝐓 = 1 + dτ(S - 𝐇)\]

where $𝐇$ is the hamiltonian and $S$ is the shift.

$𝐓$ represents the first order expansion of the exponential evolution operator of the imaginary-time Schrödinger equation (Euler step) and repeated application will project out the ground state eigenvector of the hamiltonian. It is the transition operator used in FCIQMC.

LogUpdate(ζ = 0.08) <: ShiftStrategy

Strategy for updating the shift according to the log formula with damping parameter ζ.

\[S^{n+1} = S^n -\frac{ζ}{dτ}\ln\left(\frac{\|Ψ\|_1^{n+1}}{\|Ψ\|_1^n}\right)\]

See ShiftStrategy, lomc!.

LogUpdateAfterTargetWalkers(targetwalkers, ζ = 0.08) <: ShiftStrategy

Strategy for updating the shift: After targetwalkers is reached, update the shift according to the log formula with damping parameter ζ.

See LogUpdate, ShiftStrategy, lomc!.

MultiScalar

Wrapper over a tuple that supports +, *, min, and max. Used with MPI communication because SVectors are treated as arrays by MPI.Allreduce and Tuples do not support scalar operations.

Example

Suppose you want to compute the sum of a vector dv and also get the number of positive elements it has in a single pass. You can use MultiScalar:

julia> dv = DVec(:a => 1, :b => -2, :c => 1);

julia> s, p = mapreduce(+, values(dv)) do v
    Rimu.MultiScalar(v, Int(sign(v) == 1))
end;

julia> s, p
(0, 2)

This will work with MPIData.

Note that only MultiScalars with the same types can be operated on. This is a feature, as it forces type stability.

NoStats(N=1) <: ReplicaStrategy{N}

The default ReplicaStrategy. N replicas are run, but no statistics are collected.

See also lomc!.

PostStepStrategy

Subtypes of PostStepStrategy can be used to perform arbitrary computation on a replica after an FCIQMC step is finished and report the results.

Implemented strategies:

• ProjectedEnergy
• Projector
• SignCoherence
• WalkerLoneliness
• Timer

Note: A tuple of multiple strategies can be passed to lomc!. In that case, all reported column names must be distinct.

Interface:

A subtype of this type must implement post_step(::PostStepStrategy, ::ReplicaState).

ProjectedEnergy(hamiltonian, projector; hproj=:hproj, vproj=:vproj) <: PostStepStrategy

After every step, compute hproj = dot(projector, hamiltonian, dv) and vproj = dot(projector, dv), where dv is the instantaneous coefficient vector. projector can be an AbstractDVec, or an AbstractProjector.

Reports to columns hproj and vproj, which can be used to compute projective energy, e.g. with projected_energy. The keyword arguments hproj and vproj can be used to change the names of these columns. This can be used to make the names unique when computing projected energies with different projectors in the same run.

See also projected_energy, ratio_of_means, mixed_estimator, and PostStepStrategy.

Projector(name=projector) <: PostStepStrategy

After each step, compute dot(projector, dv) and report it in the DataFrame under name. projector can be an AbstractDVec, or an AbstractProjector.

QMCState

Holds all information needed to run lomc!, except the dataframe. Holds an NTuple of ReplicaStates, the Hamiltonian, and various strategies that control the algorithm. Constructed and returned by lomc!.

ReplicaState(v, wm, pnorm, params, id)

Struct that holds all information needed for an independent run of the algorithm.

Can be advanced a step forward with advance!.

Fields

• hamiltonian: the model Hamiltonian.
• v: vector.
• pv: vector from the previous step.
• wm: working memory.
• pnorm: previous walker number (see walkernumber).
• params: the FciqmcRunStrategy.
• id: string ID appended to reported column names.

See also QMCState, ReplicaStrategy, replica_stats, lomc!.

ReplicaStrategy{N}

Supertype for strategies that can be passed to lomc! and control how many replicas are used, and what information is computed and returned. The number of replicas is N.

Concrete implementations

• NoStats: run (possibly one) replica(s), but don't report any additional info.
• AllOverlaps: report overlaps between all pairs of replica vectors.

Interface

A subtype of ReplicaStrategy{N} must implement the following function:

• Rimu.replica_stats - return a tuple of Strings or Symbols of names for replica statistics and a tuple of the values. These will be reported to the DataFrame returned by lomc!.

Report()

Internal structure that holds the temporary reported values as well as metadata.

See report!, report_metadata!, get_metadata.

ReportDFAndInfo(; reporting_interval=1, info_interval=100, io=stdout, writeinfo=false) <: ReportingStrategy

The default ReportingStrategy. Report every reporting_intervalth step to a DataFrame and write info message to io every info_intervalth reported step (unless writeinfo == false). The flag writeinfo is useful for controlling info messages in MPI codes, e.g. by setting writeinfo =is_mpi_root().

ReportToFile(; kwargs...) <: ReportingStrategy

ReportingStrategy that writes the report directly to a file in the Arrow format. Useful when dealing with long jobs or large numbers of replicas, when the report can incur a significant memory cost.

The arrow file can be read back in with load_df(filename) or using Arrow; Arrow.Table(filename).

Keyword arguments

• filename = "out.arrow": the file to report to. If the file already exists, a new file is created.
• reporting_interval = 1: interval between simulation steps that are reported.
• chunk_size = 1000: the size of each chunk that is written to the file. A DataFrame of this size is collected in memory and written to disk. When saving, an info message is also printed to io.
• save_if =is_mpi_root(): if this value is true, save the report, otherwise ignore it.
• return_df = false: if this value is true, read the file and return the data frame at the end of computation. Otherwise, an empty DataFrame is returned.
• io = stdout: The IO to print messages to. Set to devnull if you don't want to see messages printed out.
• compress = :zstd: compression algorithm to use. Can be :zstd, :lz4 or nothing.

See also load_df and save_df.

ReportingStrategy

Abstract type for strategies for reporting data in a DataFrame with report!().

Implemented strategies:

• ReportDFAndInfo
• ReportToFile

Interface:

A ReportingStrategy can define any of the following:

• refine_r_strat
• report!
• report_after_step
• finalize_report!
• reporting_interval

RunTillLastStep(step::Int = 0 # number of current/starting timestep
             laststep::Int = 100 # number of final timestep
             shiftMode::Bool = false # whether to adjust shift
             shift = 0.0 # starting/current value of shift
             dτ::Float64 = 0.01 # current value of time step
) <: FciqmcRunStrategy

Parameters for running lomc!() for a fixed number of time steps. For alternative strategies, see FciqmcRunStrategy.

Abstract type for defining the strategy for updating the shift. Passed as a parameter to lomc!.

Implemented strategies:

• DontUpdate
• DoubleLogUpdate - default in lomc!()
• LogUpdate
• LogUpdateAfterTargetWalkers - FCIQMC standard
• DoubleLogUpdateAfterTargetWalkers

SignCoherence(reference[; name=:coherence]) <: PostStepStrategy

After each step, compute the proportion of configurations that have the same sign as they do in the reference_dvec. Reports to a column named name, which defaults to coherence.

SingleParticleDensity(; save_every=1, component) <: PostStepStrategy

PostStepStrategy to compute the diagonal single_particle_density. It records a Tuple with the same eltype as the vector.

Computing the density at every time step can be expensive. This cost can be reduced by setting the save_every argument to a higher value. If the value is set, a vector of zeros is recorded when the saving is skipped.

If the address type has multiple components, the component argument can be used to compute the density on a per-component basis.

The density is not normalized, and must be divided by the vector norm(⋅,2) squared.

See also

• single_particle_density
• DensityMatrixDiagonal

TimeStepStrategy

Abstract type for strategies for updating the time step with update_dτ(). Implemented strategies:

• ConstantTimeStep

Timer <: PostStepStrategy

Record current time after every step. See Base.Libc.time for information on what time is recorded.

TripleLogUpdate(; targetwalkers = 1000, ζ = 0.08, ξ = ζ^2/4, η = 0.01) <: ShiftStrategy

Strategy for updating the shift according to the extended log formula with damping parameters ζ, ξ, and η.

\[S^{n+1} = S^n -\frac{ζ}{dτ}\ln\left(\frac{N_\mathrm{w}^{n+1}}{N_\mathrm{w}^n}\right) - \frac{ξ}{dτ}\ln\left(\frac{N_\mathrm{w}^{n+1}}{N_\mathrm{w}^\text{target}}\right) - \frac{η}{dτ}\ln\left(\frac{\|ℜ(Ψ^{n+1})\|_1^2 + \|ℑ(Ψ^{n+1})\|_1^2} {\|ℜ(Ψ^{n})\|_1^2 + \|ℑ(Ψ^{n})\|_1^2}\right),\]

where $N_\mathrm{w}$ is the walkernumber(). When ξ = ζ^2/4 this corresponds to critical damping with a damping time scale T = 2/ζ.

See ShiftStrategy, lomc!.

WalkerLoneliness(threshold=1) <: PostStepStrategy

After each step, compute the proportion of configurations that are occupied by at most threshold walkers. Reports to a column named loneliness.

_n_walkers(v, s_strat)

Returns an estimate of the expected number of walkers as an integer.

advance!(report::Report, state::QMCState, replica::ReplicaState)

Advance the replica by one step. The state is used only to access the various strategies involved. Steps, stats, and computed quantities are written to the report.

Returns true if the step was successful and calculation should proceed, false when it should terminate.

all_overlaps(operators, vectors, working_memories, vecnorm=true)

Get all overlaps between vectors and operators. This function is overloaded for MPIData. The flag vecnorm can disable the vector-vector overlap c{i}_dot_c{j}.

check_transform(r::AllOverlaps, ham)

Check that the transformation provided to r::AllOverlaps matches the given Hamiltonian ham. Used as a sanity check before starting main lomc! loop.

default_logger(args...)

Reset the global_logger to Logging.ConsoleLogger. Undoes the effect of smart_logger. Arguments are passed on to Logging.ConsoleLogger.

default_starting_vector(hamiltonian::AbstractHamiltonian; kwargs...)
default_starting_vector(
    address=starting_address(hamiltonian);
    style=IsStochasticInteger(),
    initiator=NonInitiator(),
    threading=nothing
)

Return a default starting vector for lomc!. The default choice for the starting vector is

v = PDVec(address => 10; style, initiator)

if threading is available, or otherwise

v = DVec(address => 10; style)

if initiator == NonInitiator(), and

v = InitiatorDVec(address => 10; style, initiator)

if not. See PDVec, DVec, InitiatorDVec, StochasticStyle, and [InitiatorRule].

finalize_report!(::ReportingStrategy, report)

Finalize the report. This function is called after all steps in lomc! have finished.

get_metadata(report::Report, key)

Get metadata key from report. key is converted to a String.

See also report_metadata!, Report, report!.

lomc!(ham::AbstractHamiltonian, [v]; kwargs...) -> df, state
lomc!(state::QMCState, [df]; kwargs...) -> df, state

Linear operator Monte Carlo: Perform a projector quantum Monte Carlo simulation for determining the lowest eigenvalue of ham. The details of the simulation are controlled by the optional keyword arguments and by the type of the optional starting vector v. Alternatively, a QMCState can be passed in to continue a previous simulation.

Common keyword arguments and defaults:

• laststep = 100 - controls the number of steps.
• dτ = 0.01 - time step.
• targetwalkers = 1000 - target for the 1-norm of the coefficient vector.
• address = starting_address(ham) - set starting address for default v and shift.
• style = IsStochasticInteger() - set StochasticStyle for default v; unused if v is specified.
• initiator = NonInitiator() - set InitiatorRule for default v; unused if v is specified.
• threading - default is to use multithreading and MPI if multiple threads are available. Set to true to force PDVec for the starting vector, false for serial computation; unused if v is specified.
• shift = diagonal_element(ham, address) - initial value of shift.
• post_step::NTuple{N,<:PostStepStrategy} = () - extract observables (e.g. ProjectedEnergy), see PostStepStrategy.
• replica::ReplicaStrategy = NoStats(1) - run several synchronised simulations, see ReplicaStrategy.
• r_strat::ReportingStrategy = ReportDFAndInfo() - how and when to report results, see ReportingStrategy
• name = "lomc!" - name displayed in progress bar (via ProgressLogging)
• metadata - user-supplied metadata to be added to the report df. Must be an iterable of pairs or a NamedTuple, e.g. metadata = ("key1" => "value1", "key2" => "value2"). All metadata is converted to strings.

Some metadata is automatically added to the report df including Rimu.PACKAGE_VERSION and data from state.

Return values

lomc! returns a named tuple with the following fields:

• df: a DataFrame with all statistics being reported.
• state: a QMCState that can be used for continuations.

Example

julia> add = BoseFS(1,2,3);

julia> hamiltonian = HubbardReal1D(add);

julia> df1, state = lomc!(hamiltonian; targetwalkers=500, laststep=100);

julia> df2, _ = lomc!(state, df1; laststep=200, metadata=(;info="cont")); # Continuation run

julia> size(df1)
(100, 11)

julia> size(df2)
(200, 11)

julia> using DataFrames; metadata(df2, "info") # retrieve custom metadata
"cont"

julia> metadata(df2, "hamiltonian") # some metadata is automatically added
"HubbardReal1D(BoseFS{6,3}(1, 2, 3); u=1.0, t=1.0)"

Further keyword arguments and defaults:

• τ_strat::TimeStepStrategy = ConstantTimeStep() - adjust time step or not, see TimeStepStrategy
• s_strat::ShiftStrategy = DoubleLogUpdate(; targetwalkers, ζ = 0.08, ξ = ζ^2/4) - how to update the shift, see ShiftStrategy.
• maxlength = 2 * s_strat.targetwalkers + 100 - upper limit on the length of v; when reached, lomc! will abort
• params::FciqmcRunStrategy = RunTillLastStep(laststep = 100, dτ = 0.01, shift = diagonal_element(ham, address) - basic parameters of simulation state, see FciqmcRunStrategy. Parameter values are overridden by explicit keyword arguments laststep, dτ, shift; is mutated.
• wm - working memory for re-use in subsequent calculations; is mutated.
• df = DataFrame() - when called with AbstractHamiltonian argument, a DataFrame can be passed for merging with the report df.

The default choice for the starting vector is v = default_starting_vector(; address, style, threading, initiator). See default_starting_vector, PDVec, DVec, StochasticStyle, and InitiatorRule.

post_step(::PostStepStrategy, ::ReplicaState) -> kvpairs

Compute statistics after FCIQMC step. Should return a tuple of :key => value pairs. This function is only called every reporting_interval steps, as defined by the ReportingStrategy.

See also PostStepStrategy, ReportingStrategy.

refine_r_strat(r_strat::ReportingStrategy) -> r_strat

Initialize the reporting strategy. This can be used to set up filenames or other attributes that need to be unique for a run of FCIQMC.

replica_stats(RS::ReplicaStrategy{N}, replicas::NTuple{N,ReplicaState}) -> (names, values)

Return the names and values of statistics related to N replicas consistent with the ReplicaStrategy RS. names should be a tuple of Symbols or Strings and values should be a tuple of the same length. This function will be called every reporting_interval steps from lomc!, or once per time step if reporting_interval is not defined.

Part of the ReplicaStrategy interface. See also ReplicaState.

 report!(::ReportingStrategy, step, report::Report, keys, values, id="")
 report!(::ReportingStrategy, step, report::Report, nt, id="")

Report keys and values to report, which will be converted to a DataFrame before lomc! exits. Alternatively, a nt::NamedTuple can be passed in place of keys and values. If id is specified, it is appended to all keys. This is used to differentiate between values reported by different replicas.

To overload this function for a new ReportingStrategy, overload report!(::ReportingStrategy, step, args...) and apply the report by calling report!(args...).

report!(report, keys, values, id="")
report!(report, pairs, id="")

Write keys, values pairs to report that will be converted to a DataFrame later. Alternatively, a named tuple or a collection of pairs can be passed instead of keys and values.

The value of id is appended to the name of the column, e.g. report!(report, :key, value, :_1) will report value to a column named :key_1.

report_after_step(::ReportingStrategy, step, report, state)

This function is called at the very end of a step, after reporting_interval steps. For example, it can be used to print some information to stdout.

report_metadata!(report::Report, key, value)
report_metadata!(report::Report, kvpairs)

Set metadata key to value in report. key and value are converted to Strings. Alternatively, an iterable of key-value pairs or a NamedTuple can be passed.

Throws an error if key already exists.

See also get_metadata, report!, Report.

reporting_interval(::ReportingStrategy)

Get the interval between steps for which non-essential statistics are reported. Defaults to 1 if chosen ReportingStrategy does not specify an interval.

single_particle_density(dvec; component)
single_particle_density(add; component)

Compute the diagonal single particle density of vector dvec or address add. If the component argument is given, only that component of the addresses is taken into account. The result is always normalized so that sum(result) ≈ num_particles(address).

Examples

julia> v = DVec(fs"|⋅↑⇅↓⋅⟩" => 1.0, fs"|↓↓⋅↑↑⟩" => 0.5)
DVec{FermiFS2C{2, 2, 5, 4, FermiFS{2, 5, BitString{5, 1, UInt8}}, FermiFS{2, 5, BitString{5, 1, UInt8}}},Float64} with 2 entries, style = IsDeterministic{Float64}()
  fs"|↓↓⋅↑↑⟩" => 0.5
  fs"|⋅↑⇅↓⋅⟩" => 1.0

julia> single_particle_density(v)
(0.2, 1.0, 1.6, 1.0, 0.2)

julia> single_particle_density(v; component=1)
(0.0, 1.6, 1.6, 0.4, 0.4)

See also

• SingleParticleDensity

smart_logger(args...)

Enable terminal progress bar during interactive use (i.e. unless running on CI or HPC). Arguments are passed on to the logger. This is run once during Rimu startup. Undo with default_logger or by setting Base.global_logger().

update_dτ(s<:TimeStepStrategy, dτ, tnorm) -> new dτ

Update the time step according to the strategy s.

update_shift(s <: ShiftStrategy, shift, shiftMode, tnorm, pnorm, dτ, step, df, v_new, v_old)

Update the shift according to strategy s. See ShiftStrategy.