U

#### Udyant Wig

answer by Ron Maimon*, and I saw that I could implement the simulation

described therein. I first did a version in Emacs Lisp, on reviewing

which I felt that I could develop a version in C, using the Emacs Lisp

version as a rough guide.

I quote from the README:

Briefly, the Ising model is a theoretical description of

ferromagnetism, which arises when the magnetic moments of atomic

spins align in the same direction. It is used to study phase

transitions and the emergence of equilibrium configurations.

I have not gone into the physics at all. I only model the changing

configurations of a lattice of bits.

A copy of the program is at Bitbucket at this URL:

https://bitbucket.org/udyant/simplified-square-lattice-ising-model/src

I have used the implementation of linked lists from kazlib written by

Kaz Kylheku, and found at this URL:

http://www.kylheku.com/~kaz/kazlib.html

The source code is used under the BSD license mentioned within them.

The license is reproduced below:

Copyright 1996-2012

Kaz Kylheku <[email protected]>

Vancouver, Canada

All rights reserved.

BSD License:

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions

are met:

1. Redistributions of source code must retain the above copyright

notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright

notice, this list of conditions and the following disclaimer in

the documentation and/or other materials provided with the

distribution.

3. The name of the author may not be used to endorse or promote

products derived from this software without specific prior

written permission.

THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR

IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED

WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.

I would like to request the readers to go over the source and suggest

any improvements to layout, style, structure, etc. I would be highly

appreciative and very grateful.

The amalgamated source is about 628 lines; I have spaced the modules by

two empty lines.

/* common.h begins */

#ifndef COMMON_H

#define COMMON_H

#include <stdio.h>

#include <stdlib.h>

#include <stdbool.h>

#include <ctype.h>

#include <string.h>

#include <errno.h>

extern double beta;

extern int dimension;

extern int scale;

extern int minority_size;

#endif

/* common.h ends */

/* common.c begins */

#include "common.h"

double beta = 0.5;

int dimension = 3;

int scale = 100;

int minority_size = 1;

/* common.c ends */

/* utilities.h begins */

#ifndef UTILITIES_H

#define UTILITIES_H

#include "common.h"

bool random_bit (void);

bool is_all_zeroes (char *string);

bool is_positive_integer (char *string);

int how_many (const char *s, int c);

int find (const char string [], char character, int offset);

char *copy_substring (char *substring, \

const char *string, \

size_t start, \

size_t end);

#endif

/* utilities.h ends */

/* utilities.c begins */

#include "utilities.h"

bool random_bit (void)

{

return rand () % 2;

}

bool is_all_zeroes (char *string)

{

char *sp;

for (sp = string; *sp != '\0'; sp++) {

if (*sp != '0') {

return false;

}

}

return true;

}

bool is_positive_integer (char *string)

{

char *sp;

for (sp = string; *sp != '\0'; sp++) {

if (!isdigit (*sp)) {

return false;

}

}

return (is_all_zeroes (string) ? false : true);

}

/* From /C: A Reference Manual/, 5th ed., by Harbison and Steele

* page 352

*/

int how_many (const char *s, int c)

{

int n = 0;

if (c == 0) return 0;

while (s) {

s = strchr (s, c);

if (s) n++, s++;

}

return n;

}

int find (const char string [], char character, int offset)

{

int length = strlen (string);

int i;

for (i = offset; i < length; i++) {

if (string

*== character)*

return i;

}

return -1;

}

char *copy_substring (char *substring, \

const char *string, \

size_t start, \

size_t end)

{

size_t length = end - start + 1;

errno = 0;

substring = malloc (1 + length);

if (substring == NULL) {

fprintf (stderr, "copy_substring: %s\n", strerror (errno));

exit (1);

}

strncpy (substring, string + start, length - 1);

substring [length] = '\0';

return substring;

}

/* utilities.c ends */

/* bits.h begins */

#ifndef BITS_H

#define BITS_H

#include "common.h"

struct bit {

int x;

int y;

int energy;

bool bitvalue;

};

typedef struct bit bit;

/* Bits */

bool bits_equal (bit one, bit two);

bool negate (bit *b);

bit complement (bit *b);

void flip (bit *b);

void dump_bit (bit b);

/* Lattices */

bit **allocate_lattice (bit **lattice);

void free_lattice (bit **lattice);

void print_lattice (bit **lattice);

bool lattices_equal (bit **a, bit **b);

#endif

/* bits.h ends */

/* bits.c begins */

#include "bits.h"

/* Bits */

bool bits_equal (bit one, bit two)

{

return one.bitvalue == two.bitvalue;

}

