Example 6: Molecular Hamiltonian Calculation

This is an example calculation of the water molecule (H₂O). It requires ElemCo.jl to be installed. To install it, in Julia REPL type ] to enter the Pkg REPL mode and run:

pkg> add ElemCo

Or via the Pkg package

julia> import Pkg; Pkg.add("ElemCo")

A runnable script for this example is located here. Run it with julia Molecule-example.jl.

First, we load the required packages. Rimu for the calculation and ElemCo for loading the FCIDUMP file.

using Rimu
using ElemCo

━━━━━━━━━━━━━━━
   ElemCo.jl
━━━━━━━━━━━━━━━
Version: 0.14.1
Git hash: unknown
Website: elem.co.il
Julia version: 1.11.7
BLAS threads: 2
OpenMP threads: 1
Hostname: runnervmf4ws1
Scratch directory: /tmp
Date: 2025-09-21 14:31:45
╭──────────────────────────────╮ 
│        ╭─────────────╮       ├─╮
│ Electron Correlation methods │ │
│        ╰─────────────╯       │ │
╰─┬────────────────────────────╯ │
  ╰──────────────────────────────╯

Setting up the model

We specify the path to the FCIDUMP file describing the H₂O molecule.

fcidump = joinpath(pkgdir(Rimu), "test/examples/h2o.FCIDUMP");

Next we construct the Hamiltonian by the constructor MolecularHamiltonian. The Hartree-Fock ground state is generated automatically at the same time.

h = MolecularHamiltonian(fcidump)
a = starting_address(h);

Running the calculation

Exact Diagonalization calculation

We first define the problem via ExactDiagonalizationProblem interface. Since H₂O has a rather large Hilbert space, a iterative solver is recomended. We use KrylovKitSolver in this example. It requires KrylovKit.jl to be installed. To install it, in Julia REPL type ] to enter the Pkg REPL mode and run:

pkg> add KrylovKit

Or via the Pkg package

julia> import Pkg; Pkg.add("KrylovKit")

using KrylovKit

p = ExactDiagonalizationProblem(h, algorithm=KrylovKitSolver(true))

ExactDiagonalizationProblem(
  MolecularHamiltonian("/home/runner/work/Rimu.jl/Rimu.jl/test/examples/h2o.FCIDUMP", starting_address=fs"|⇅⇅⇅⇅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⟩"),
  nothing;
  (algorithm = KrylovKitSolver(matrix_free = true,),)...
)

Then it can be solved by executing

s = solve(p)

Projector Monte Carlo / FCIQMC

Set up the FCIQMC parameters.

steps_equilibrate = 1_000
steps_measure = 2_000
target_walkers = 1_000
time_step = 0.001;

Define the problem via ProjectorMonteCarloProblem interface.

p = ProjectorMonteCarloProblem(h;
    time_step,
    last_step = steps_equilibrate + steps_measure,
    target_walkers,
    initiator=true,
)

ProjectorMonteCarloProblem with 1 replica(s) and 1 spectral state(s):
  algorithm = FCIQMC(DoubleLogUpdate{Int64}(1000, 0.08, 0.0016), ConstantTimeStep())
  hamiltonian = MolecularHamiltonian("/home/runner/work/Rimu.jl/Rimu.jl/test/examples/h2o.FCIDUMP", starting_address=fs"|⇅⇅⇅⇅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⟩")
  style = IsDynamicSemistochastic{Float64,ThresholdCompression,DynamicSemistochastic}()
  initiator = Initiator{Float64}(1.0)
  threading = false
  simulation_plan = SimulationPlan(starting_step=0, last_step=3000, wall_time=Inf)
  replica_strategy = NoStats{1}()
  reporting_strategy = ReportDFAndInfo
  reporting_interval: Int64 1
  info_interval: Int64 100
  io: IOContext{Base.PipeEndpoint}
  writeinfo: Bool false
  post_step_strategy = ()
  spectral_strategy = GramSchmidt(1; orthogonalization_interval=1)
  max_length = 2100
  metadata = OrderedCollections.LittleDict("display_name" => "PMCSimulation")
  random_seed = 8848213650822725940

Run the calculation.

result = solve(p);

Store the result into DataFrame.

df = DataFrame(result);

To analyse the energy shift, we can use shift_estimator.

se = shift_estimator(df; skip=steps_equilibrate)

BlockingResult{Float64}
  mean = -76.128 ± 0.065
  with uncertainty of ± 0.005900556412098986
  from 62 blocks after 5 transformations (k = 6).

This page was generated using Literate.jl.