/* Monte-Carlo simulation code for statistical physics Copyright (C) 2001-2004 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.model; import fr.ensea.chart.*; import fr.ensea.math.*; import fr.ensea.montecarlo.canonical.*; import fr.ensea.montecarlo.multicanonical.*; import fr.ensea.montecarlo.data.*; import java.util.*; /** * This class represents a lattice spin model with nearest-neighbor interactions, and helical boundary * conditions. Each spin can take on Q integer values b/w 0 and Q-1 (Potts model). The special value * Q=2 gives the Ising model. * * This class holds the state of each spin, the lattice energy, the number of states of the model, * and the number of spins per phase (this is used to compute the order parameter). * * Helical conditions are implemented by storing spins in a one-dimensional array, from left to right * and from top to bottom. Thus the "right" neighbor (resp. "left") of a spin belonging to the last * column is the first (resp. last) one of the next row. However, vertical boundary conditions are * implemented in the standard way (PBC's, aka Born-Von-Karman). Helical boundary conditions yield * a lighter implementation, and are known to reduce sensitivity to finite lattice size. * [FR] * Cette classe représente un modèle de spins sur réseau avec interactions entre plus proches voisins, * et des conditions aux limites hélicoïdales. Chaque spin peut prendre Q valeurs entre 0 et Q-1 (modèle de Potts). * On retrouve le modèle d'Ising en imposant Q=2. * * Sont stockés dans cette classe : l'état de chaque spin et l'énergie totale du réseau, le nombre d'état Q * du modèle, et le nombre de spins par phase, cette dernière variable étant nécessaire au calcul du paramètre d'ordre. * * Les conditions aux limites hélicoïdales sont implémentées en rangeant l'ensemble des spins du réseau, * de gauche à droite et de haut en bas, dans un tableau unidimensionnel. Ainsi, le voisin "de droite" * (respectivement "de gauche") d'un spin situé sur la dernière colonne est le premier (resp. le dernier) * spin de la ligne suivante (resp. précédente). En revanche, les conditions aux limites verticales sont les * conditions périodiques habituelles (à la Born-Van-Karman). L'avantage des conditions aux limites hélicoïdales * est simplement d'alléger l'implémentation de l'algorithme. */ public class PottsLattice implements ILattice2D { // package access for use inside PottsWolffCluster etc... public int pottsQ; // nb of Q states public int size; // lattice size public int size2; // L*L public int[] spins; // array holding the state of each spin; helical boundary conditions are used public int[] nbSpinsQ; // nbSpins[s]=nb of spins in phase "s" (from 0 to Q-1) public double energy; // lattice energy public double externalField; // external magnetic field /** * Init a 2D square lattice with L=10 (i.e. 100 spins), with each spin inited at random b/w 0 and 2 (Q=3). */ public PottsLattice() { init(3,10); } /** * Init a 2D square lattice with the given size L, and each spin inited at random b/w 0 and Q-1 */ public PottsLattice(int Q, int L) { init(Q,L); } /** * Reinit lattice properties. * @param Q nb of states of the model * @param L lattice size */ public synchronized void init(int Q, int L){ pottsQ = Q; size = L; size2 = L*L; spins = new int[size2]; createRandomConfiguration(); updateEnergy(); } /** * Reinit lattice properties and spin values from the given source. * @param Q nb of states of the model * @param L lattice size */ public synchronized void init(PottsLattice src){ pottsQ = src.pottsQ; size = src.size; size2 = src.size2; spins = new int[src.size2]; System.arraycopy(src.spins,0, spins,0, spins.length); } /** * Creates a random configuration */ public synchronized void createRandomConfiguration(){ for(int i=0; i= size2) nn -= size2; return nn; } /** * Returns the index of the nearest-neighbor below the given spin */ public int neighborBelow(int idx){ int nn=idx-size; if (nn < 0) nn += size2; return nn; } /** * Returns the index of the nearest-neighbor above the given spin */ public int neighborAbove(int idx){ int nn=idx+size; if (nn >= size2) nn -= size2; return nn; } ///////////////////////////////////////////////////////////////// //// Properties ///////////////////////////////////////////////////////////////// /** * Return the lattice size (=side length) */ public int getLatticeSize(){ return size; } /** * Return the nb of spins (=side length squared) */ public int getSpinCount(){ return size2; } /** * Return the number of states of the model (2 for the Ising model) */ public int getPottsQ(){ return this.pottsQ; } /** * Returns the exact critical temperature */ public double getExactTc(){ return 1.0/Math.log(1+Math.sqrt(pottsQ)); } /** * Returns the lattice energy */ public double getEnergy(){ return energy; } /** * Computes then returns the magnetization */ public double getMagnetization(){ int maxQ=nbSpinsQ[0]; for (int q=1; q maxQ) maxQ = nbSpinsQ[q]; } double mag = (pottsQ * (double)maxQ / (double)size2 - 1.0)/(pottsQ - 1.0); return mag; } /** * Returns the value of the ith spin (with helical PBC's) */ public int getSpinValue(int i){ if (i>=spins.length) return 0; return spins[i]; } /** * Returns the value of the spin at (i,j) */ public int getSpinValue(int i, int j){ return spins[i*size+j]; } /** * Changes the value of the ith spin to val, then updates the nb of spins per phase. */ public void setSpinValue(int i, int val){ if (val >= this.pottsQ) val=pottsQ-1; nbSpinsQ[spins[i]]--; spins[i] = val; nbSpinsQ[val]++; } /** * Changes the value of the external magnetic field (forces condensation in phase s=0) */ public void setMagneticField(double fieldValue){ externalField = fieldValue; } /** * Returns the current value of the external magnetic field */ public double getMagneticField(){ return externalField; } /** * Returns the energy of the ground-state */ public double getGroundstateEnergy(){ return -this.size2*2; } ///////////////////////////////////////////////////////////////// //// Renormalization ///////////////////////////////////////////////////////////////// /** * returns the power of 2 that is immediately below (or equal to) the given number */ public static int floorToPowerOf2(final int x){ for(int i=1,y=1<<30;i<32;i++,y>>=1) { // generate 2^31,...,4,2,1 int z=x&y; // isolate one bit, starting from MSB if(z!=0) { // check if this bit is non-null => this yields the power of 2 just below x return y; } } return x; } /** * Generates a renormalized lattice using the majority rule and 2x2 plaquettes. * Note that the lattice size must be a power of 2. * @param result renormalized lattice (may be null, in which case it's allocated here) * @return result, for convenience. */ public PottsLattice renormalize(PottsLattice result){ if (size==1 || (((size>>1)<<1)!=size)) { // can't renormalize if size=1 or not even if (result==null) result = new PottsLattice(pottsQ,size); result.init(this); return result; } if (result==null || result.size != (size >>1)) result = new PottsLattice(pottsQ,size>>1); int spinDest=0; int spinSrc; int valRenorm; switch (pottsQ){ // --- Q=2 --- case 2: // from bottom to top and from left to right: // val3 val4 // val1 val2 // with val=0 or 1 // * valRenorm = sum(val1 à val4) = 0 to 4 // * if valRenorm==2, draw at random // * if valRenorm<2, set it to 0 // * if valRenorm>2, set it to 1 for (int row=0; row< size; row+=2){ for(int col=0; col < size; col+=2){ spinSrc=row*size+col; // majority rule: valRenorm = this.spins[spinSrc]; valRenorm += this.spins[neighborRight(spinSrc)]; valRenorm += this.spins[neighborAbove(spinSrc)]; valRenorm += this.spins[neighborRight(neighborAbove(spinSrc))]; if (valRenorm > 2) valRenorm=1; else if (valRenorm < 2) valRenorm=0; else valRenorm = (int)(2*Math.random()); result.spins[spinDest] = valRenorm; spinDest++; } } break; // --- Q>2 --- default: //System.out.println("========== RG:Q>2 =========="); // from bottom to top and from left to right: // val3 val4 // val1 val2 // with val=0 or Q-1 // * valRenorm = maj(val1,val2,val3,val4), and we draw at random if undetermined. int[] nbSpinsQRG = new int[pottsQ]; for(int row=0; row < size; row+=2){ for (int col=0; col< size; col+=2){ //System.out.println("[col,row]=["+col+","+row+"]"); // reset array for (int i=0; i -1 ! } } if (valRenorm==-1){ // "valRenorm" not set yet // {X X , X X} or {X X , X Y} or {X X , Z Y} // note that {X X , Y Y} was handled in the for() loop just above if (maxSpinCountPerPhase > 1) valRenorm=majorityPhase; else { // {X Y, Z T} => draw one spin value at random int b = (int)(4.0*Math.random()); // 0, 1, 2 or 3 switch (b){ case 0: valRenorm = this.spins[spinSrc]; break; case 1: valRenorm = this.spins[neighborRight(spinSrc)]; break; case 2: valRenorm = this.spins[neighborAbove(spinSrc)]; break; case 3: valRenorm = this.spins[neighborRight(neighborAbove(spinSrc))]; break; default: } } } //System.out.println("valRenorm="+valRenorm); result.spins[spinDest] = valRenorm; spinDest++; } } break; } result.updateEnergy(); result.updateNbSpinsPerPhase(); //System.out.println("lattice = " + toString()); //System.out.println("RG lattice = " + result.toString()); return result; } ///////////////////////////////////////////////////////// //// Metropolis ///////////////////////////////////////////////////////// /** * Compute the energy difference b/w this lattice and the same lattice with spin "spinIdx" moved to "newValue" */ public double getEnergyChange(int spinIdx, int newVal){ double dE=0; int oldVal = getSpinValue(spinIdx); if (newVal==oldVal) return 0; int neighborVal = spins[neighborLeft(spinIdx)]; if (oldVal==neighborVal) dE++; if (newVal==neighborVal) dE--; neighborVal = spins[neighborRight(spinIdx)]; if (oldVal==neighborVal) dE++; if (newVal==neighborVal) dE--; neighborVal = spins[neighborAbove(spinIdx)]; if (oldVal==neighborVal) dE++; if (newVal==neighborVal) dE--; neighborVal = spins[neighborBelow(spinIdx)]; if (oldVal==neighborVal) dE++; if (newVal==neighborVal) dE--; // external field if (oldVal==0) dE += externalField; else if (newVal==0) dE -= externalField; return dE; } ////////////////////////////////////// //// algorithms ////////////////////////////////////// public Montecarlo createMetropolis(SamplesBag bag, double kT){ return new PottsMetropolis(this, bag, kT); } public Montecarlo createWolffCluster(SamplesBag bag, double kT){ return new PottsWolffCluster(this, bag, kT); } public WangLandau createWangLandau(SamplesBag bag, double weight){ if (this.getPottsQ()==2) return new PottsWL(this,bag,2.0, weight); else return new PottsWL(this,bag,1.0, weight); } ///////////////////////////////////////////////////////// //// debug ///////////////////////////////////////////////////////// public String toString(){ String s = "Modèle de Potts bidimensionnel : \n"; s += "* size de réseau=" + Integer.toString(size) + "\n"; s += "* M=" + getMagnetization() + "\n"; s += "* E=" + getEnergy() + "\n"; // affiche le contenu de chaque ligne du réseau : int i=0; for (int row = 0; row < size; row++){ for(int col = 0; col < size; col++){ s += Integer.toString(spins[i++]) + " "; } s += "\n"; } return s; } // for testing purpose public static void main(String[] args){ PottsLattice lat = new PottsLattice(4, 2); for (int i=0; i