/* Monte-Carlo simulation code for statistical physics Copyright (C) 2001-2005 Sylvain Reynal Département de Physique Ecole Nationale Supérieure de l'Electronique et de ses Applications (ENSEA) 6, avenue du Ponceau, F-95014 CERGY CEDEX et Laboratoire de Physique Théorique et Modélisation (LPTM) Université de Cergy-Pontoise - Site de Neuville F-95031 CERGY CEDEX Tel : 00 +33 130 736 245 Fax : 00 +33 130 736 667 e-mail : reynal@ensea.fr web page : http://www.ensea.fr/staff/reynal/ */ package fr.ensea.montecarlo.canonical; import fr.ensea.montecarlo.data.*; import fr.ensea.montecarlo.model.*; import fr.ensea.montecarlo.canonical.*; import fr.ensea.montecarlo.misc.*; import java.io.*; import java.awt.Color; import java.text.*; import java.util.Locale; /** * Single-cluster Wolff algorithm for the 2D Potts model. */ public class PottsWolffCluster extends Montecarlo { private PottsLattice lattice; // --- RG --- private PottsLattice[] latticeRG; // non-null if RG is active private int factorRG; private Stack fifo; private static final int BUFFER_EMPTY=-1; // used by Stack private double bondProba; // 1 - exp(-beta) private int bonds; // nb of active bonds for the last flipped cluster private long totActiveBondsCumul; // incremented by "bonds" at each invokation of sweep() private int mcstepCounter; // incremented at each invokation of sweep() private int newQState, oldQState; // used by "tryNewBond()" private NumberFormat nfPercent; /** * Initialize a new cluster algorithm. */ public PottsWolffCluster(PottsLattice lattice, SamplesBag bag, double kT) { super(bag, kT); this.lattice = lattice; this.fifo = new Stack(); nfPercent = NumberFormat.getPercentInstance(Locale.US); nfPercent.setGroupingUsed(false); nfPercent.setMaximumFractionDigits(0); } public void activateRG(){ setFactorRG(1); } public void setFactorRG(int f){ if (f<1) f=1; this.factorRG=f; latticeRG = new PottsLattice[factorRG]; latticeRG[factorRG-1] = lattice.renormalize(latticeRG[factorRG-1]); for (int s=factorRG-2; s >= 0; s--){ latticeRG[s] = latticeRG[s+1].renormalize(latticeRG[s]); } } public PottsLattice getRGLattice(){ return latticeRG[0]; } public synchronized void setTemperature(double kT){ this.kT = kT; this.bondProba = 1.0 - Math.exp(-1.0/kT); reinitCounters(); } /** * Return a freshly generated random number b/w 0 and 1. */ public double createRandomNumber(){ return Math.random(); // [pending] use MersenneTwister } /** * Build then "flips" a cluster of spins */ public void sweep(){ final int qmax = lattice.pottsQ-1; mcstepCounter++; int spinIdx; // *) pick a spin at random spinIdx = (int)(lattice.size2 * createRandomNumber()); // b/w 0 and L^2-1 (raw() => slow) // *) save old q-state this.oldQState = lattice.getSpinValue(spinIdx); // *) pick a new spin state this.newQState = lattice.getSpinValue(spinIdx) + 1 + (int)(createRandomNumber() * qmax); if (this.newQState > qmax) this.newQState -= lattice.pottsQ; // assurons-nous que la nouvelle valeur est dans [0,Q-1] // *) flip the cluster initial seed: lattice.spins[spinIdx] = newQState; this.bonds=0; fifo.reset(); fifo.push(spinIdx); while((spinIdx=fifo.pop()) != BUFFER_EMPTY){ // -1 tryNewBond(lattice.neighborLeft(spinIdx)); tryNewBond(lattice.neighborRight(spinIdx)); tryNewBond(lattice.neighborAbove(spinIdx)); tryNewBond(lattice.neighborBelow(spinIdx)); } int flippedSpins = bonds+1; lattice.nbSpinsQ[oldQState] -= flippedSpins; lattice.nbSpinsQ[newQState] += flippedSpins; lattice.updateEnergy(); this.totActiveBondsCumul += bonds; samplesBag.addEnergySample(lattice.getEnergy()/lattice.size2); samplesBag.addMagnetizationSample(lattice.getMagnetization()); if (latticeRG != null){ // means RG has been activated latticeRG[factorRG-1] = lattice.renormalize(latticeRG[factorRG-1]); for (int s=factorRG-2; s >= 0; s--){ latticeRG[s] = latticeRG[s+1].renormalize(latticeRG[s]); } samplesBag.addEnergyRGSample(latticeRG[0].getEnergy()/latticeRG[0].size2); } } /** try and add a new spin to the cluster with probability "bondProba" if spins are identical */ private void tryNewBond(int nn){ if (oldQState==lattice.getSpinValue(nn) && createRandomNumber() < bondProba) { // oldQState = initial value of cluster seed lattice.spins[nn]=newQState; // flip spin bonds++; // increment B = number of active bonds fifo.push(nn); } } /** * Implementation of a FIFO buffer used to stores indices of spins belonging to the cluster */ class Stack { private int[] fifoBuffer; // FIFO buffer for cluster construction private int fifoPointer; // index of last added spin index ; FIFO grows upwards Stack(){ fifoBuffer = new int[400]; // 20x20 reset(); } void reset(){ fifoPointer=-1; } void push(int i){ if (fifoPointer==fifoBuffer.length-1){ // stack overflow ! int[] _fifoBuffer = new int[fifoBuffer.length*2]; System.arraycopy(fifoBuffer,0,_fifoBuffer,0,fifoBuffer.length); fifoBuffer = _fifoBuffer; } fifoBuffer[++fifoPointer]=i; } /** returns BUFFER_EMPTY whenever buffer is empty */ int pop(){ if (fifoPointer < 0) return BUFFER_EMPTY; // -1 else { int i = fifoBuffer[fifoPointer]; fifoPointer--; return i; } } public String toString(){ StringBuffer s = new StringBuffer(); s.append('['); for (int i=0; i<=fifoPointer; i++){ s.append(fifoBuffer[i]); s.append(' '); } s.append(']'); return s.toString(); } } /////////////////////////////////////////////////////////// /// Statistics /////////////////////////////////////////////////////////// public double getFractionActiveBonds(){ return (double)totActiveBondsCumul / (double)(mcstepCounter * lattice.size2 * 2); } public synchronized void reinitCounters(){ this.totActiveBondsCumul = this.mcstepCounter = 0; } public String getCountersRateLabel(){ return "Active bonds="+nfPercent.format(getFractionActiveBonds()); } ///////////////////////////////////////////////////////////////////////////////// // TEST ///////////////////////////////////////////////////////////////////////////////// public String toString(){ String s = "MC Wolff algorithm :\n"; s += "\tmean acceptance rate = " + nfPercent.format(getFractionActiveBonds()) + "\n"; return s; } public static void main(String args[]){ int L=10, Q=3; PottsLattice lattice = new PottsLattice(L,Q); SamplesBag bag = new SamplesBag(0); PottsWolffCluster algo = new PottsWolffCluster(lattice, bag, 0.5); for (double kT=0.5; kT<2.0; kT+=0.2){ System.out.println("---------- kT="+kT+" ----------"); algo.setTemperature(kT); final int nbPas = 10000; for (int i=0; i<1000; i++) algo.sweep(); for (int i=0; i="+meanE); System.out.println("\t\tCv="+(varE/algo.kT/algo.kT/L)); System.out.print("="+meanM); System.out.println("\tChi="+(L*varM/algo.kT)); System.out.println("Fraction of active bonds =" + (100.0*algo.getFractionActiveBonds())+" %"); } } }