Update (July 9 2025)

 

Update 7/9/25: By The End Of Summer 2025

Main Updates: Looking In The Near Future

(1)Detailed Ising model and spin glass system replications

(2) Noisy Coin

(3)Original Research, New Simulation technique on spin glass systems with (2)

I want to get to (3) before my birthday on the 13th although that's unlikely :)

More Updates:

I've got a social media account dedicated to this page: https://bsky.app/profile/gaurav-kamath.bsky.social

I've abandoned by personal codex approach due to it's lack of replication and small scale experimentation 

needed to approach new ideas. It has however helped me generate new ideas that I wish to build in the

near future.

I also apologize for any indentations problems that may arise in my released code as copying it from

 Jupyter Notebook into blogspot was done a little poorly.

Replicating Sandpile Research: Nature Fig 1a) From Greedy Control Of Cascading Failures In Interdependent Networks

Mu Controlled Sandpile Model 

#I switched to python from C to make the sandpile model to make my life easier

#replicating Nature fig 1a: https://www.nature.com/articles/s41598-021-82843-8/figures/1

import numpy as np

import matplotlib.pyplot as plt

import random

def make_area(x, y):

    return [[0]*x for _ in range(y)]

def find_at_capacity_sites(height, XMAX, YMAX):

    return [(x,y) for x in range(XMAX) for y in range(YMAX) if height[y][x] >= 3] 

def find_not_at_capacity_sites(height, XMAX, YMAX):

    return [(x,y) for x in range(XMAX) for y in range(YMAX) if height[y][x] < 3] 

def flip(mu):

    return random.random() < mu

def add_grain(area, mu, XMAX, YMAX):

    at_capacity = find_at_capacity_sites(area, XMAX, YMAX)

    not_at_capacity = find_not_at_capacity_sites(area, XMAX, YMAX)

    if not at_capacity:

        rx, ry = random.choice(not_at_capacity)

    else:

        if flip(mu):

            rx, ry = random.choice(at_capacity)

        else:

            rx, ry = random.choice(not_at_capacity)

    area[ry][rx] += 1

    return (rx, ry)

def topple(area, XMAX, YMAX): 

    avalanche_size = 0

    unstable_sites = []

    # Find initial unstable sites

    for y in range(YMAX):

        for x in range(XMAX):

            if area[y][x] >= 4:

                unstable_sites.append((x, y))

    

    while unstable_sites:

        x, y = unstable_sites.pop()

        if area[y][x] < 4:

            continue

        avalanche_size += 1

        area[y][x] -= 4

        # Update neighbors

        for dx, dy in [(1,0), (-1,0), (0,1), (0,-1)]:

            nx, ny = x + dx, y + dy

            if 0 <= nx < XMAX and 0 <= ny < YMAX:

                area[ny][nx] += 1

                if area[ny][nx] >= 4:  # Check if neighbor becomes unstable

                    unstable_sites.append((nx, ny))

    return avalanche_size

def run_sandpile(mu, grains_to_add, XMAX, YMAX):

    area = make_area(XMAX, YMAX)

    avalanche_sizes = []

    for _ in range(grains_to_add):

        add_grain(area, mu, XMAX, YMAX)

        aval_size = topple(area, XMAX, YMAX)

        if aval_size > 0:

            avalanche_sizes.append(aval_size)

    return avalanche_sizes

def plot_loglog_cascade_size(avalanche_data, mu_values):

    plt.figure(figsize=(8,6))

    

    # Use markers/colors from paper

    markers = {0.20: 'o', 0.37: 's', 0.90: '^'}  # circle, square, triangle

    colors = {0.20: 'blue', 0.37: 'green', 0.90: 'red'}

    

    for mu, avalanches in zip(mu_values, avalanche_data):

        if len(avalanches) == 0:

            continue

            

        max_s = max(avalanches)

        bins = np.logspace(0, np.log10(max_s), 30)

        hist, bin_edges = np.histogram(avalanches, bins=bins, density=True)

        bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

        

        # Plot only non-zero points

        non_zero = hist > 0

        plt.loglog(bin_centers[non_zero], hist[non_zero], 

                   marker=markers[mu], linestyle='-', lw=1.5, ms=6,

                   color=colors[mu], label=f'μ = {mu}')

    

    plt.xlabel("Cascade size, s", fontsize=14)

    plt.ylabel("Probability, P(s)", fontsize=14)

    plt.title("Avalanche Size Distribution", fontsize=16)

    plt.legend(fontsize=12)

    plt.grid(True, which="both", ls="--", alpha=0.7)

    plt.xlim(1, 1e5)  # Match paper's axis limits

    plt.ylim(1e-6, 1)

    plt.tight_layout()

    plt.show()

