14949967345

Committed 10 May 2025 11:08PM UTC coverage: 94.621% (-0.09%) from 94.706%

Build # 14949967345

Build Type

Pull #314

github

Committed by

joachimbrand

Commit Message

VMSize and get_vmsize

Pull Request Pull Request #314: Provide guidance for exact diagonalization algorithms

Run Details

76 of 84 new or added lines in 10 files covered. (90.48%)

1 existing line in 1 file now uncovered.

6983 of 7380 relevant lines covered (94.62%)

11371990.96 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.0

/src/Hamiltonians/abstract.jl

"""
    check_address_type(h::AbstractObservable, addr_or_type)
Throw an `ArgumentError` if `addr_or_type` is not compatible with `h`, otherwise return
`true`. Acceptable arguments are either an address or an address type, or a tuple or
array thereof.

See also [`allows_address_type`](@ref).
"""
@inline function check_address_type(h, ::Type{A}) where {A}
    if !allows_address_type(h, A)
        throw(ArgumentError("address type $A not allowed for operator $h"))
    end
    return true
end
@inline check_address_type(h, addr) = check_address_type(h, typeof(addr))
@inline function check_address_type(h::AbstractObservable, v::Union{AbstractArray,Tuple})
    all(check_address_type(h, a) for a in v)
end

(h::AbstractObservable)(v) = h * v
(h::AbstractObservable)(w, v) = mul!(w, h, v)

BitStringAddresses.num_modes(h::AbstractHamiltonian) = num_modes(starting_address(h))

"""
    dimension(h::AbstractHamiltonian, addr=starting_address(h))
    dimension(h::AbstractObservable, addr)
    dimension(addr::AbstractFockAddress)
    dimension(::Type{<:AbstractFockAddress})

Return the estimated dimension of Hilbert space. May return a `BigInt` number.

When called on an address or address type, the dimension of the linear space spanned by the
address type is returned. When called on an `AbstractHamiltonian`, an upper bound on the
dimension of the matrix representing the Hamiltonian is returned.

# Examples

```jldoctest
julia> dimension(OccupationNumberFS(1,2,3))
16777216

julia> dimension(HubbardReal1D(OccupationNumberFS(1,2,3)))
28

julia> dimension(BoseFS{200,100})
1386083821086188248261127842108801860093488668581216236221011219101585442774669540

julia> Float64(ans)
1.3860838210861882e81
```

Part of the [`AbstractHamiltonian`](@ref) interface. See also
[`BasisSetRepresentation`](@ref).
# Extended Help

The default fallback for `dimension` called on an [`AbstractHamiltonian`](@ref) is to return
the dimension of the address space, which provides an upper bound. For new Hamiltonians a
tighter bound can be provided by defining a custom method.

When extending [`AbstractHamiltonian`](@ref), define a method for the two-argument form
`dimension(h::MyNewHamiltonian, addr)`. For number-conserving Hamiltonians, the function
[`Hamiltonians.number_conserving_dimension`](@ref) may be useful.

When extending [`AbstractFockAddress`](@ref), define a method for
`dimension(::Type{MyNewFockAddress})`.
"""
dimension(h::AbstractHamiltonian) = dimension(h, starting_address(h))
dimension(::AbstractObservable, addr) = dimension(addr)
dimension(addr::AbstractFockAddress) = dimension(typeof(addr))
dimension(::T) where {T<:Number} = typemax(T) # e.g. integer addresses
function dimension(m::AbstractMatrix)
    r,c = size(m)
    if r == c
        return r
    else
        throw(ArgumentError("Matrix is not square"))
    end
end

function dimension(::Type{<:BoseFS{N,M}}) where {N,M}
    return number_conserving_bose_dimension(N,M)
end
function dimension(::Type{<:OccupationNumberFS{M,T}}) where {M,T}
    n = typemax(T)
    return BigInt(n + 1)^BigInt(M)
end
function dimension(::Type{<:FermiFS{N,M}}) where {N,M}
    return number_conserving_fermi_dimension(N, M)
end
function dimension(::Type{<:CompositeFS{<:Any,<:Any,<:Any,T}}) where {T}
    return prod(dimension, T.parameters)
    # This relies on an implementation detail of the Tuple type and may break in future
    # julia versions.