bool negate (bit *b)

{

return !b->bitvalue;

}

bit complement (bit *b)

{

bit complement;

complement.x = b->x;

complement.y = b->y;

complement.energy = b->energy;

complement.bitvalue = negate (b);

return complement;

}

void flip (bit *b)

{

b->bitvalue = negate (b);

}

void dump_bit (bit b)

{

printf ("x := %d\n", b.x);

printf ("y := %d\n", b.y);

printf ("energy := %d\n", b.energy);

printf ("bitvalue := %d\n", b.bitvalue);

}

/* Lattices */

bit **allocate_lattice (bit **lattice)

{

/* bit **lattice; */

int i;

errno = 0;

lattice = malloc (dimension * sizeof (bit *));

if (lattice == NULL)

{

fprintf (stderr, "%s\n", strerror (errno));

exit (1);

}

for (i = 0; i < dimension; i++)

{

lattice

return i;

}

return -1;

}

char *copy_substring (char *substring, \

const char *string, \

size_t start, \

size_t end)

{

size_t length = end - start + 1;

errno = 0;

substring = malloc (1 + length);

if (substring == NULL) {

fprintf (stderr, "copy_substring: %s\n", strerror (errno));

exit (1);

}

strncpy (substring, string + start, length - 1);

substring [length] = '\0';

return substring;

}

/* utilities.c ends */

/* bits.h begins */

#ifndef BITS_H

#define BITS_H

#include "common.h"

struct bit {

int x;

int y;

int energy;

bool bitvalue;

};

typedef struct bit bit;

/* Bits */

bool bits_equal (bit one, bit two);

bool negate (bit *b);

bit complement (bit *b);

void flip (bit *b);

void dump_bit (bit b);

/* Lattices */

bit **allocate_lattice (bit **lattice);

void free_lattice (bit **lattice);

void print_lattice (bit **lattice);

bool lattices_equal (bit **a, bit **b);

#endif

/* bits.h ends */

/* bits.c begins */

#include "bits.h"

/* Bits */

bool bits_equal (bit one, bit two)

{

return one.bitvalue == two.bitvalue;

}

bool negate (bit *b)

{

return !b->bitvalue;

}

bit complement (bit *b)

{

bit complement;

complement.x = b->x;

complement.y = b->y;

complement.energy = b->energy;

complement.bitvalue = negate (b);

return complement;

}

void flip (bit *b)

{

b->bitvalue = negate (b);

}

void dump_bit (bit b)

{

printf ("x := %d\n", b.x);

printf ("y := %d\n", b.y);

printf ("energy := %d\n", b.energy);

printf ("bitvalue := %d\n", b.bitvalue);

}

/* Lattices */

bit **allocate_lattice (bit **lattice)

{

/* bit **lattice; */

int i;

errno = 0;

lattice = malloc (dimension * sizeof (bit *));

if (lattice == NULL)

{

fprintf (stderr, "%s\n", strerror (errno));

exit (1);

}

for (i = 0; i < dimension; i++)

{

lattice

*= malloc (dimension * sizeof (bit));*

if (latticeif (lattice

*== NULL)*

{

fprintf (stderr, "%s\n", strerror (errno));

exit (1);

}

}

return lattice;

}

void free_lattice (bit **lattice)

{

int i;

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

free (lattice{

fprintf (stderr, "%s\n", strerror (errno));

exit (1);

}

}

return lattice;

}

void free_lattice (bit **lattice)

{

int i;

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

free (lattice

*);*

}

free (lattice);

}

void print_lattice (bit **lattice)

{

int i, j;

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

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

printf ("%d", (lattice}

free (lattice);

}

void print_lattice (bit **lattice)

{

int i, j;

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

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

printf ("%d", (lattice

*[j]).bitvalue);*

/*

dump_bit (lattice/*

dump_bit (lattice

*[j]);*

*/

}

putchar ('\n');

}

}

bool lattices_equal (bit **a, bit **b)

{

int i, j;

bool equal = true;

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

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

if (!bits_equal (a*/

}

putchar ('\n');

}

}

bool lattices_equal (bit **a, bit **b)

{

int i, j;

bool equal = true;

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

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

if (!bits_equal (a

*[j], b**[j]))*

equal = false;

}

}

return equal;

}

/* bits.c ends */

/* ising.h begins */

#ifndef ISING_H

#define ISING_H

#include <math.h>

#include <time.h>

#include "common.h"

#include "bits.h"

#include "utilities.h"

#include "list.h"

/* Core */

void check_and_collect_neighbor (list_t *neighborhood, int x, int y);

list_t *neighbors (list_t *neighborhood, bit b);