if __name__ == '__main__':

    XMAX, YMAX = 100, 100

    grains_to_add = 50000

    mu_values = [0.20, 0.37, 0.90]

    avalanche_data = []

    for mu in mu_values:

        print(f"Running sandpile for μ = {mu}...")

        avalanches = run_sandpile(mu, grains_to_add, XMAX, YMAX)

        avalanche_data.append(avalanches)

    

    plot_loglog_cascade_size(avalanche_data, mu_values)

Frenkel Kontorova Model On A 2-D Lattice

 # Frenkel-Kontorova Model

                                   

#Please keep in mind to run the model multiple times. Sometimes the restoration force or time steps

#aren't enough to restore the lattice back into the lattice back.

import random

import matplotlib.pyplot as plt

XMAX = 15

YMAX = 15

delta = 0.01  

dt = 0.1

tmax = 1500

def actualx(x): return x % XMAX

def actualy(y): return y % YMAX

def metric(x, x2): return min(abs(x - x2), XMAX - abs(x - x2))

class Particle:

    def __init__(self, x, y, vx, vy, q, m):

        self.x = x

        self.y = y

        self.vx = vx

        self.vy = vy

        self.q = q

        self.m = m

    def potential(self, x, y):

        dx = metric(x, self.x)

        dy = metric(y, self.y)

        r = (dx**2 + dy**2)**0.5

        return self.q / r if r != 0 else 0

class Defect:

    def __init__(self, x, y, strength):

        self.x = x

        self.y = y

        self.strength = strength

    def force_on(self, px, py):

        dx = metric(px, self.x)

        dy = metric(py, self.y)

        r2 = dx**2 + dy**2

        if r2 == 0: return 0, 0

        return (-0.5*self.strength * dx / r2, -0.5*self.strength * dy / r2)

def chargelistgenerator(m, q, xmax, ymax, vxmax, vymax, numberofparticles):

    return [Particle(random.uniform(0, xmax), random.uniform(0, ymax),

                     random.uniform(-vxmax, vxmax), random.uniform(-vymax, vymax), q, m)

            for _ in range(numberofparticles)]

def generate_defects(num_defects, strength):

    return [Defect(random.uniform(0, XMAX), random.uniform(0, YMAX), strength)

            for _ in range(num_defects)]

def gridPotentialAtPoint(x, y, chargelist):

    return sum(charge.potential(x, y) for charge in chargelist)

def gridForceAtPointX(x, y, q, chargelist, Edrive):

    return q*Edrive - (q * (gridPotentialAtPoint(actualx(x + delta), actualy(y), chargelist) -

                            gridPotentialAtPoint(actualx(x - delta), actualy(y), chargelist)) / (2 * delta))

def gridForceAtPointY(x, y, q, chargelist):

    return - (q * (gridPotentialAtPoint(actualx(x), actualy(y + delta), chargelist) -

                   gridPotentialAtPoint(actualx(x), actualy(y - delta), chargelist)) / (2 * delta))

def particleVelocityUpdate(particle, chargelist, defects, Edrive):

    fx = gridForceAtPointX(particle.x, particle.y, particle.q, chargelist, Edrive)

    fy = gridForceAtPointY(particle.x, particle.y, particle.q, chargelist)

#account for force from defects 

    for defect in defects:

        dfx, dfy = defect.force_on(particle.x, particle.y)

        fx += dfx

        fy += dfy

    particle.vx += fx / particle.m

    particle.vy += fy / particle.m

    particle.vx *= 0.999

    particle.vy *= 0.999

def particlePositionUpdate(particle, deltat):

    particle.x = actualx(particle.x + particle.vx * deltat)

    particle.y = actualy(particle.y + particle.vy * deltat)

