/* 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) CNRS/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.multicanonical; import fr.ensea.montecarlo.data.*; import fr.ensea.chart.*; /** * A factory that produces reweighted statistical averages from a sample set and a weight w(E) obtained from S(E). */ public class MomentsReweighter { private WLDOSHistogram SE; private int VOL; private SamplesBag samplesBag; private Samples energySamples; private Samples magnetizationSamples; private ThermalAveragesSet set; public DefaultCurve2D meanEnergyVsKTCurve; // public DefaultCurve2D meanMagVsKTCurve; // public DefaultCurve2D susceptibilityVsKTCurve; // chi /** * reweighting w/o using flat histogram. */ public MomentsReweighter(WLDOSHistogram SE, SamplesBag samplesBag, ThermalAveragesSet set, double[] kTarray, int VOL){ this.SE = SE; this.samplesBag = samplesBag; this.energySamples = samplesBag.getEnergySamples(); this.magnetizationSamples = samplesBag.getMagnetizationSamples(); this.VOL = VOL; this.set=set; this.meanEnergyVsKTCurve = set.meanEnergyVsKTCurve; this.meanEnergyVsKTCurve.clear(); this.meanMagVsKTCurve = set.meanMagVsKTCurve; this.meanMagVsKTCurve.clear(); this.susceptibilityVsKTCurve = set.susceptibilityVsKTCurve; this.susceptibilityVsKTCurve.clear(); for (int i=0; i = 1/Z * sum_E (E) N(E) w(E), double sumE, sumM, sumM2; sumE=sumM=sumM2=0; // 3a°) accumulate sum_states A(state) w(E(state)) ... for (int imeas=0; imeas < samplesBag.size(); imeas++){ double E = energySamples.getSample(imeas)*VOL; // lattice energy (energySamples stores E/N !) double W = mw.getWeight(E); // w(E(state)) // Energy Per Spin: (E_ID=0) double e = energySamples.getSample(imeas); double e2 = e * e; sumE += e * W; // Magnetization (per spin, as usual): (M_ID=1) double M = magnetizationSamples.getSample(imeas); double M2 = M * M; sumM += M * W; sumM2 += M2 * W; } // 3b°) normalize, then fills moments[i], i.e. for this temperature point : double mean,variance; //System.out.print("kT="+kT); // E: mean = sumE/Z; //variance = sumE2/Z - sumE/Z * sumE/Z; if (!Double.isNaN(mean) && !Double.isInfinite(mean)) meanEnergyVsKTCurve.addXY(kT, mean); //System.out.print("\t="+mean); // M and chi mean = sumM/Z; if (!Double.isNaN(mean) && !Double.isInfinite(mean)) meanMagVsKTCurve.addXY(kT, mean); variance = sumM2/Z - sumM/Z * sumM/Z; if (!Double.isNaN(variance) && !Double.isInfinite(variance)) susceptibilityVsKTCurve.addXY(kT, VOL * variance / kT);// chi //System.out.print("\t="+mean); //System.out.println("\tchi="+(VOL * variance / kT)); set.updateGUI(); } ////////////////// PRIVATE METHODS ///////////////////// /** * Compute Z * exp(-K) for the given temperature w/o flat-histogram, where K is here to avoid overflows. * @param kT temperature point * @param mw Emu(E) reweighter */ private double computeZ(double kT, MuCaWeighter mw){ // compute Z = sum_states w(E(state)) ; actually what we compute is Z * exp(-K), where K is the normalization // constant, since otherwise Z would be too prone to overflowing ! double Z = 0.0; for (int imeas=0; imeas < samplesBag.size(); imeas++){ Z += mw.getWeight(energySamples.getSample(imeas)*VOL); // Z += w(E) for each state } return Z; } }