The hard core lattice gas
Hard-core lattice gas problem
Hard-core lattice gas refers to a model used in statistical physics to study the behavior of particles on a lattice, where the particles are subject to an exclusion principle known as the "hard-core" interaction that characterized by a blockade radius. Distances between two particles can not be smaller than this radius.
- Nath T, Rajesh R. Multiple phase transitions in extended hard-core lattice gas models in two dimensions[J]. Physical Review E, 2014, 90(1): 012120.
- Fernandes H C M, Arenzon J J, Levin Y. Monte Carlo simulations of two-dimensional hard core lattice gases[J]. The Journal of chemical physics, 2007, 126(11).
Let define a $10 \times 10$ triangular lattice, with unit vectors
\[\begin{align*} \vec a &= \left(\begin{matrix}1 \\ 0\end{matrix}\right)\\ \vec b &= \left(\begin{matrix}\frac{1}{2} \\ \frac{\sqrt{3}}{2}\end{matrix}\right) \end{align*}\]
a, b = (1, 0), (0.5, 0.5*sqrt(3))
Na, Nb = 10, 10
sites = vec([50 .* (a .* i .+ b .* j) for i=1:Na, j=1:Nb])100-element Vector{Tuple{Float64, Float64}}:
(75.0, 43.30127018922193)
(125.0, 43.30127018922193)
(175.0, 43.30127018922193)
(225.0, 43.30127018922193)
(275.0, 43.30127018922193)
(325.0, 43.30127018922193)
(375.0, 43.30127018922193)
(425.0, 43.30127018922193)
(475.0, 43.30127018922193)
(525.0, 43.30127018922193)
(100.0, 86.60254037844386)
(150.0, 86.60254037844386)
(200.0, 86.60254037844386)
(250.0, 86.60254037844386)
(300.0, 86.60254037844386)
(350.0, 86.60254037844386)
(400.0, 86.60254037844386)
(450.0, 86.60254037844386)
(500.0, 86.60254037844386)
(550.0, 86.60254037844386)
(125.0, 129.9038105676658)
(175.0, 129.9038105676658)
(225.0, 129.9038105676658)
(275.0, 129.9038105676658)
(325.0, 129.9038105676658)
(375.0, 129.9038105676658)
(425.0, 129.9038105676658)
(475.0, 129.9038105676658)
(525.0, 129.9038105676658)
(575.0, 129.9038105676658)
(150.0, 173.20508075688772)
(200.0, 173.20508075688772)
(250.0, 173.20508075688772)
(300.0, 173.20508075688772)
(350.0, 173.20508075688772)
(400.0, 173.20508075688772)
(450.0, 173.20508075688772)
(500.0, 173.20508075688772)
(550.0, 173.20508075688772)
(600.0, 173.20508075688772)
(175.0, 216.50635094610965)
(225.0, 216.50635094610965)
(275.0, 216.50635094610965)
(325.0, 216.50635094610965)
(375.0, 216.50635094610965)
(425.0, 216.50635094610965)
(475.0, 216.50635094610965)
(525.0, 216.50635094610965)
(575.0, 216.50635094610965)
(625.0, 216.50635094610965)
(200.0, 259.8076211353316)
(250.0, 259.8076211353316)
(300.0, 259.8076211353316)
(350.0, 259.8076211353316)
(400.0, 259.8076211353316)
(450.0, 259.8076211353316)
(500.0, 259.8076211353316)
(550.0, 259.8076211353316)
(600.0, 259.8076211353316)
(650.0, 259.8076211353316)
(225.0, 303.1088913245535)
(275.0, 303.1088913245535)
(325.0, 303.1088913245535)
(375.0, 303.1088913245535)
(425.0, 303.1088913245535)
(475.0, 303.1088913245535)
(525.0, 303.1088913245535)
(575.0, 303.1088913245535)
(625.0, 303.1088913245535)
(675.0, 303.1088913245535)
(250.0, 346.41016151377545)
(300.0, 346.41016151377545)
(350.0, 346.41016151377545)
(400.0, 346.41016151377545)
(450.0, 346.41016151377545)
(500.0, 346.41016151377545)
(550.0, 346.41016151377545)
(600.0, 346.41016151377545)
(650.0, 346.41016151377545)
(700.0, 346.41016151377545)
(275.0, 389.71143170299734)
(325.0, 389.71143170299734)
(375.0, 389.71143170299734)
(425.0, 389.71143170299734)
(475.0, 389.71143170299734)
(525.0, 389.71143170299734)
(575.0, 389.71143170299734)
(625.0, 389.71143170299734)
(675.0, 389.71143170299734)
(725.0, 389.71143170299734)
(300.0, 433.0127018922193)
(350.0, 433.0127018922193)
(400.0, 433.0127018922193)
(450.0, 433.0127018922193)
(500.0, 433.0127018922193)
(550.0, 433.0127018922193)
(600.0, 433.0127018922193)
(650.0, 433.0127018922193)
(700.0, 433.0127018922193)
(750.0, 433.0127018922193)There exists blockade interactions between hard-core particles. We connect two lattice sites within blockade radius by an edge. Two ends of an edge can not both be occupied by particles.
blockade_radius = 55
using GenericTensorNetworks: show_graph, unit_disk_graph
using GenericTensorNetworks.Graphs: edges, nv
graph = unit_disk_graph(vec(sites), blockade_radius)
show_graph(graph, sites; texts=fill("", length(sites)))
These constraints defines an independent set problem that characterized by the following energy based model. Let $G = (V, E)$ be a graph, where $V$ is the set of vertices and $E$ is the set of edges. The energy model for the hard-core lattice gas problem is
\[E(\mathbf{n}) = -\sum_{i \in V}w_i n_i + U \sum_{(i, j) \in E} n_i n_j\]
where $n_i \in \{0, 1\}$ is the number of particles at site $i$, and $w_i$ is the weight associated with it. For unweighted graphs, the weights are uniform. $U$ is the repulsive interaction strength between two particles. To represent the independence constraint, we let $U = \infty$, i.e. coexitence of two particles at two sites connected by an edge is completely forbidden. The solution space hard-core lattice gas is equivalent to that of an independent set problem. The independent set problem involves finding a set of vertices in a graph such that no two vertices in the set are adjacent (i.e., there is no edge connecting them). One can create a tensor network based modeling of an independent set problem with package GenericTensorNetworks.jl.
using GenericTensorNetworks
problem = IndependentSet(graph)GenericTensorNetworks.IndependentSet{GenericTensorNetworks.UnitWeight}(Graphs.SimpleGraphs.SimpleGraph{Int64}(261, [[2, 11], [1, 3, 11, 12], [2, 4, 12, 13], [3, 5, 13, 14], [4, 6, 14, 15], [5, 7, 15, 16], [6, 8, 16, 17], [7, 9, 17, 18], [8, 10, 18, 19], [9, 19, 20], [1, 2, 12, 21], [2, 3, 11, 13, 21, 22], [3, 4, 12, 14, 22, 23], [4, 5, 13, 15, 23, 24], [5, 6, 14, 16, 24, 25], [6, 7, 15, 17, 25, 26], [7, 8, 16, 18, 26, 27], [8, 9, 17, 19, 27, 28], [9, 10, 18, 20, 28, 29], [10, 19, 29, 30], [11, 12, 22, 31], [12, 13, 21, 23, 31, 32], [13, 14, 22, 24, 32, 33], [14, 15, 23, 25, 33, 34], [15, 16, 24, 26, 34, 35], [16, 17, 25, 27, 35, 36], [17, 18, 26, 28, 36, 37], [18, 19, 27, 29, 37, 38], [19, 20, 28, 30, 38, 39], [20, 29, 39, 40], [21, 22, 32, 41], [22, 23, 31, 33, 41, 42], [23, 24, 32, 34, 42, 43], [24, 25, 33, 35, 43, 44], [25, 26, 34, 36, 44, 45], [26, 27, 35, 37, 45, 46], [27, 28, 36, 38, 46, 47], [28, 29, 37, 39, 47, 48], [29, 30, 38, 40, 48, 49], [30, 39, 49, 50], [31, 32, 42, 51], [32, 33, 41, 43, 51, 52], [33, 34, 42, 44, 52, 53], [34, 35, 43, 45, 53, 54], [35, 36, 44, 46, 54, 55], [36, 37, 45, 47, 55, 56], [37, 38, 46, 48, 56, 57], [38, 39, 47, 49, 57, 58], [39, 40, 48, 50, 58, 59], [40, 49, 59, 60], [41, 42, 52, 61], [42, 43, 51, 53, 61, 62], [43, 44, 52, 54, 62, 63], [44, 45, 53, 55, 63, 64], [45, 46, 54, 56, 64, 65], [46, 47, 55, 57, 65, 66], [47, 48, 56, 58, 66, 67], [48, 49, 57, 59, 67, 68], [49, 50, 58, 60, 68, 69], [50, 59, 69, 70], [51, 52, 62, 71], [52, 53, 61, 63, 71, 72], [53, 54, 62, 64, 72, 73], [54, 55, 63, 65, 73, 74], [55, 56, 64, 66, 74, 75], [56, 57, 65, 67, 75, 76], [57, 58, 66, 68, 76, 77], [58, 59, 67, 69, 77, 78], [59, 60, 68, 70, 78, 79], [60, 69, 79, 80], [61, 62, 72, 81], [62, 63, 71, 73, 81, 82], [63, 64, 72, 74, 82, 83], [64, 65, 73, 75, 83, 84], [65, 66, 74, 76, 84, 85], [66, 67, 75, 77, 85, 86], [67, 68, 76, 78, 86, 87], [68, 69, 77, 79, 87, 88], [69, 70, 78, 80, 88, 89], [70, 79, 89, 90], [71, 72, 82, 91], [72, 73, 81, 83, 91, 92], [73, 74, 82, 84, 92, 93], [74, 75, 83, 85, 93, 94], [75, 76, 84, 86, 94, 95], [76, 77, 85, 87, 95, 96], [77, 78, 86, 88, 96, 97], [78, 79, 87, 89, 97, 98], [79, 80, 88, 90, 98, 99], [80, 89, 99, 100], [81, 82, 92], [82, 83, 91, 93], [83, 84, 92, 94], [84, 85, 93, 95], [85, 86, 94, 96], [86, 87, 95, 97], [87, 88, 96, 98], [88, 89, 97, 99], [89, 90, 98, 100], [90, 99]]), GenericTensorNetworks.UnitWeight())There are plenty of discussions related to solution space properties in the GenericTensorNetworks documentaion page. In this example, we show how to use TensorInference to use probabilistic inference for understand the finite temperature properties of this statistical model. We use TensorNetworkModel to convert a combinatorial optimization problem to a probabilistic model. Here, we let the inverse temperature be $\beta = 3$.
Probabilistic modeling correlation functions
using TensorInference
β = 3.0
pmodel = TensorNetworkModel(problem, β)TensorNetworkModel{Int64, OMEinsum.DynamicNestedEinsum{Int64}, Array{Float64}}
variables: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100
contraction time = 2^17.727, space = 2^12.0, read-write = 2^15.72The partition function of this statistical model can be computed with the probability function.
partition_func = probability(pmodel)exp(107.422678904294) * fill(1.0)The default return value is a log-rescaled tensor. Use indexing to get the real value.
partition_func[]4.4985927541460835e46The marginal probabilities can be computed with the marginals function, which measures how likely a site is occupied.
mars = marginals(pmodel)
show_graph(graph, sites; vertex_colors=[(b = mars[[i]][2]; (1-b, 1-b, 1-b)) for i in 1:nv(graph)], texts=fill("", nv(graph)))
The can see the sites at the corner is more likely to be occupied. To obtain two-site correlations, one can set the variables to query marginal probabilities manually.
pmodel2 = TensorNetworkModel(problem, β; mars=[[e.src, e.dst] for e in edges(graph)])
mars = marginals(pmodel2);We show the probability that both sites on an edge are not occupied
show_graph(graph, sites; edge_colors=[(b = mars[[e.src, e.dst]][1, 1]; (1-b, 1-b, 1-b)) for e in edges(graph)], texts=fill("", nv(graph)), config=GraphDisplayConfig(; edge_line_width=5))
The most likely configuration
The MAP and MMAP can be used to get the most likely configuration given an evidence. The relavant function is most_probable_config. If we fix the vertex configuration at one corner to be one, we get the most probably configuration as bellow.
pmodel3 = TensorNetworkModel(problem, β; evidence=Dict(1=>1))
mars = marginals(pmodel3)
logp, config = most_probable_config(pmodel3)(102.0, [1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1])The log probability is 102. Let us visualize the configuration.
show_graph(graph, sites; vertex_colors=[(1-b, 1-b, 1-b) for b in config], texts=fill("", nv(graph)))
The number of particles is
sum(config)34Otherwise, we will get a suboptimal configuration.
pmodel3 = TensorNetworkModel(problem, β; evidence=Dict(1=>0))
logp2, config2 = most_probable_config(pmodel)(102.0, [1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1])The log probability is 99, which is much smaller.
show_graph(graph, sites; vertex_colors=[(1-b, 1-b, 1-b) for b in config2], texts=fill("", nv(graph)))
The number of particles is
sum(config2)
# Sampling configurations34One can ue sample to generate samples from hard-core lattice gas at finite temperature. The return value is a matrix, with the columns correspond to different samples.
configs = sample(pmodel3, 1000)
sizes = sum.(configs)
[count(==(i), sizes) for i=0:34] # counting sizes35-element Vector{Int64}:
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
2
4
22
60
128
185
200
200
126
60
13
0This page was generated using Literate.jl.