Monte Carlo Sampling#
sample draws configurations from \(|\psi(x)|^2\) with the Metropolis–Hastings
algorithm. At each step, a move proposes a new configuration \(x'\) on every
chain, and the chain accepts it with probability
where \(q(x' \mid x)\) is the probability that the move proposes \(x'\) from \(x\).
from tachys.montecarlo import sample
key, subkey = jax.random.split(key)
mc_keys = jax.random.split(subkey, N_mc) # one key per chain
state, log_amps, acceptance = sample(nsweeps, state, action, mc_keys, wf)
A sweep is Ns steps, so each chain receives nsweeps × Ns proposals. sample
returns the new configurations, their log-amplitudes, which
compute_expectation reuses, and the fraction of accepted proposals.
Choosing a move#
A move must let the chains reach every configuration of the sector, and keep the quantum numbers the Hamiltonian conserves — otherwise the chains leave the sector set by the initial configurations.
Move |
Configurations |
Proposal |
Conserves |
|---|---|---|---|
|
spins |
swap two antiparallel spins on a bond |
\(S^z\) |
|
fermions |
move an electron along a bond, keeping its spin |
\(N_\uparrow\), \(N_\downarrow\) |
|
fermions |
swap the spins of two singly occupied sites on a bond |
\(N_\uparrow\), \(N_\downarrow\) |
|
spins |
flip one spin |
— |
|
spins |
flip both spins on a bond |
parity of \(N_\uparrow\) |
BondExchange lives in tachys.lattice.bond_exchange, SpinFlip and
BondFlip in tachys.lattice.spins.spin_action, FermionSpinExchange in
tachys.lattice.fermions.fermion_action.
BondExchange and FermionSpinExchange draw their bonds from the distance
shells of the lattice (Lattices): max_dist=1, the default, uses
nearest neighbours, and max_dist=2 adds the next shell, which helps when the
Hamiltonian couples further neighbours. BondFlip takes cell displacements, as
lattice.bonds does.
For fermions at strong coupling and close to half filling, an electron that
moves usually creates a doubly occupied site and is rejected, so the spins
rearrange slowly.
FermionSpinExchange exchanges spins without moving charge. CompositeAction
mixes moves, choosing one on each chain at each step with the given
probabilities:
from tachys.montecarlo import CompositeAction
from tachys.lattice.bond_exchange import BondExchange
from tachys.lattice.fermions.fermion_action import FermionSpinExchange
action = CompositeAction(
actions=(
BondExchange.create(lattice, Nbands=2),
FermionSpinExchange.create(lattice),
),
probs=(0.5, 0.5),
)
With a composite move, acceptance has one entry per move.
Sweeps, burn-in and acceptance#
Successive configurations of a chain are correlated, and nsweeps sets how
many proposals separate two samples. The chains carry over from one call to the
next, so a single sweep per optimization step is often enough; the quickstart
uses nsweeps=1.
The initial configurations are not distributed according to \(|\psi|^2\). Before measuring anything, and at the start of training, a burn-in call brings the chains to equilibrium:
key, subkey = jax.random.split(key)
mc_keys = jax.random.split(subkey, N_mc)
state, log_amps, acceptance = sample(20, state, action, mc_keys, wf)
A low acceptance means the chains move slowly: increase nsweeps, or choose a
move better suited to the wavefunction.
Writing a move#
A move is a subclass of _BaseAction (tachys.montecarlo). Its
__call__(key, state) receives one PRNG key per chain and the batch of
configurations, and returns three values:
new_state, the proposed configurations;allowed_move, one boolean per chain:Falserejects the proposal without evaluating the wavefunction, for proposals that leave the configuration unchanged or violate a constraint;log_prob_correction, the logarithm of \(q(x \mid x')/q(x' \mid x)\), or0.0for a symmetric proposal.
This move swaps two sites drawn anywhere on the lattice. It conserves \(S^z\),
like BondExchange, but also connects distant sites in a single step:
import jax
import jax.numpy as jnp
from tachys.montecarlo import _BaseAction
class RandomExchange(_BaseAction):
def __call__(self, key, state):
spins = state.spins # (N_mc, Ns)
chains = jnp.arange(spins.shape[0])
Ns = spins.shape[1]
pick = lambda k: jax.random.choice(k, Ns, (2,), replace=False)
i, j = jax.vmap(pick)(key).T # two distinct sites
si, sj = spins[chains, i], spins[chains, j]
new_spins = spins.at[chains, i].set(sj).at[chains, j].set(si)
allowed = si != sj # swapping equal spins changes nothing
return state.replace(spins=new_spins), allowed, 0.0
action = RandomExchange()
The pair is drawn uniformly, so the proposal is symmetric and the correction vanishes.