Excited states #269
Excited states #269
Conversation
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
Architecture: x86_64
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
Architecture: x86_64
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
Architecture: x86_64
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
Architecture: x86_64
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
Architecture: x86_64
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rimu.jl/Rimu.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
Architecture: x86_64 |
There was a problem hiding this comment.
Thanks for the nice PR!
I like the idea of orthogonalising the state before the spawning step. It makes life easier for the norm control machinery. It does mean, however, that we cannot rely on the instantaneous states being exactly orthogonal.
It would be great to have this work smoothly with ProjectorMonteCarloProblem and init! to create an appropriate set of starting vectors. I think currently the number of permissible spectral states is still hardcoded to be one (see line 39 of "pmc_simulation.jl").
It would further be good to make it easy to pass suitable initial state. It would also be super convenient to have suitable initial state generated automatically if no starting vectors (or an insufficient number of them) are provided.
A possible strategy is the following: Say only a single address is available (e.g. from starting_address(ham). Generate further addresses with offdiagonals until you have at least num_spectral_states of them and then sort them by the value of diagonal_element. Set as starting vectors vectors with one of these addresses each. This should work well for diagonally dominated Hamiltonians.
Yes it would be good to be able to automatically generate starting vectors, thanks, I'll look into doing that. |
Rewrite orthogonalization of vectors, make orthogonalization_interval a keyword argument, fix spelling Co-authored-by: Joachim Brand <j.brand@massey.ac.nz> Co-authored-by: mtsch <matijacufar@gmail.com>
Pull Request Test Coverage Report for Build 11991921089Details
💛 - Coveralls |
There was a problem hiding this comment.
Looks great! The next step would be to add some tests. It doesn't have to be anything complicated, just something that makes sure the code works correctly. I suggest something like setting up a small Hamiltonian (maybe HubbardReal1D(BoseFS(1,1,1,1,1))), solving it with FCIQMC with the style set to IsDeterministic(), then comparing the resulting energies and vectors with the result that comes out of solveing an ExactDiagonalizationProblem.
If you run into issues getting this to work feel free to ask here or on Zulip!
src/pmc_simulation.jl
Outdated
| basis = build_basis(ham; cutoff=diagonal_element(ham, starting_address(ham)) + 5) | ||
| basis = sort!(basis; by=x -> diagonal_element(ham, x))[1:2*n_spectral] | ||
| trunc = BasisSetRepresentation(ham, basis; filter=Returns(false), sizelim=Inf) | ||
| vecs = eigvecs(Matrix(trunc)) | ||
| if threading | ||
| eigs = [PDVec(zip(basis, 10*vecs[:,i]); style, initiator) for i in 1:n_spectral] | ||
| else | ||
| if initiator == NonInitiator() | ||
| eigs = [DVec(zip(basis, 10*vecs[:,i]); style) for i in 1:n_spectral] | ||
| else | ||
| eigs = [InitiatorDVec(zip(basis,10*vecs[:,i]);style,initiator) for i in 1:n_spectral] | ||
| end | ||
| end |
There was a problem hiding this comment.
I like this idea a lot, but we should be careful. It's possible to setup a Hamiltonian that's huge, but still has all diagonal elements below diagonal_element(ham, starting_address(ham)) + 5. Running this on it will not work.
There was a problem hiding this comment.
As discussed, I have the same concern. I would suggest to not use an energy cutoff, but instead implement a new option for build_basis where you specify a lower bound, e.g. called mininum_size, that allows build_basis to stop generating new addresses once the mininum_size has been reached.
This will require modifying the code of build_basis, but should be fairly easy to implement.
There was a problem hiding this comment.
I've tried to implement this and it seems to work for me, but one of the tests failed and I'm not understanding what happened - could somebody please take a look?
There was a problem hiding this comment.
The failing test has nothing to do with your (or our) code. A different package which is a dependency of our code (here REPL.jl) failed to precompile on the nightly version of Julia. This can happen when things are changed in the bleeding-edge (and unreleased) development version of Julia. Usually these problems will disappear by themselves.
There was a problem hiding this comment.
Very nice PR! I've found two issues:
- The new functionality is not yet well documented. I suggest documenting the option of passing a
spectral_strategyin the docstring ofProjectorMonteCarloProblem, and adding a reference toGramSchmidtin the docstring ofSpectralStrategy. I'm not completely sure whether the latter should be exported. The used would not need to use it directly in most cases, but I'd like to properly document the available options (which currently is onlyGramSchmidt. - I've checked that this works nicely when generating excited state with a single replica. However, combining replicas with excited states throws an error:
julia> h = HubbardReal1D(near_uniform(BoseFS{10,10}))
HubbardReal1D(fs"|1 1 1 1 1 1 1 1 1 1⟩"; u=1.0, t=1.0)
julia> spectral_strategy = GramSchmidt(3)
GramSchmidt(3; orthogonalization_interval=1)
julia> p = ProjectorMonteCarloProblem(h; spectral_strategy, replica_strategy=AllOverlaps(2), last_step=2_000)
ProjectorMonteCarloProblem with 2 replica(s) and 3 spectral state(s):
Algorithm: FCIQMC(DoubleLogUpdate{Int64}(1000, 0.08, 0.0016), ConstantTimeStep())
Hamiltonian: HubbardReal1D(fs"|1 1 1 1 1 1 1 1 1 1⟩"; u=1.0, t=1.0)
Style: IsDynamicSemistochastic{Float64,ThresholdCompression,DynamicSemistochastic}()
Initiator: NonInitiator()
Threading: false
Simulation Plan: SimulationPlan(starting_step=0, last_step=2000, walltime=Inf)
Replica Strategy: AllOverlaps{2, 0, Tuple{}, true}(())
Reporting Strategy: ReportDFAndInfo
reporting_interval: Int64 1
info_interval: Int64 100
io: Base.TTY
writeinfo: Bool false
Post Step Strategy: ()
Spectral Strategy: GramSchmidt(3; orthogonalization_interval=1)
Maxlength: 2100
Metadata: OrderedCollections.LittleDict("display_name" => "PMCSimulation")
Random Seed: 1943348664599472552
julia> sim = solve(p)
ERROR: ArgumentError: Collection has multiple elements, must contain exactly 1 element
Stacktrace:
[1] only
@ ./iterators.jl:1527 [inlined]
[2] #34
@ ~/git/code/Rimu/src/strategies_and_params/replicastrategy.jl:136 [inlined]
[3] ntuple
@ ./ntuple.jl:49 [inlined]
[4] replica_stats(rs::AllOverlaps{…}, spectral_states::Tuple{…})
@ Rimu ~/git/code/Rimu/src/strategies_and_params/replicastrategy.jl:136
[5] step!(sm::Rimu.PMCSimulation)
@ Rimu ~/git/code/Rimu/src/pmc_simulation.jl:283
[6] macro expansion
@ ~/git/code/Rimu/src/pmc_simulation.jl:428 [inlined]
[7] macro expansion
@ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
[8]
@ Rimu ~/git/code/Rimu/src/pmc_simulation.jl:422
[9] solve!
@ ~/git/code/Rimu/src/pmc_simulation.jl:337 [inlined]
[10] solve(args::ProjectorMonteCarloProblem{2, 3})
@ CommonSolve ~/.julia/packages/CommonSolve/JfpfI/src/CommonSolve.jl:23
[11] top-level scope
@ REPL[77]:1
Some type information was truncated. Use `show(err)` to see complete types.
julia> Co-authored-by: Joachim Brand <joachim.brand@gmail.com>
There was a problem hiding this comment.
Nice work!
See comments below.
To make AllOverlaps work with spectral states requires more thoughts about reporting. I'd be happy to leave that for another PR.
By the way, why were the matrix dimensions swapped (n_replicas and n_spectral)? Was that inconsistent in the earlier code?
src/pmc_simulation.jl
Outdated
| default_starting_vector(sa; style, initiator, threading) for sa in start_at | ||
| ) | ||
| end | ||
| elseif n_spectral > 1 | ||
| basis = build_basis(ham; minimum_size=n_spectral*5) | ||
| basis = sort!(basis; by=x -> diagonal_element(ham, x)) |
There was a problem hiding this comment.
Why are we sorting the basis here? If we later do a full diagonalisation then there is no point, as the eigenvectors returned by eigvecs will be sorted regardless of what the basis looks like.
There was a problem hiding this comment.
Thanks, that's now been removed.
Yes, the dimensions should now be consistent. |
There was a problem hiding this comment.
Very nice!
Can we add a test for the minimum_size option for build_basis?
Using AllOverlaps with more than 1 spectral state is not yet implemented and leads to an unhelpful error message. Maybe a check can be added to the constuctor for ProjectorMonteCarloProblem to throw an error with a helpful message for this combination?
… for excited states Co-authored-by: Joachim Brand <joachim.brand@gmail.com>
…and/Rimu.jl into feature/excited-states
Co-authored-by: Joachim Brand <joachim.brand@gmail.com>
Added a new method for the
advancefunction that takes multiple states and orthogonalises them. Theorthogonalisation_intervalparameter of theGramSchmidtSpectralStrategydetermines how often the states are orthogonalised.