Numerical Simulation of Classical Ising Model using the Core Interface

Numerical Simulation of Classical Ising Model using the Core Interface#

In this tutorial, we will use OpenJij’s core interface (core python interface) to perform numerical simulations of the Ising model with random interactions and random longitudinal magnetic fields.

First, we define Graph and $J_{ij}, h_i$ for the system you wish to simulate numerically.

import openjij.cxxjij.graph as G
# Set the problem size to 100
N = 100

graph = G.Dense(N)
# Use below for sparse
#graph = G.Sparse(N)

Next, we set $J_{ij}, h_i$ . We set the values generated from a Gaussian distribution with a mean of 0 and a standard deviation of 1.

# Use numpy for random number generation
# !pip install numpy

import numpy as np
mu, sigma = 0, 1

for i in range(N):
    for j in range(N):
        # Jij value would be too large, so we use 1/N for standardization.
        graph[i,j] = 0 if i == j else np.random.normal()/N

for i in range(N):
    graph[i] = np.random.normal()/N

The longitudinal magnetic field can be accessed either by graph[i] or by graph[i,i]. $J_{ij}$ and $J_{ji}$ are automatically the same value by definition of the Ising model. Let us try the following output.

graph[20] = 0.5
print(graph[20,20])
print(graph[20])
graph[12,34] = -0.6
print(graph[12,34])
print(graph[34,12])

0.5
0.5
-0.6
-0.6

System#

Next, define the system to perform the calculations.

Here we want to perform a numerical simulation of the classical Ising model, so we create a system for that model with system.make_classical_ising.

import openjij.cxxjij.system as S

mysystem = S.make_classical_ising(graph.gen_spin(), graph)

Here, we assign a randomly generated spin to the first argument and the Graph itself to the second. This creates a system of classical Ising models whose initial spin configuration is graph.gen_spin().

We can also access the system directly and read the values.

print(mysystem.spin)

[-1.  1. -1.  1.  1.  1. -1. -1.  1. -1.  1.  1. -1. -1. -1.  1. -1.  1.
 -1. -1. -1. -1.  1.  1. -1. -1.  1. -1.  1. -1. -1.  1.  1. -1. -1. -1.
  1. -1. -1.  1. -1.  1.  1.  1. -1.  1. -1. -1.  1.  1.  1. -1. -1.  1.
  1.  1.  1.  1.  1.  1. -1.  1.  1. -1.  1. -1. -1. -1. -1.  1. -1. -1.
  1. -1.  1. -1.  1.  1.  1. -1. -1.  1. -1. -1. -1.  1.  1.  1. -1. -1.
  1.  1.  1. -1.  1.  1.  1.  1.  1. -1.  1.]

Run the Algorithm with Updater#

After defining the System, we select Updater and run the Algorithm.

Updater#

The Updater for the classical Ising model is mainly:

SingleSpinFlip (update one spin at a time using the Metropolis-Hasting method)
SwendsenWang (cluster updates by the SwendsenWang method)

To run the Algorithm, you will need a schedule list. Let us start by creating a schedule list.

Algorithm#

Schedule list#

The schedule list is given as a list of (parameters, number of Monte Carlo steps). The value to be entered in the parameter depends on the System. For example, for the classical Ising model, the parameter is the inverse temperature $\beta$ . Therefore, this refers to a list of temperature schedules. Let us set the following example.

schedule_list = [(0.01, 10),(10, 80),(0.1, 30)]

This means that 10 Monte Carlo steps are performed at inverse temperature $\beta=0.01$ , 80 steps at $\beta=10$ , and 30 steps at $\beta=0.1$ , for a total of 120 Monte Carlo steps.
In annealing, the inverse temperature is often increased proportionally, so you can easily create a schedule by using the make_classical_schedule_list in the utility. The temperature is changed in 10 steps from $\beta=0.1$ to $\beta=50$ , with 20 Monte Carlo steps calculated at each temperature. A total of 200 Monte Carlo steps are calculated.

import openjij.cxxjij.utility as U
schedule_list = U.make_classical_schedule_list(0.1, 50, 20, 10)
print(schedule_list)

