/* 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) 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; /** * Returns exp(Emu(E)/kT0 - E/kT - antiOverflowK), where "antiOverflowK" enables us to avoid overflows. * * This is based on the better estimate of n(E), which is N(E)/PI_mu(E), where: * - PI_mu(E) is exp(-Emu(E)/kT0), * - N(E) is the energy histogram. * * Hence n(E) ~ N(E) exp(Emu(E)/kT0). Once reweighted at temperature kT, * we thus get : N(E) exp(Emu(E)/kT0 - E/kT). Hence the reweighting factor w(E) is given by exp(Emu(E)/kT0 - E/kT). */ public class MuCaWeighter { private double kT; private double antiOverflowK; private WLDOSHistogram hmue; // store Emu(E) /** * Construct a new MuCa weighter from the given S(E) */ public MuCaWeighter(double kT, WLDOSHistogram hmue){ this.kT = kT; this.hmue = hmue; this.antiOverflowK=0; computeAntiOverflowConstant(); } /** * Overflow work-around : max=Math.exp(709.78) and min=Math.exp(-744.44) * The idea is to first determine max{w(E)} (actually max{Emu(E)/kT0 - E/kT}), that is, log w(E), * and possibly renormalize to avoid overflow */ private void computeAntiOverflowConstant(){ double maxEmuCa = -Double.MAX_VALUE; //double dE = muca.getDeltaE(); [pending] double minEnergyFlat = hmue.xAxis().binCentre(0); double maxEnergyFlat = hmue.xAxis().binCentre(hmue.xAxis().bins()-1); double dE = (maxEnergyFlat-minEnergyFlat)/100.0; for (double ene=minEnergyFlat+dE; ene < maxEnergyFlat; ene+=dE){ if (getExponent(ene) == Double.POSITIVE_INFINITY) continue; // i.e. if S(E)==REJECT_MOVE maxEmuCa = Math.max(maxEmuCa, getExponent(ene)); } //System.out.println("Max Emuca = " + maxEmuCa); if (maxEmuCa >= Math.log(Double.MAX_VALUE)-40){ this.antiOverflowK = maxEmuCa - Math.log(Double.MAX_VALUE) + 40; // e.g. if maxEmuCa=800, antiOverflowK=121 } else this.antiOverflowK=0.0; // hence, now there can't be any overflow (yet underflows are still allowed!) } /** * return S(E) - E/kT ; allow computing max value for overflow workaround */ private double getExponent(double energy){ return this.hmue.getSE(energy)-energy/this.kT; } ////////////////////////////////////////////////////////////////////////////////////////// /// IWeighter inteface ////////////////////////////////////////////////////////////////////////////////////////// /** * Returns w(E,kT) = exp(S(E) - E/kT - antiOverflowK), except if E > Emax or E < Emin (i.e. flatness boundaries) * in which case 0.0 is returned. */ public double getWeight(double energy){ double SE = this.hmue.getSE(energy); return Math.exp(SE-energy/this.kT-this.antiOverflowK); } /** * return the normalization constant */ public double getAntiOverflowK(){ return antiOverflowK; } }