end

# for backward compatibility
function dimension(::Type{T}, h, addr=starting_address(h)) where {T}
    return T(dimension(h, addr))
end
dimension(::Type{T}, addr::AbstractFockAddress) where {T} = T(dimension(addr))

Base.isreal(h::AbstractObservable) = eltype(h) <: Real
LinearAlgebra.isdiag(h::AbstractObservable) = LOStructure(h) ≡ IsDiagonal()
LinearAlgebra.ishermitian(h::AbstractObservable) = LOStructure(h) ≡ IsHermitian()
LinearAlgebra.issymmetric(h::AbstractObservable) = ishermitian(h) && isreal(h)

BitStringAddresses.near_uniform(h::AbstractHamiltonian) = near_uniform(typeof(starting_address(h)))

"""
    number_conserving_bose_dimension(n, m)

Return the dimension of the number-conserving Fock space for `n` bosons in `m` modes:
`binomial(n + m - 1, n)`.

See also [`number_conserving_fermi_dimension`](@ref), [`number_conserving_dimension`](@ref).
"""
number_conserving_bose_dimension(n, m) = binomial(BigInt(n + m - 1), BigInt(n))
"""
    number_conserving_fermi_dimension(n, m)

Return the dimension of the number-conserving Fock space for `n` fermions in `m` modes:
`binomial(m, n)`.

See also [`number_conserving_bose_dimension`](@ref), [`number_conserving_dimension`](@ref).
"""
number_conserving_fermi_dimension(n, m) = binomial(BigInt(m), BigInt(n))

"""
    number_conserving_dimension(address <: AbstractFockAddress)

Return the dimension of the Fock space spanned by the address type assuming particle number
conservation.

See also [`number_conserving_bose_dimension`](@ref),
[`number_conserving_fermi_dimension`](@ref), [`dimension`](@ref).
"""
function number_conserving_dimension(address::Union{BoseFS,OccupationNumberFS})
    m = num_modes(address)
    n = num_particles(address)
    return number_conserving_bose_dimension(n, m)
end
function number_conserving_dimension(address::FermiFS)
    m = num_modes(address)
    n = num_particles(address)
    return number_conserving_fermi_dimension(n, m)
end
function number_conserving_dimension(address::CompositeFS)
    return prod(number_conserving_dimension, address.components)
end

"""
    rayleigh_quotient(H, v)

Return the Rayleigh quotient of the linear operator `H` and the vector `v`:

```math
\\frac{⟨ v | H | v ⟩}{⟨ v|v ⟩}
```
"""
rayleigh_quotient(lo, v) = dot(v, lo, v) / norm(v)^2

"""
    TwoComponentHamiltonian{T} <: AbstractHamiltonian{T}

Abstract type for representing interacting two-component Hamiltonians in a Fock space with
two different species. At least the following fields should be present:

* `ha` Hamiltonian for species A
* `hb` Hamiltonian for species B

See [`AbstractHamiltonian`](@ref) for a list of methods that need to be defined.

Provides an implementation of [`dimension`](@ref).
"""
abstract type TwoComponentHamiltonian{T} <: AbstractHamiltonian{T} end

function dimension(h::TwoComponentHamiltonian)
    return dimension(h.ha) * dimension(h.hb)
end

"""
    momentum(ham::AbstractHamiltonian)

Momentum as a linear operator in Fock space. Pass a Hamiltonian `ham` in order to convey
information about the Fock basis. Returns an [`AbstractHamiltonian`](@ref) that represents
the momentum operator.

Note: `momentum` is currently only defined on [`HubbardMom1D`](@ref).

# Example

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


julia> ham = HubbardMom1D(add; u = 2.0, t = 1.0);


julia> mom = momentum(ham);


julia> diagonal_element(mom, add) # calculate the momentum of a single configuration
-1.5707963267948966

julia> v = DVec(add => 10; capacity=1000);


julia> rayleigh_quotient(mom, v) # momentum expectation value for state vector `v`
-1.5707963267948966
```
Part of the [`AbstractHamiltonian`](@ref) interface.
"""
momentum

