package etomica.data.meter;

import etomica.EtomicaInfo;
import etomica.action.MoleculeActionTranslateTo;
import etomica.api.IBox;
import etomica.api.IMolecule;
import etomica.api.IRandom;
import etomica.api.ISpecies;
import etomica.box.RandomPositionSource;
import etomica.box.RandomPositionSourceRectangular;
import etomica.data.DataSourceScalar;
import etomica.integrator.IntegratorBox;
import etomica.space.ISpace;
import etomica.units.Null;

/**
 * Meter to measure the chemical potential (as its exponent: exp(-mu/kT)) of a
 * species via the Widom insertion method. Call to getDataAsScalar returns a
 * Widom-insertion average (i.e., sum of exp(-energy/kT) using a configurable
 * number of insertion trials) for a molecule of the desired species in the
 * given box in its current configuration. It is expected that repeated calls
 * will be averaged to give an overall insertion average for many configurations
 * of the box.
 * <br>
 * Meter can be configured to measure the residual chemical potential -- above
 * the ideal-gas value -- or the full chemical potential.  If not configured for
 * the residual chemical potential, the insertion average will be multiplied by
 * the species number density in the box in its current configuration.  Default
 * behavior will give residual chemical potential.
 * 
 * @author David Kofke
 */
public class MeterWidomInsertion extends DataSourceScalar {

    public MeterWidomInsertion(ISpace space, IRandom random) {
        super("exp(-\u03BC/kT)", Null.DIMENSION);//"\u03BC" is Unicode for greek "mu"
        setNInsert(100);
        setResidual(true);
        atomTranslator = new MoleculeActionTranslateTo(space);
        positionSource = new RandomPositionSourceRectangular(space, random);
    }

    public static EtomicaInfo getEtomicaInfo() {
        EtomicaInfo info = new EtomicaInfo(
            "Chemical potential via Widom's ghost-particle insertion method");
        return info;
    }

    /**
     * Sets flag specifying if full or residual chemical potential is computed
     * Default is <code>true</code> (only residual is computed)
     */
    public void setResidual(boolean b) {
        residual = b;
    }

    /**
     * Accessor for flag specifying if full or residual chemical potential is
     * computed
     */
    public boolean isResidual() {
        return residual;
    }

    /**
     * Sets the species, takes a prototype molecule, and gets handle to
     * appropriate species agent in box
     */
    public void setSpecies(ISpecies s) {
        species = s;
        testMolecule = s.makeMolecule();
    }

    /**
     * Accessor for the species for which chemical potential is evaluated
     */
    public ISpecies getSpecies() {
        return species;
    }

    /**
     * Number of Widom insertions attempted with each call to currentValue
     */
    public void setNInsert(int n) {
        nInsert = n;
    }

    /**
     * Accessor to number of Widom insertions attempted with each call to
     * currentValue
     */
    public int getNInsert() {
        return nInsert;
    }

    /**
     * Performs a Widom insertion average, doing nInsert insertion attempts
     * Temperature used to get exp(-uTest/kT) is that of the integrator for the
     * box
     * 
     * @return the sum of exp(-uTest/kT)/nInsert, multiplied by V/N if 
     * <code>residual</code> is false
     */
    public double getDataAsScalar() {
        if (integrator == null) throw new IllegalStateException("must call setBox before using meter");
        IBox box = integrator.getBox();
        double sum = 0.0; //sum for local insertion average
        box.addMolecule(testMolecule);
        energyMeter.setTarget(testMolecule);
        for (int i = nInsert; i > 0; i--) { //perform nInsert insertions
            atomTranslator.setDestination(positionSource.randomPosition());
            atomTranslator.actionPerformed(testMolecule);
            double u = energyMeter.getDataAsScalar();
            sum += Math.exp(-u / integrator.getTemperature());
        }

        box.removeMolecule(testMolecule);

        if (!residual) {
            // multiply by V/N
            sum *= box.getBoundary().volume() / box.getNMolecules(species);
        }
        return sum / nInsert; //return average
    }

    /**
     * Returns the integrator associated with this class.  The box, potentialMaster
     * and temperature are taken from the integrator.
     */
    public IntegratorBox getIntegrator() {
        return integrator;
    }

    /**
     * Sets the integrator associated with this class.  The box, potentialMaster
     * and temperature are taken from the integrator.
     */
    public void setIntegrator(IntegratorBox newIntegrator) {
        integrator = newIntegrator;
        energyMeter = new MeterPotentialEnergy(integrator.getPotentialMaster());
        energyMeter.setBox(integrator.getBox());
        positionSource.setBox(integrator.getBox());
    }
    
    /**
     * Sets a new RandomPositionSource for this meter to use.  By default, a
     * position source is used which assumes rectangular boundaries.
     */
    public void setPositionSource(RandomPositionSource newPositionSource) {
        positionSource = newPositionSource;
        if (integrator != null) {
            positionSource.setBox(integrator.getBox());
        }
    }
    
    /**
     * Returns the RandomPositionSource used by this meter.
     */
    public RandomPositionSource getPositionSource() {
        return positionSource;
    }

    private static final long serialVersionUID = 1L;
    private IntegratorBox integrator;

    /**
     * Number of insertions attempted in each call to currentValue Default is
     * 100
     */
    private int nInsert;
    private ISpecies species;
    private IMolecule testMolecule;// prototype insertion molecule
    private boolean residual; // flag to specify if total or residual chemical
                              // potential evaluated. Default true
    private MoleculeActionTranslateTo atomTranslator;
    protected RandomPositionSource positionSource;
    private MeterPotentialEnergy energyMeter;
}

@@ -1311,16 +1311,139 @@
 al.%0A * %0A
+ * The actual chemical potential can be calculated as -kT ln(%3Cx%3E) where x is%0A * the value returned by getDataAsScalar.%0A * %0A
  * @auth

add a bit about how to actually get mu