[((beta: 0.100000) mcs: 20), ((beta: 0.199474) mcs: 20), ((beta: 0.397897) mcs: 20), ((beta: 0.793701) mcs: 20), ((beta: 1.583223) mcs: 20), ((beta: 3.158114) mcs: 20), ((beta: 6.299605) mcs: 20), ((beta: 12.566053) mcs: 20), ((beta: 25.065966) mcs: 20), ((beta: 50.000000) mcs: 20)]

Run the Algorithm#

By writing Algorithm_[Updater]_run, the calculation is performed with the specified Updater. In the following example, SingleSpinFlip is executed.

import openjij.cxxjij.algorithm as A
A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list)

The process is over instantly, but a total of 200 Monte Carlo steps have been computed during this time.

A.Algorithm_SingleSpinFlip_run(mysystem, seed, schedule_list) allows you to run the calculation with the seed fixed. This can be used to make the results reproducible.

We use the callback to get the system for each 1 Monte Carlo step during the execution of the Algorithm. For the classical Ising model, we create a function with the system and parameters (inverse temperature) as arguments.
As an example, below we create a callback that records the value of the energy of the system.

energies = []

def callback_log_energy(system, beta):
    #graph is the object defined previously in the Graph module
    energies.append(graph.calc_energy(system.spin))

The same Algorithm is executed using this callback.

#Take a longer schedule (total of 20000 Monte Carlo steps)
schedule_list = U.make_classical_schedule_list(0.1, 50, 200, 100)
A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list, callback_log_energy)

The recorded system energy is plotted with Monte Carlo steps on the horizontal axis and energy on the vertical axis as follows:

# !pip install matplotlib

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(range(len(energies)), energies)
plt.xlabel('Monte Carlo step')
plt.ylabel('energy')
plt.show()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[13], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
      2 import matplotlib.pyplot as plt
      3 plt.plot(range(len(energies)), energies)

File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/IPython/core/interactiveshell.py:2456, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2454     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2455 with self.builtin_trap:
-> 2456     result = fn(*args, **kwargs)
   2458 # The code below prevents the output from being displayed
   2459 # when using magics with decorator @output_can_be_silenced
   2460 # when the last Python token in the expression is a ';'.
   2461 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/IPython/core/magics/pylab.py:99, in PylabMagics.matplotlib(self, line)
     97     print("Available matplotlib backends: %s" % backends_list)
     98 else:
---> 99     gui, backend = self.shell.enable_matplotlib(args.gui.lower() if isinstance(args.gui, str) else args.gui)
    100     self._show_matplotlib_backend(args.gui, backend)

File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/IPython/core/interactiveshell.py:3630, in InteractiveShell.enable_matplotlib(self, gui)
   3609 def enable_matplotlib(self, gui=None):
   3610     """Enable interactive matplotlib and inline figure support.
   3611 
   3612     This takes the following steps:
   (...)
   3628         display figures inline.
   3629     """
-> 3630     from matplotlib_inline.backend_inline import configure_inline_support
   3632     from IPython.core import pylabtools as pt
   3633     gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)

File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib_inline/__init__.py:1
----> 1 from . import backend_inline, config  # noqa
      2 __version__ = "0.1.7"  # noqa

File /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages/matplotlib_inline/backend_inline.py:6
      1 """A matplotlib backend for publishing figures via display_data"""
      3 # Copyright (c) IPython Development Team.
      4 # Distributed under the terms of the BSD 3-Clause License.
----> 6 import matplotlib
      7 from matplotlib import colors
      8 from matplotlib.backends import backend_agg

ModuleNotFoundError: No module named 'matplotlib'

The energy is gradually decreasing as the annealing proceeds. Thus, the callback is useful to know how the system looks like while the Algorithm is running.

Result#

With result.get_solution we get the spin configuration that is the result of the calculation. For the classical Ising model, we can also directly refer to mysystem.spin to get the spin configuration. However, result.get_solution is a convenient method to do so for other systems as well.

import openjij.cxxjij.result as R
print(R.get_solution(mysystem))

[-1, -1, 1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, -1, 1, -1, 1, 1, -1, -1, -1, -1, -1, 1]

This spin configuration is the answer obtained by annealing. It is expected to be the ground state (close to the ground state) of the Hamiltonian.