def run_simulation(Edrive):

    particles = chargelistgenerator(m=1, q=0.5, xmax=XMAX, ymax=YMAX,

                                     vxmax=0.05, vymax=0.05, numberofparticles=20)

    defects = generate_defects(num_defects=11, strength=0.01)

    t = 0

    while t < tmax:

        for p in particles:

            particleVelocityUpdate(p, particles, defects, Edrive)

        for p in particles:

            particlePositionUpdate(p, dt)

        t += dt

    return particles, defects

particles_before, defects = run_simulation(Edrive=0)

particles_after, _ = run_simulation(Edrive=450)

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

for i, (particles, title) in enumerate(zip([particles_before, particles_after],

                                           ['Before Edrive', 'After Edrive'])):

    x_vals = [p.x for p in particles]

    y_vals = [p.y for p in particles]

    dx_vals = [d.x for d in defects]

    dy_vals = [d.y for d in defects]

    axs[i].scatter(x_vals, y_vals, c='blue', label='Particles')

    axs[i].scatter(dx_vals, dy_vals, c='red', marker='x', s=100, label='Defects')

    axs[i].set_xlim(0, XMAX)

    axs[i].set_ylim(0, YMAX)

    axs[i].set_title(title)

    axs[i].set_aspect('equal')

    axs[i].grid(True)

plt.suptitle("Frenkel-Kontorova: Before and After Applying Driving Force", fontsize=16)

plt.show()

Lorenz Driven Operationally Deterministic Ising Model

#Lorenz Driven Operationally Deterministic Ising Model

#Please cite me for credit if you are planning to use my code

#Operation:

# Only the initial state is randomized. all else (core functions) operates deterministically. 

# Here we deviate from the traditional Ising Model which are probabilistic, stochastic in operation: 

#X(t),Y(t),Z(t) are coupled through a system of continuous non-linear differential equations.

# raw Y(t) -> used to select sites.

# normalized X(t) -> used to define spin flip conditions. 

# raw Z(t) + forcing -> equation where we add our forcing component (forcing = 0 if you want pattern formation to occur)

# The free parameters of the Lorenz equation(r,s,b) are material specific and temperature independent properties of the material 

# which I aptly call "spin-sensivity parameters" of a material on which the experiment is being conducted. 

# Model Motivation & Model Outcome:

#The goal here is to bridge the gap between theory and practice of solid state physics by abandoning random.random() 

#for the internal operations of the Ising model. 

#The Ising Model is an ontologically deterministic spin system, but we are limited in

#knowing the initial conditions due to finite precision from our measurement instruments. Different potential patterns that emerge vary. 

#depending on perturbing our inputs which we simulate due to said limitations in the knowledge of the initial condition of

 #the system.These are extrinsic limitations. In addition the model has built structural limitations that it acknowledges:

 #computers have finite resolution and time steps despite the Lorenz model being innately continuous rather than discrete

 #.So our model reflects this intrinsic limitation and should provide fruitful study for how extrinsic and intrinsic

 #limitations interact with one another in limiting predictability of the exact patterns formed in experiment.

#the coupled equations in real life. The coupled nature of the equations in addition to this intrinsic limitation with delta #time and fidelity could also give the model more utility to the modeler over simplistic substitution of random.random with #one dimensional discrete chaotic map where time steps are discrete  f: x(t) -> x(t+1)  whereas my model can do outputs for t+0.001.

# The consequence being that different pattern formation appears "stochastic" due to limitations in knowledge of the initial conditions.

#At the same time we have a "know how" that some kind of pattern of spin 'up' and 'down' will eventually emerge due to #dropping temperature. 

#when you add a forcing onto Z you are altering the geometry of the attractor and thus altering the effect of the spin-#sensitivity parameters on the spin of the material on it. You consider it akin to chemical or physical alteration of a #material. This means we may be able to study the Ising with time dependent forcing with reversible alterations of the  material 

#Potential Numerical Advantages:

# pre-packaged inputs would offer and speed and predictability in computation over random.random() operations.

import numpy as np

import matplotlib.pyplot as plt

from decimal import Decimal, getcontext, ROUND_HALF_UP