function Base.getindex(ham::AbstractObservable{T}, address1, address2) where {T}
    # calculate the matrix element when only two bitstring addresses are given
    # this is NOT used for the QMC algorithm and is currently not used either
    # for building the matrix for conventional diagonalisation.
    # Only used for verifying matrix.
    # This will be slow and inefficient. Avoid using for larger Hamiltonians!
    address1 == address2 && return diagonal_element(ham, address1) # diagonal
    res = zero(T)
    for (add, val) in offdiagonals(ham, address2) # off-diag column as iterator
        if add == address1
            res += val # found address1
        end
    end
    return res
end

LinearAlgebra.adjoint(op::AbstractObservable) = adjoint(LOStructure(op), op)

"""
    adjoint(::LOStructure, op::AbstractObservable)

Represent the adjoint of an [`AbstractObservable`](@ref). Extend this method to define
custom adjoints.
"""
function LinearAlgebra.adjoint(::S, op) where {S<:LOStructure}
    throw(ArgumentError(
        "`adjoint()` is not defined for `AbstractObservable`s with `LOStructure` `$(S)`. "*
        " Is your Hamiltonian hermitian?"
    ))
end

LinearAlgebra.adjoint(::IsHermitian, op) = op # adjoint is known
function LinearAlgebra.adjoint(::IsDiagonal, op)
    if eltype(op) <: Real || eltype(eltype(op)) <: Real # for op's that return vectors
        return op
    else
        throw(ArgumentError("adjoint() is not implemented for complex diagonal Hamiltonians"))
    end
end

"""
    TransformUndoer(transform::AbstractHamiltonian, op::AbstractObservable) <: AbstractHamiltonian

Create a new operator for the purpose of calculating overlaps of transformed
vectors, which are defined by some transformation `transform`. The new operator should
represent the effect of undoing the transformation before calculating overlaps, including
with an optional operator `op`.

Not exported; transformations should define all necessary methods and properties,
see [`AbstractHamiltonian`](@ref). An `ArgumentError` is thrown if used with an
unsupported transformation.

# Example

A similarity transform ``\\hat{G} = f \\hat{H} f^{-1}`` has eigenvector
``d = f \\cdot c`` where ``c`` is an eigenvector of ``\\hat{H}``. Then the
overlap ``c' \\cdot c = d' \\cdot f^{-2} \\cdot d`` can be computed by defining all
necessary methods for `TransformUndoer(G)` to represent the operator ``f^{-2}`` and
calculating `dot(d, TransformUndoer(G), d)`.

Observables in the transformed basis can be computed by defining `TransformUndoer(G, A)`
to represent ``f^{-1} A f^{-1}``.

# Supported transformations

* [`GutzwillerSampling`](@ref)
* [`GuidingVectorSampling`](@ref)
"""
struct TransformUndoer{
    T,K<:AbstractHamiltonian,O<:Union{AbstractObservable,Nothing}
} <: AbstractHamiltonian{T}
    transform::K
    op::O
end

function TransformUndoer(k::AbstractHamiltonian, op)
    # default check
    throw(ArgumentError("Unsupported transformation: $k"))
end
TransformUndoer(k::AbstractHamiltonian) = TransformUndoer(k::AbstractHamiltonian, nothing)

# common methods
starting_address(s::TransformUndoer) = starting_address(s.transform)
dimension(s::TransformUndoer, addr) = dimension(s.transform, addr)
function Base.:(==)(a::TransformUndoer, b::TransformUndoer)
    return a.transform == b.transform && a.op == b.op
end