int energy (bit b);

void calculate_energies (void);

void ising_flip (bit *b);

void iterate (void);

void print_iteration (const char *format_string, int count);

/* Housekeeping */

static void dump_neighbors (list_t *neighbors);

static int count_ones (void);

static int count_zeroes (void);

/* Error checking */

void minority_error (void);

void beta_error (void);

void dimension_error (void);

#endif

/* ising.h ends */

/* ising.c begins */

#include "ising.h"

static bit **ising_lattice;

const char initial_format [] = "Initial configuration %6d";

const char iteration_format [] = "Iteration %6d";

const char final_format [] = "Final configuration %6d";

/* Core */

void check_and_collect_neighbor (list_t *neighborhood, int x, int y)

{

lnode_t *neighbor;

bit *payload;

if (x >= 0 && y >= 0 && x < dimension && y < dimension) {

payload = &ising_lattice [x] [y];

neighbor = lnode_create (payload);

if (neighbor == NULL) {

fprintf (stderr, "%s\n", strerror (errno));

lnode_destroy (neighbor);

free (payload);

exit (1);

}

list_append (neighborhood, neighbor);

}

}

list_t *neighbors (list_t *neighborhood, bit b)

{

neighborhood = list_create (-1);

check_and_collect_neighbor (neighborhood, b.x, b.y - 1);

check_and_collect_neighbor (neighborhood, b.x, b.y + 1);

check_and_collect_neighbor (neighborhood, b.x - 1, b.y);

check_and_collect_neighbor (neighborhood, b.x + 1, b.y);

return neighborhood;

}

int energy (bit b)

{

int e = 0;

lnode_t *neighbor;

bit *rover;

list_t *neighborhood = NULL;

neighborhood = neighbors (neighborhood, b);

for (neighbor = list_first (neighborhood); neighbor != NULL; \

neighbor = list_next (neighborhood, neighbor)) {

rover = lnode_get (neighbor);

if (b.bitvalue != rover->bitvalue) {

e++;

}

}

list_destroy_nodes (neighborhood);

list_destroy (neighborhood);

return e;

}

void calculate_energies (void)

{

int i, j;

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

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

(ising_latticeequal = false;

}

}

return equal;

}

/* bits.c ends */

/* ising.h begins */

#ifndef ISING_H

#define ISING_H

#include <math.h>

#include <time.h>

#include "common.h"

#include "bits.h"

#include "utilities.h"

#include "list.h"

/* Core */

void check_and_collect_neighbor (list_t *neighborhood, int x, int y);

list_t *neighbors (list_t *neighborhood, bit b);

int energy (bit b);

void calculate_energies (void);

void ising_flip (bit *b);

void iterate (void);

void print_iteration (const char *format_string, int count);

/* Housekeeping */

static void dump_neighbors (list_t *neighbors);

static int count_ones (void);

static int count_zeroes (void);

/* Error checking */

void minority_error (void);

void beta_error (void);

void dimension_error (void);

#endif

/* ising.h ends */

/* ising.c begins */

#include "ising.h"

static bit **ising_lattice;

const char initial_format [] = "Initial configuration %6d";

const char iteration_format [] = "Iteration %6d";

const char final_format [] = "Final configuration %6d";

/* Core */

void check_and_collect_neighbor (list_t *neighborhood, int x, int y)

{

lnode_t *neighbor;

bit *payload;

if (x >= 0 && y >= 0 && x < dimension && y < dimension) {

payload = &ising_lattice [x] [y];

neighbor = lnode_create (payload);

if (neighbor == NULL) {

fprintf (stderr, "%s\n", strerror (errno));

lnode_destroy (neighbor);

free (payload);

exit (1);

}

list_append (neighborhood, neighbor);

}

}

list_t *neighbors (list_t *neighborhood, bit b)

{

neighborhood = list_create (-1);

check_and_collect_neighbor (neighborhood, b.x, b.y - 1);

check_and_collect_neighbor (neighborhood, b.x, b.y + 1);

check_and_collect_neighbor (neighborhood, b.x - 1, b.y);

check_and_collect_neighbor (neighborhood, b.x + 1, b.y);

return neighborhood;

}

int energy (bit b)

{

int e = 0;

lnode_t *neighbor;

bit *rover;

list_t *neighborhood = NULL;

neighborhood = neighbors (neighborhood, b);

for (neighbor = list_first (neighborhood); neighbor != NULL; \

neighbor = list_next (neighborhood, neighbor)) {

rover = lnode_get (neighbor);

if (b.bitvalue != rover->bitvalue) {

e++;

}

}

list_destroy_nodes (neighborhood);

list_destroy (neighborhood);

return e;

}

void calculate_energies (void)

{

int i, j;

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

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

(ising_lattice

*[j]).energy = energy (ising_lattice**[j]);*

}

}

}

void ising_flip (bit *b)

{

int old_energy = b->energy;

int new_energy = energy (complement (b));

int delta_energy = new_energy - old_energy;

if (delta_energy > 0) {

double ising_threshold = exp (-(beta * delta_energy));

double random_value = rand () % scale / (double) scale;

if (ising_threshold < random_value) {

flip (b);

}

}

else {

flip (b);

}

}

void iterate (void)

{

int x = (int) rand () % dimension;

int y = (int) rand () % dimension;

ising_flip (&ising_lattice [x] [y]);

calculate_energies ();

}

void print_iteration (const char *format_string, int count)

{

size_t i;

size_t length = strlen (format_string);

putchar ('\n');

printf (format_string, count);

putchar ('\n');

for (i = 0; i < length + 3; i++)

putchar ('-');

putchar ('\n');

print_lattice (ising_lattice);

}

/* Housekeeping */

static void dump_neighbors (list_t *neighbors)

{

bit *rover;

lnode_t *neighbor;

for (neighbor = list_first (neighbors); neighbor != NULL; \

neighbor = list_next (neighbors, neighbor)) {

rover = lnode_get (neighbor);

dump_bit (*rover);

}

}

static int count_ones (void)

{

int count = 0;

int i, j;

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

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

if ((ising_lattice}

}

}

void ising_flip (bit *b)

{

int old_energy = b->energy;

int new_energy = energy (complement (b));

int delta_energy = new_energy - old_energy;

if (delta_energy > 0) {

double ising_threshold = exp (-(beta * delta_energy));

double random_value = rand () % scale / (double) scale;

if (ising_threshold < random_value) {

flip (b);

}

}

else {

flip (b);

}

}

void iterate (void)

{

int x = (int) rand () % dimension;

int y = (int) rand () % dimension;

ising_flip (&ising_lattice [x] [y]);

calculate_energies ();

}

void print_iteration (const char *format_string, int count)

{

size_t i;

size_t length = strlen (format_string);

putchar ('\n');

printf (format_string, count);

putchar ('\n');

for (i = 0; i < length + 3; i++)

putchar ('-');

putchar ('\n');

print_lattice (ising_lattice);

}

/* Housekeeping */

static void dump_neighbors (list_t *neighbors)

{

bit *rover;

lnode_t *neighbor;

for (neighbor = list_first (neighbors); neighbor != NULL; \

neighbor = list_next (neighbors, neighbor)) {

rover = lnode_get (neighbor);

dump_bit (*rover);

}

}

static int count_ones (void)

{

int count = 0;

int i, j;

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

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

if ((ising_lattice

*[j]).bitvalue == true) {*

count++;

}

}

}

return count;

}

static int count_zeroes (void)

{

int count = 0;

int i, j;

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

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

if ((ising_latticecount++;

}

}

}

return count;

}

static int count_zeroes (void)

{

int count = 0;

int i, j;

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

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

if ((ising_lattice

*[j]).bitvalue == false) {*

count++;

}

}

}

return count;

}

/* Error checking */

void minority_error (void)

{

fputs ("minority_size <argv [3]> is not a nonnegative integer\n", stderr);

exit (2);

}

void beta_error (void)

{

fputs ("beta <argv [2]> is not positive floating-point number.\n", stderr);

exit (2);

}

void dimension_error (void)

{

fputs ("dimension <argv [1]> is not a nonnegative integer.\n", stderr);

exit (2);

}

/* Return codes:

* 0 Successful execution

* 1 Memory allocation failed

* 2 Wrong form of argument

*/

int main (int argc, char *argv [])

{

if (argc > 3) {

if (is_positive_integer (argv [3]) || is_all_zeroes (argv [3])) {

minority_size = atoi (argv [3]);

}

else {

minority_error ();

}

if (minority_size < 0) {

minority_size = 0;

}

else if (minority_size > dimension * dimension) {

minority_size = dimension * dimension;

}

else if (minority_size > (dimension * dimension) / 2) {

minority_size = (dimension * dimension) - minority_size;

}

}

if (argc > 2) {

char *beta_string = argv [2];

size_t length = strlen (beta_string);

if (how_many (beta_string, '.') == 1) {

int position_decimal = find (beta_string, '.', 0);

char *integral_part = NULL;

char *fractional_part = NULL;

integral_part = copy_substring (integral_part, \

beta_string, \

0, \

position_decimal - 1);

if (is_positive_integer (integral_part) || \

is_all_zeroes (integral_part)) {

fractional_part = copy_substring (fractional_part, \

beta_string, \

position_decimal + 1, \

length - 1);

if (is_positive_integer (fractional_part) || \

is_all_zeroes (fractional_part)) {

beta = atof (argv [2]);

}

else {

beta_error ();

}

free (fractional_part);

}

else {

beta_error ();

}

free (integral_part);

}

else {

beta_error ();

}

}

if (argc > 1) {

if (is_positive_integer (argv [1]))

dimension = atoi (argv [1]);

else {

dimension_error ();

}

}

srand (time (0));

int i, j;

bit b;

/* Allocate the square Ising lattice */

ising_lattice = allocate_lattice (ising_lattice);

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

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

b.x = i;

b.y = j;

b.energy = 0;

b.bitvalue = random_bit ();

ising_latticecount++;

}

}

}

return count;

}

/* Error checking */

void minority_error (void)

{

fputs ("minority_size <argv [3]> is not a nonnegative integer\n", stderr);

exit (2);

}

void beta_error (void)

{

fputs ("beta <argv [2]> is not positive floating-point number.\n", stderr);

exit (2);

}

void dimension_error (void)

{

fputs ("dimension <argv [1]> is not a nonnegative integer.\n", stderr);

exit (2);

}

/* Return codes:

* 0 Successful execution

* 1 Memory allocation failed

* 2 Wrong form of argument

*/

int main (int argc, char *argv [])

{

if (argc > 3) {

if (is_positive_integer (argv [3]) || is_all_zeroes (argv [3])) {

minority_size = atoi (argv [3]);

}

else {

minority_error ();

}

if (minority_size < 0) {

minority_size = 0;

}

else if (minority_size > dimension * dimension) {

minority_size = dimension * dimension;

}

else if (minority_size > (dimension * dimension) / 2) {

minority_size = (dimension * dimension) - minority_size;

}

}

if (argc > 2) {

char *beta_string = argv [2];

size_t length = strlen (beta_string);

if (how_many (beta_string, '.') == 1) {

int position_decimal = find (beta_string, '.', 0);

char *integral_part = NULL;

char *fractional_part = NULL;

integral_part = copy_substring (integral_part, \

beta_string, \

0, \

position_decimal - 1);

if (is_positive_integer (integral_part) || \

is_all_zeroes (integral_part)) {

fractional_part = copy_substring (fractional_part, \

beta_string, \

position_decimal + 1, \

length - 1);

if (is_positive_integer (fractional_part) || \

is_all_zeroes (fractional_part)) {

beta = atof (argv [2]);

}

else {

beta_error ();

}

free (fractional_part);

}

else {

beta_error ();

}

free (integral_part);

}

else {

beta_error ();

}

}

if (argc > 1) {

if (is_positive_integer (argv [1]))

dimension = atoi (argv [1]);

else {

dimension_error ();

}

}

srand (time (0));

int i, j;

bit b;

/* Allocate the square Ising lattice */

ising_lattice = allocate_lattice (ising_lattice);

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

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

b.x = i;

b.y = j;

b.energy = 0;

b.bitvalue = random_bit ();

ising_lattice

*[j] = b;*

}

}

calculate_energies ();

/* Begin simulation */

print_iteration (initial_format, 0);

long k = 1;

while (true) {

if ((count_zeroes () == minority_size) || \

(count_ones () == minority_size)) {

break;

}

iterate ();

print_iteration (iteration_format, k++);

}

print_iteration (final_format, k);

/* End simulation */

free_lattice (ising_lattice);

exit (0);

}

/* ising.c ends */

The answer mentioned earlier:

* Ron Maimon (http://physics.stackexchange.com/users/4864/ron-maimon),

How do you start self-learning physics, URL (version: 2012-07-31):

http://physics.stackexchange.com/q/33223

Udyant Wig, programming novice}

}

calculate_energies ();

/* Begin simulation */

print_iteration (initial_format, 0);

long k = 1;

while (true) {

if ((count_zeroes () == minority_size) || \

(count_ones () == minority_size)) {

break;

}

iterate ();

print_iteration (iteration_format, k++);

}

print_iteration (final_format, k);

/* End simulation */

free_lattice (ising_lattice);

exit (0);

}

/* ising.c ends */

The answer mentioned earlier:

* Ron Maimon (http://physics.stackexchange.com/users/4864/ron-maimon),

How do you start self-learning physics, URL (version: 2012-07-31):

http://physics.stackexchange.com/q/33223

Udyant Wig, programming novice