def lorenz_step(state, dt=0.01, s=15, r=23, b=2):

    x, y, z = state

    x_dot = s * (y - x)

    y_dot = r * x - y - x * z

    z_dot = x * y - b * z

    new_state = np.array([x + x_dot * dt, y + y_dot * dt, z + z_dot * dt])

    return new_state

class LorenzRNG:

    def __init__(self, initial_state=(0.0, 1.0, 1.05), dt=0.01, buffer_size=1000, skip_transient=500):

        self.state = np.array(initial_state)

        self.dt = dt

        self.x_buffer = []

        self.transient_counter = 0

        self.skip_transient = skip_transient

        self.buffer_size = buffer_size

    def next(self):

        self.state = lorenz_step(self.state, self.dt)

        x, y, _ = self.state

        if self.transient_counter < self.skip_transient:

            self.transient_counter += 1

            return 0.5, y  

        self.x_buffer.append(x)

        if len(self.x_buffer) > self.buffer_size:

            self.x_buffer.pop(0)

            

        x_min = min(self.x_buffer)

        x_max = max(self.x_buffer)

        if x_max != x_min:

            x_norm = (x - x_min) / (x_max - x_min)

        else:

            x_norm = 0.5 

        x_norm = np.clip(x_norm, 0, 1)

        return x_norm, y

XMAX, YMAX = 50, 50

N = XMAX * YMAX

def delta_energy(x, y, grid):

    return grid[x, y] * (

        grid[(x + 1) % XMAX, y] + grid[(x - 1) % XMAX, y] +

        grid[x, (y + 1) % YMAX] + grid[x, (y - 1) % YMAX]

    )

def switch_spin(x, y, grid):

    grid[x, y] *= -1

def grid_with_finder(rows, columns, fidelity, target_num):

    getcontext().prec = 28

    getcontext().rounding = ROUND_HALF_UP

    fidelity_decimal = Decimal(str(fidelity))

    target_decimal = Decimal(str(target_num))

    grid_max = rows * columns * fidelity_decimal

    if grid_max != 0:

        adjusted_target = (target_decimal % grid_max) if target_decimal >= 0 else (grid_max - (abs(target_decimal) % grid_max))

    else:

        adjusted_target = Decimal('0')

    grid = []

    current_value = Decimal('0')

    

    for y in range(rows):

        row = []

        for x in range(columns):

            wrapped_value = current_value % grid_max if grid_max != 0 else current_value

            quantized_value = wrapped_value.quantize(fidelity_decimal)

            row.append(quantized_value)

            current_value += fidelity_decimal

        

        if y % 2 == 1:

            row.reverse()

        grid.insert(0, row)

    closest = None

    min_diff = float('inf')

    closest_x, closest_y = -1, -1

    for y in range(rows):

        for x in range(columns):

            diff = abs(grid[y][x] - adjusted_target)

            if diff < min_diff:

                min_diff = diff

                closest = grid[y][x]

                closest_x, closest_y = x, y

                

    return grid, closest, (closest_x, rows - 1 - closest_y), adjusted_target

def metropolis_step(grid, temperature, lorenz_rng, scaling_factor=0.153):

    x_norm, raw_y = lorenz_rng.next()

    _, _, (x, y), _ = grid_with_finder(XMAX, YMAX, scaling_factor, raw_y)

    dE = delta_energy(x, y, grid)

    if dE <= 0:

        switch_spin(x, y, grid)

    else:

        if x_norm < np.exp(-2 * dE / temperature):

            switch_spin(x, y, grid)

def temperature_schedule(step, total_steps, T_start, T_end):

    return T_end + (T_start - T_end) * np.exp(-5 * step / total_steps)

def run_simulation_with_snapshots(T_start, T_end, snapshot_steps, lorenz_initial_state, scaling_factor=0.153, initial_grid=None):

    lorenz_rng = LorenzRNG(initial_state=lorenz_initial_state, dt=0.01)

    if initial_grid is None:

        grid = np.random.choice([-1, 1], size=(XMAX, YMAX))

    else:

        grid = np.copy(initial_grid)

    snapshots = {}

    total_steps = max(snapshot_steps)

    for step in range(1, total_steps + 1):

        temperature = temperature_schedule(step, total_steps, T_start, T_end)

        for _ in range(XMAX * YMAX):

            metropolis_step(grid, temperature, lorenz_rng, scaling_factor=scaling_factor)

        if step in snapshot_steps:

            snapshots[step] = np.copy(grid)

    return snapshots

np.random.seed(42)

initial_grid = np.random.choice([-1, 1], size=(XMAX, YMAX))

fixed_initial_condition = (0.15, 2.32, 6.04)

scaling_factors = [0.0001]  

snapshot_steps = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

T_start = 100

T_end = 0.1

fig, axes = plt.subplots(len(scaling_factors), len(snapshot_steps), figsize=(16, 3 * len(scaling_factors)))

for i, scale in enumerate(scaling_factors):

    snapshots = run_simulation_with_snapshots(

        T_start, T_end, snapshot_steps, fixed_initial_condition,

        scaling_factor=scale,

        initial_grid=initial_grid

    )

    for j, step in enumerate(snapshot_steps):

        ax = axes[i, j] if len(scaling_factors) > 1 else axes[j]

        ax.imshow(snapshots[step], cmap='gray', interpolation='nearest')

        ax.set_title(f"Step {step}")

        ax.axis('off')

plt.suptitle(f"Lorenz Driven Operationally Deterministic Ising Model. Scale: {scale}", fontsize=18)

plt.tight_layout(rect=[0, 0, 1, 0.96])

plt.show()

Internship Documentation & My Personal Codex

 

The internship's subject matter is solid state physics & complexity sciences adjacent.

 

(1)Internship:https://docs.google.com/document/d/1g1H2tHqhZvM84_8PzWZJH5ote5eglgGKFTA6D07hVDg/edit#heading=h.vga3z4m33bei

 (2) Personal Codex: https://docs.google.com/document/d/1ZagN2UDr82n5hrnfReOri63jSu4GErI3V94h_qH5nkI/edit?usp=sharing

The documents should be accessible to anyone from UCDAVIS or by request.

                                                                                                    

Sandpile On Different Topologies

How would one program a sandpile on different topologies? Well I've thought about it and drawn it out:

Link: https://docs.google.com/document/d/1bsV9MWaCSkGUy4q77-S1i6x4JXGbyEoisEIlnS015xY/edit

I was thinking of how sandpiles would look on different topologies. For a double torus and so on, you would just do something similar for the single torus except you would connect like boundaries as shown bellow:

Anyways this post should make it easier to adjust the boundary conditions of my prototype program to the desired topologies.

Sandpile Model Prototype

As promised; I recently whipped up a Non-Abelian Sandpile model written in C. It's a fairly simple code where there is a grid of cells and if the sand grains in a cell is greater than the threshold, then the sand in that cell topples to the neighbors. In an abelian model this sand is distributed equally to its surrounding neighbors. In the model I have chosen: this is distributed randomly to neighboring cells. There's still much to be done. For one I need to test the power law for this code. And play around to see what kind of patterns need to be generated. This will be done in Part 2 if ever.     

#include <stdio.h>

#include <math.h>

#include <stdbool.h>

#include <stdlib.h>

const int size1 = 50;

const int size2 = 50;

int amount = 300;

int xi = 25;

int yi = 25;

int threshold = 4;

int timeFinal = 500000;

int grid[size1][size2];

void generategrid(){

int i = 0, j = 0;

for(int i = 0; i < size1; i++){

for(int j = 0; j < size2; j++){

grid[i][j] = 0;

}

}

}

void printmatrix(){

int i = 0,j = 0;

for(int i = 0; i < size1; i++){

for(int j = 0; j < size2; j++){

printf(" %4d ", grid[i][j]);

}

printf("\n");

}

}

void Topple(int i, int j) {

int total;

if( (i < 0 || i >= size1) || (j < 0 && j >= size2) ) {// check for invalid grid coorinates

return;

}

if(grid[i][j] < threshold){ // do we really need to topple?

return;

}

total = threshold;

if(i==0 || j==0 || i==size1-1 || j==size2-1){

//on borders, sand may drop off the edges/corners

if(i==0 && j==0 ){ //bottom left corner

int r1 = rand()%(total);

grid[0][1] += r1;

total -= r1;

grid[1][0] += rand()%total;

}

else if(i==size1-1 && j==size2-1 ){ //top right corner

int r1 = rand()%(total);

grid[size1-2][size2-1] += r1;

total -= r1;

grid[size1-1][size2-2] +=rand()% total;

}

else if(i==size1-1 && j==0 ){ // bottom right corner

int r1 = rand()%(total);

grid[size1-2][0] += r1;

total -= r1;

grid[size1-1][1] += rand()%total;

}

else if(i==0 && j==size2-1){ // top left corner

int r1 =rand()%(total);

grid[0][size2-2] +=r1;

total -= r1;

grid[1][size2-1] += rand()%total;

}

else if(i==0 && (j!=0 || j!=size2-1)){ //left edge but not corners

int r1 = rand()%(total);

grid[i][j+1] += r1;

total -= r1;

int r2 = rand()%(total);

grid[i][j-1] += r2;

total -= r2;

grid[i+1][j] += rand()%total;

}

else if(j==0 && (i!=0 || i!=size1-1)){ //bottom edge but not corners

int r1 = rand()%(total);

grid[i+1][j] += r1;

total -=r1;

int r2 = rand()%(total);

grid[i-1][j] += r2;

total -= r2;

grid[i][j+1] += rand()%total;

}

else if(i==size1-1 && (j!=0 || j!=size2-1)){ //right edge but not corners

int r1 = rand()%(total);

grid[i][j+1] += r1;

total -= r1;

int r2 = rand()%(total);

grid[i][j-1] += r2;

total -= r2;

grid[i-1][j] += rand()%total;

}

else if(j== size2-1 && (i!=0 || i!=size2-1)){ //top edge but not corners

int r1 = rand()%(total);

grid[i+1][j] += r1;

total -= r1;

int r2 = rand()%(total);

grid[i-1][j] += r2;

total -= r2;

grid[i][j-1] += rand()%total;

}

}

else {

int r1 = rand()%total;

grid[i+1][j] += r1;

total -= r1;

int r2 = rand()%total;

grid[i-1][j] += r2;

total -= r2;

int r3 = rand()%total;

grid[i][j+1] += r3;

total -= r3;

grid[i][j-1] += total; // no sand lost in this case

}

grid[i][j] = 0;

//debuggers:

//printmatrix();

// printf("Press any key to continue\n");

//getchar();

return;

}

void CheckThreshold(){

int i, j = 0;

    int gridfinal;

for( i = 0; i < size1; i++){

for( j = 0; j < size2; j++){

if(grid[i][j]>=threshold){

                gridfinal = grid[i][j]-threshold;

do{

Topple(i,j);

Topple(i+1,j);

Topple(i,j+1);

Topple(i,j-1);

}while(grid[i][j]>gridfinal);

}

}

}

}

void addgrain(){

int count;

int t;

do{

grid[xi][yi] += amount;

CheckThreshold();

t += 1;

}while(t < timeFinal);

}

int main(){

generategrid();

addgrain();

printmatrix();

}

Update July 13: I'll test these when I get the chance. But I plan on improving the random addition of sand grains as well. Here are the updates to the code bellow. An important limitation to the new add grain feature is that is assumes uniform distribution.  

                             int rx = 4; // new var for boundary of sand addition

                int ry = 4; // new var for boundary of sand addition

grid[xi+rand()%(rx)][yi+rand()%(ry)] += amount; //mod for add grain

Update Dec 3. I realized there was an issue with the way I defined the toppling conditions. In alignment with sandpile I modded the code (but beware its not been tested yet).

I caught this when reading the rules for the basic abelian sandpile model with the assumed threshold of 4 grains of sand.

Any site ${\displaystyle (x,y)}$ with

${\displaystyle z(x,y)\ge 4}$

is unstable and can topple (or fire), sending one of its chips to each of its four neighbors:

${\displaystyle z(x,y)\to z(x,y)-4,}$

${\displaystyle z(x\pm 1,y)\to z(x\pm 1,y)+1,}$

${\displaystyle z(x,y\pm 1)\to z(x,y\pm 1)+1.}$

Obviously 4 was the threshold for the tile so the while condition should operate until the grid value is gridfinal= grid[i][j]-threshold.

and then the threshold value is distributed to its neighbors as seen with the fact that the topple function's total (which is the max number of grains it can distribute) = to the global variable: threshold.