#!/usr/bin/python -Wall # ================================================================ # John Kerl # kerl at math dot arizona dot edu # 2010-03-02 # ================================================================ # This software is released under the terms of the GNU GPL. # Please see LICENSE.txt in the same directory as this file. # ================================================================ from __future__ import division from math import * import sackmat_m import re, os, sys, random import Image # ================================================================ def print_spins_pm(spins, L): s = '' for i in range(0, L): for j in range(0, L): if (spins[i][j] == 1): #print '+', print '@', else: #print '-', print '.', print '' # ---------------------------------------------------------------- def print_spins(spins, L): print_spins_pm(spins, L) # ---------------------------------------------------------------- def spins_to_image_file(spins, L, image_file_name, pps=1): im = Image.new('RGB', (L*pps, L*pps)) up_rgb = (0x99, 0x33, 0x33) down_rgb = (0x40, 0x40, 0x40) for si in xrange(0, L): for sj in xrange(0, L): if spins[si][sj] > 0: rgb = up_rgb else: rgb = down_rgb for di in range(0, pps): for dj in range(0, pps): im.putpixel((si*pps+di,sj*pps+dj), rgb) im.save(image_file_name) print 'Wrote %s' % (image_file_name) # ================================================================ def set_random_initial_spins(spins, L, p_of_up): for i in range(0, L): for j in range(0, L): u = random.uniform(0.0, 1.0) if (u < p_of_up): spins[i][j] = 1 else: spins[i][j] = -1 # ================================================================ # State energy for the Ising model. def compute_energy(spins, L, c, h): Lm1 = L-1 E = 0.0 for i in range(0, L): for j in range(0, L): if i < Lm1: E -= c * spins[i][j] * spins[i+1][ j ] else: E -= c * spins[i][j] * spins[ 0 ][ j ] # PBC if j < Lm1: E -= c * spins[i][j] * spins[ i ][j+1] else: E -= c * spins[i][j] * spins[ i ][ 0 ] # PBC E += h * spins[i][j] return E # ---------------------------------------------------------------- def compute_mag(spins, L): M = 0.0 for i in xrange(0, L): for j in xrange(0, L): M += spins[i][j] return M # ---------------------------------------------------------------- def compute_site_corr(spins, L, site_corrs): for i in xrange(0, L): site_corrs[i] = spins[0][0] * spins[i][0] # ---------------------------------------------------------------- # DE = -2 sij * [ c * (L + R + D + U) + h]. def compute_energy_delta(spins, L, i, j, c, h): Lm1 = L-1 DE = 0.0 if i > 0: DE -= spins[i-1][ j ] # Left else: DE -= spins[Lm1][ j ] # PBC if i < Lm1: DE -= spins[i+1][ j ] # Right else: DE -= spins[ 0 ][ j ] # PBC if j > 0: DE -= spins[ i ][j-1] # Down else: DE -= spins[ i ][Lm1] # PBC if j < Lm1: DE -= spins[ i ][j+1] # Up else: DE -= spins[ i ][ 0 ] # PBC DE *= c DE += h DE *= -2.0 * spins[i][j] return DE # ================================================================ def metro_step(spins, L, i, j, c, h, verbosity, EM_list): action = 'reject' Delta_E = compute_energy_delta(spins, L, i, j, c, h) p = min(1.0, exp(-Delta_E)) u = random.uniform(0.0, 1.0) if (u < p): if spins[i][j] == 1: Delta_M = -2 else: Delta_M = 2 EM_list[0] += Delta_E EM_list[1] += Delta_M spins[i][j] = -spins[i][j] action = 'accept' if (verbosity >= 3): print ' metro step at site (%2d,%2d):' % (i, j) print_spins(spins, L) print 'DE=%11.7f' % (Delta_E), print 'E=%11.7f %s' % (EM_list[0], action) # ---------------------------------------------------------------- def metro_sweep(spins, L, c, h, verbosity, EM_list): for i in range(0, L): for j in range(0, L): metro_step(spins, L, i, j, c, h, verbosity, EM_list) # ================================================================ def usage(): print >> sys.stderr, 'Usage: %s [options ...]' % (sys.argv[0]) print >> sys.stderr, 'Options:' print >> sys.stderr, ' L=...: box length' print >> sys.stderr, ' h=...: magnetic-field strength' print >> sys.stderr, ' c=...: coupling strength' print >> sys.stderr, ' nsweep=...: number of metropolis sweeps' print >> sys.stderr, ' pup=[0.0-1.0]: probability of up for metro init state' print >> sys.stderr, ' mv={0-3}: metropolis verbosity' sys.exit(1) # ================================================================ # Start of script L = 200 h = 0.0 c = 0.35 nsweep = 500 pup = 1.0 mv = 0 pix_per_site = 1 argc = len(sys.argv) for argi in range(1, argc): arg = sys.argv[argi] if arg == '-h' or arg == '--help': usage() elif re.match(r'^L=', arg): L = int (arg[2:]) elif re.match(r'^h=', arg): h = float (arg[2:]) elif re.match(r'^c=', arg): c = float (arg[2:]) elif re.match(r'^nsw=', arg): nsweep = int (arg[4:]) elif re.match(r'^pup=', arg): pup = float (arg[4:]) elif re.match(r'^mv=', arg): mv = int (arg[3:]) elif re.match(r'^pps=', arg): pix_per_site = int (arg[4:]) else: usage() spins = sackmat_m.make_zero_matrix(L, L) set_random_initial_spins(spins, L, pup) site_corrs = [0.0] * L N = L*L # ---------------------------------------------------------------- # Print header '# c = %11.7f' % (c) '# h = %11.7f' % (h) E = compute_energy(spins, L, c, h) M = compute_mag(spins, L) EM_list = [E, M] taps = [[0, 0], [0,2], [0,4], [0,6], [0, 10], [0, 20], [1, 0], [1, nsweep-1]] #taps=[] # Pass 1: thermalization # Pass 2: data accumulation descs = ['thermalization', 'data accumulation'] for iter in [0, 1]: e_mean = 0.0 m_mean = 0.0 sc_means = [0.0] * L for sweepi in range(0, nsweep): if (sweepi % 50) == 0: print '# iter %d sweepi %6d' % (iter, sweepi) if (mv >= 2): print ' metro sweep %2d (%s):' % (sweepi, descs[iter]) os.system('clear') print_spins(spins, L) if 1: img_file_name = '%5.3f_%d_%06d.png' % (c, iter, sweepi) spins_to_image_file(spins, L, img_file_name, pix_per_site) else: for [tapiter, tapsweepi] in taps: if [tapiter, tapsweepi] == [iter, sweepi]: img_file_name = '%5.3f_%d_%06d.png' % (c, iter, sweepi) spins_to_image_file(spins, L, img_file_name, pix_per_site) #E = compute_energy(spins, L, c, h) #M = compute_mag(spins) #m = M/N e = EM_list[0]/N m = EM_list[1]/N compute_site_corr(spins, L, site_corrs) e_mean += e m_mean += m for i in xrange(0, L): sc_means[i] += site_corrs[i] if (mv >= 1): print '%11.7f %11.7f # e m' % (e, m) metro_sweep(spins, L, c, h, mv, EM_list) if (mv >= 3): print if (mv >= 2): print e_mean /= nsweep m_mean /= nsweep for i in xrange(0, L): sc_means[i] /= nsweep if iter == 1: print '#' print '# ebar = %11.7f' % (e_mean) print '# mbar = %11.7f' % (m_mean) for i in xrange(0, L): print '# sc %2d = %11.7f' % (i, sc_means[i])