1	"""
2	check_address_type(h::AbstractObservable, addr_or_type)
3	Throw an `ArgumentError` if `addr_or_type` is not compatible with `h`, otherwise return
4	`true`. Acceptable arguments are either an address or an address type, or a tuple or
5	array thereof.
6
7	See also [`allows_address_type`](@ref).
8	"""
9	@inline function check_address_type(h, ::Type{A}) where {A}	1,118✔
10	if !allows_address_type(h, A)	1,134✔
11	throw(ArgumentError("address type $A not allowed for operator $h"))	2✔
12	end
13	return true	1,116✔
14	end
15	@inline check_address_type(h, addr) = check_address_type(h, typeof(addr))	467✔
16	@inline function check_address_type(h::AbstractObservable, v::Union{AbstractArray,Tuple})	81✔
17	all(check_address_type(h, a) for a in v)	195✔
18	end
19
20	(h::AbstractObservable)(v) = h * v	1,775✔
21	(h::AbstractObservable)(w, v) = mul!(w, h, v)	39✔
22
23	BitStringAddresses.num_modes(h::AbstractHamiltonian) = num_modes(starting_address(h))	24,972,979✔
24
25	"""
26	dimension(h::AbstractHamiltonian, addr=starting_address(h))
27	dimension(h::AbstractObservable, addr)
28	dimension(addr::AbstractFockAddress)
29	dimension(::Type{<:AbstractFockAddress})
30
31	Return the estimated dimension of Hilbert space. May return a `BigInt` number.
32
33	When called on an address or address type, the dimension of the linear space spanned by the
34	address type is returned. When called on an `AbstractHamiltonian`, an upper bound on the
35	dimension of the matrix representing the Hamiltonian is returned.
36
37	# Examples
38
39	```jldoctest
40	julia> dimension(OccupationNumberFS(1,2,3))
41	16777216
42
43	julia> dimension(HubbardReal1D(OccupationNumberFS(1,2,3)))
44	28
45
46	julia> dimension(BoseFS{200,100})
47	1386083821086188248261127842108801860093488668581216236221011219101585442774669540
48
49	julia> Float64(ans)
50	1.3860838210861882e81
51	```
52
53	Part of the [`AbstractHamiltonian`](@ref) interface. See also
54	[`BasisSetRepresentation`](@ref).
55	# Extended Help
56
57	The default fallback for `dimension` called on an [`AbstractHamiltonian`](@ref) is to return
58	the dimension of the address space, which provides an upper bound. For new Hamiltonians a
59	tighter bound can be provided by defining a custom method.
60
61	When extending [`AbstractHamiltonian`](@ref), define a method for the two-argument form
62	`dimension(h::MyNewHamiltonian, addr)`. For number-conserving Hamiltonians, the function
63	[`Hamiltonians.number_conserving_dimension`](@ref) may be useful.
64
65	When extending [`AbstractFockAddress`](@ref), define a method for
66	`dimension(::Type{MyNewFockAddress})`.
67	"""
68	dimension(h::AbstractHamiltonian) = dimension(h, starting_address(h))	57,114✔
69	dimension(::AbstractObservable, addr) = dimension(addr)	181✔
70	dimension(addr::AbstractFockAddress) = dimension(typeof(addr))	235✔
71	dimension(::T) where {T<:Number} = typemax(T) # e.g. integer addresses	3✔
72	function dimension(m::AbstractMatrix)	10✔
73	r,c = size(m)	10✔
74	if r == c	10✔
75	return r	10✔
76	else
NEW 77	throw(ArgumentError("Matrix is not square"))	×
78	end
79	end
80
81	function dimension(::Type{<:BoseFS{N,M}}) where {N,M}	221✔
82	return number_conserving_bose_dimension(N,M)	221✔
83	end
84	function dimension(::Type{<:OccupationNumberFS{M,T}}) where {M,T}	5✔
85	n = typemax(T)	5✔
86	return BigInt(n + 1)^BigInt(M)	5✔
87	end
88	function dimension(::Type{<:FermiFS{N,M}}) where {N,M}	31✔
89	return number_conserving_fermi_dimension(N, M)	31✔
90	end
91	function dimension(::Type{<:CompositeFS{<:Any,<:Any,<:Any,T}}) where {T}	16✔
92	return prod(dimension, T.parameters)	16✔
93	# This relies on an implementation detail of the Tuple type and may break in future
94	# julia versions.
95	end
96
97	# for backward compatibility
98	function dimension(::Type{T}, h, addr=starting_address(h)) where {T}	×
99	return T(dimension(h, addr))	×
100	end
101	dimension(::Type{T}, addr::AbstractFockAddress) where {T} = T(dimension(addr))	1✔
102
103	Base.isreal(h::AbstractObservable) = eltype(h) <: Real	109✔
104	LinearAlgebra.isdiag(h::AbstractObservable) = LOStructure(h) ≡ IsDiagonal()	505,921✔
105	LinearAlgebra.ishermitian(h::AbstractObservable) = LOStructure(h) ≡ IsHermitian()	490✔
106	LinearAlgebra.issymmetric(h::AbstractObservable) = ishermitian(h) && isreal(h)	136✔
107
108	BitStringAddresses.near_uniform(h::AbstractHamiltonian) = near_uniform(typeof(starting_address(h)))	1✔
109
110	"""
111	number_conserving_bose_dimension(n, m)
112
113	Return the dimension of the number-conserving Fock space for `n` bosons in `m` modes:
114	`binomial(n + m - 1, n)`.
115
116	See also [`number_conserving_fermi_dimension`](@ref), [`number_conserving_dimension`](@ref).
117	"""
118	number_conserving_bose_dimension(n, m) = binomial(BigInt(n + m - 1), BigInt(n))	647✔
119	"""
120	number_conserving_fermi_dimension(n, m)
121
122	Return the dimension of the number-conserving Fock space for `n` fermions in `m` modes:
123	`binomial(m, n)`.
124
125	See also [`number_conserving_bose_dimension`](@ref), [`number_conserving_dimension`](@ref).
126	"""
127	number_conserving_fermi_dimension(n, m) = binomial(BigInt(m), BigInt(n))	694✔
128
129	"""
130	number_conserving_dimension(address <: AbstractFockAddress)
131
132	Return the dimension of the Fock space spanned by the address type assuming particle number
133	conservation.
134
135	See also [`number_conserving_bose_dimension`](@ref),
136	[`number_conserving_fermi_dimension`](@ref), [`dimension`](@ref).
137	"""
138	function number_conserving_dimension(address::Union{BoseFS,OccupationNumberFS})	426✔
139	m = num_modes(address)	426✔
140	n = num_particles(address)	426✔
141	return number_conserving_bose_dimension(n, m)	426✔
142	end
143	function number_conserving_dimension(address::FermiFS)	663✔
144	m = num_modes(address)	663✔
145	n = num_particles(address)	663✔
146	return number_conserving_fermi_dimension(n, m)	663✔
147	end
148	function number_conserving_dimension(address::CompositeFS)	230✔
149	return prod(number_conserving_dimension, address.components)	230✔
150	end
151
152	"""
153	rayleigh_quotient(H, v)
154
155	Return the Rayleigh quotient of the linear operator `H` and the vector `v`:
156
157	```math
158	\\frac{⟨ v \| H \| v ⟩}{⟨ v\|v ⟩}
159	```
160	"""
161	rayleigh_quotient(lo, v) = dot(v, lo, v) / norm(v)^2	6✔
162
163	"""
164	TwoComponentHamiltonian{T} <: AbstractHamiltonian{T}
165
166	Abstract type for representing interacting two-component Hamiltonians in a Fock space with
167	two different species. At least the following fields should be present:
168
169	* `ha` Hamiltonian for species A
170	* `hb` Hamiltonian for species B
171
172	See [`AbstractHamiltonian`](@ref) for a list of methods that need to be defined.
173
174	Provides an implementation of [`dimension`](@ref).
175	"""
176	abstract type TwoComponentHamiltonian{T} <: AbstractHamiltonian{T} end
177
178	function dimension(h::TwoComponentHamiltonian)	×
179	return dimension(h.ha) * dimension(h.hb)	×
180	end
181
182	"""
183	momentum(ham::AbstractHamiltonian)
184
185	Momentum as a linear operator in Fock space. Pass a Hamiltonian `ham` in order to convey
186	information about the Fock basis. Returns an [`AbstractHamiltonian`](@ref) that represents
187	the momentum operator.
188
189	Note: `momentum` is currently only defined on [`HubbardMom1D`](@ref).
190
191	# Example
192
193	```jldoctest
194	julia> add = BoseFS((1, 0, 2, 1, 2, 1, 1, 3));
195
196
197	julia> ham = HubbardMom1D(add; u = 2.0, t = 1.0);
198
199
200	julia> mom = momentum(ham);
201
202
203	julia> diagonal_element(mom, add) # calculate the momentum of a single configuration
204	-1.5707963267948966
205
206	julia> v = DVec(add => 10; capacity=1000);
207
208
209	julia> rayleigh_quotient(mom, v) # momentum expectation value for state vector `v`
210	-1.5707963267948966
211	```
212	Part of the [`AbstractHamiltonian`](@ref) interface.
213	"""
214	momentum
215
216	function Base.getindex(ham::AbstractObservable{T}, address1, address2) where {T}	121✔
217	# calculate the matrix element when only two bitstring addresses are given
218	# this is NOT used for the QMC algorithm and is currently not used either
219	# for building the matrix for conventional diagonalisation.
220	# Only used for verifying matrix.
221	# This will be slow and inefficient. Avoid using for larger Hamiltonians!
222	address1 == address2 && return diagonal_element(ham, address1) # diagonal	121✔
223	res = zero(T)	110✔
224	for (add, val) in offdiagonals(ham, address2) # off-diag column as iterator	220✔
225	if add == address1	400✔
226	res += val # found address1	40✔
227	end
228	end	690✔
229	return res	110✔
230	end
231
232	LinearAlgebra.adjoint(op::AbstractObservable) = adjoint(LOStructure(op), op)	192✔
233
234	"""
235	adjoint(::LOStructure, op::AbstractObservable)
236
237	Represent the adjoint of an [`AbstractObservable`](@ref). Extend this method to define
238	custom adjoints.
239	"""
240	function LinearAlgebra.adjoint(::S, op) where {S<:LOStructure}	11✔
241	throw(ArgumentError(	11✔
242	"`adjoint()` is not defined for `AbstractObservable`s with `LOStructure` `$(S)`. "*
243	" Is your Hamiltonian hermitian?"
244	))
245	end
246
247	LinearAlgebra.adjoint(::IsHermitian, op) = op # adjoint is known	153✔
248	function LinearAlgebra.adjoint(::IsDiagonal, op)	28✔
249	if eltype(op) <: Real \|\| eltype(eltype(op)) <: Real # for op's that return vectors	28✔
250	return op	28✔
251	else
252	throw(ArgumentError("adjoint() is not implemented for complex diagonal Hamiltonians"))	×
253	end
254	end
255
256	"""
257	TransformUndoer(transform::AbstractHamiltonian, op::AbstractObservable) <: AbstractHamiltonian
258
259	Create a new operator for the purpose of calculating overlaps of transformed
260	vectors, which are defined by some transformation `transform`. The new operator should
261	represent the effect of undoing the transformation before calculating overlaps, including
262	with an optional operator `op`.
263
264	Not exported; transformations should define all necessary methods and properties,
265	see [`AbstractHamiltonian`](@ref). An `ArgumentError` is thrown if used with an
266	unsupported transformation.
267
268	# Example
269
270	A similarity transform ``\\hat{G} = f \\hat{H} f^{-1}`` has eigenvector
271	``d = f \\cdot c`` where ``c`` is an eigenvector of ``\\hat{H}``. Then the
272	overlap ``c' \\cdot c = d' \\cdot f^{-2} \\cdot d`` can be computed by defining all
273	necessary methods for `TransformUndoer(G)` to represent the operator ``f^{-2}`` and
274	calculating `dot(d, TransformUndoer(G), d)`.
275
276	Observables in the transformed basis can be computed by defining `TransformUndoer(G, A)`
277	to represent ``f^{-1} A f^{-1}``.
278
279	# Supported transformations
280
281	* [`GutzwillerSampling`](@ref)
282	* [`GuidingVectorSampling`](@ref)
283	"""
284	struct TransformUndoer{
285	T,K<:AbstractHamiltonian,O<:Union{AbstractObservable,Nothing}
286	} <: AbstractHamiltonian{T}
287	transform::K	111✔
288	op::O
289	end
290
291	function TransformUndoer(k::AbstractHamiltonian, op)	6✔
292	# default check
293	throw(ArgumentError("Unsupported transformation: $k"))	6✔
294	end
295	TransformUndoer(k::AbstractHamiltonian) = TransformUndoer(k::AbstractHamiltonian, nothing)	24✔
296
297	# common methods
298	starting_address(s::TransformUndoer) = starting_address(s.transform)	7✔
299	dimension(s::TransformUndoer, addr) = dimension(s.transform, addr)	3✔
300	function Base.:(==)(a::TransformUndoer, b::TransformUndoer)	1✔
301	return a.transform == b.transform && a.op == b.op	1✔
302	end

RimuQMC / Rimu.jl / 14949967345

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous