package etomica.virial;

import etomica.Phase;
import etomica.potential.Potential0;

/**
 * @author David Kofke
 *
 * Pair potential given according to the Mayer bonds in a cluster integral.
 * Does not require that the value of the cluster is non-negative.
 */

/* History
 * 08/20/03 (DAK) small changes to energy method (check for g = 0; abs(g)->g in
 * log argument
 * 08/21/03 (DAK) invoke resetPairs for pairSet in pi method
 * 12/16/03 (DAK) added field to hold indication of whether cluster giving value
 * of potential is positive or negative
 */
public class P0Cluster extends Potential0 {

	protected ClusterAbstract cluster;
	private boolean isPositive;
//	protected double g = 0;
	/**
	 * Constructor for P0Cluster.
	 * @param parent
	 */
	public P0Cluster(SimulationElement parent) {
		super(parent);
	}
	public P0Cluster(SimulationElement parent, ClusterAbstract cluster) {
		this(parent);
		setCluster(cluster);
	}
	
	public void setCluster(ClusterAbstract cluster) {
		this.cluster = cluster;
	}
	/**
	 * Returns the cluster.
	 * @return ClusterAbstract
	 */
	public ClusterAbstract getCluster() {
		return cluster;
	}
	
	public double pi(PhaseCluster phase) {
		double pi = cluster.value(phase.getPairSet().resetPairs(), 1.0/phase.integrator().temperature());
		isPositive = (pi > 0);
		return isPositive ? pi : -pi;
	}
//	public double pi(PairSet pairSet, double beta) {
//		double pi = cluster.value(pairSet, beta);
//		return (pi>0) ? pi : -pi;
//	}

//	public void calculate(Phase phase, IteratorDirective id, PotentialCalculation pc) {
//	   if(!enabled) return;
//	   double beta = 1.0/phase.integrator().temperature();
//	   PairSet pairs = ((PhaseCluster)phase).getPairSet();
//	   switch(id.atomCount()) {
//	   		case 0: g = cluster.value(pairs, beta); break;
//	   		case 1: g = cluster.value(id.atom1(), pairs, beta); break;
//	   		default: throw new RuntimeException();
//	   }
//	   ((PotentialCalculationEnergySum)pc).set(this).actionPerformed(phase);
//	}//end of calculate

	public double energy(Phase phase) {
		double g = pi((PhaseCluster)phase);
		return (g==1.0) ? 0.0 : ( (g==0) ? Double.POSITIVE_INFINITY : -phase.integrator().temperature()*Math.log(g) );
	}

	/**
	 * Returns true if the cluster value calculated during most recent call to
	 * the pi method was positive; returns false otherwise.
	 * @return boolean
	 */
	public boolean mostRecentIsPositive() {return isPositive;}

}

@@ -33,20 +33,63 @@
 etomica.
-Phas
+Atom;%0Aimport etomica.Phase;%0Aimport etomica.Spac
 e;%0Aimpor
@@ -641,98 +641,77 @@
  %7B%0A%0A
-%09protected ClusterAbstract cluster;%0A%09private boolean isPositive;%0A//%09protected double g = 0
+    private double temperature;%0A    private PhaseCluster phaseCluster
 ;%0A%09/
@@ -748,26 +748,8 @@
 er.%0A
-%09 * @param parent%0A
 %09 */
@@ -768,39 +768,26 @@
 luster(S
-imulationElement parent
+pace space
 ) %7B%0A%09%09su
@@ -794,28 +794,29 @@
 per(
+s
 pa
-rent
+ce
 );%0A%09%7D%0A
+%09%0A
 %09public 
 P0Cl
@@ -815,397 +815,90 @@
 lic 
-P0Cluster(SimulationElement parent, ClusterAbstract cluster) %7B%0A%09%09this(parent);%0A%09%09setCluster(cluster);%0A%09%7D%0A%09%0A%09public void setCluster(ClusterAbstract cluster) %7B%0A%09%09this.cluster = cluster;%0A%09%7D%0A%09/**%0A%09 * Returns the cluster.%0A%09 * @return ClusterAbstract%0A%09 */%0A%09public ClusterAbstract getCluster() %7B%0A%09%09return cluster;%0A%09%7D%0A%09%0A%09public double pi(PhaseCluster phase) %7B%0A%09%09double pi = cluster.value(phase
+double weight() %7B%0A%09%09return phaseCluster.getSampleCluster().value(phaseCluster
 .get
+C
 Pair
@@ -906,46 +906,12 @@
 et()
-.resetPairs(), 1.0/phase.integrator().
+, 1/
 temp
@@ -921,746 +921,256 @@
 ture
-()
 );%0A%09
-%09isPositive = (pi %3E 0);%0A%09%09return isPositive ? pi : -pi;%0A%09%7D%0A//%09public double pi(PairSet pairSet, double beta) %7B%0A//%09%09double pi = cluster.value(pairSet, beta);%0A//%09%09return (pi%3E0) ? pi : -pi;%0A//%09%7D%0A%0A//%09public void calculate(Phase phase, IteratorDirective id, PotentialCalculation pc) %7B%0A//%09   if(!enabled) return;%0A//%09   double beta = 1.0/phase.integrator().temperature();%0A//%09   PairSet pairs = ((PhaseCluster)phase).getPairSet();%0A//%09   switch(id.atomCount()) %7B%0A//%09   %09%09case 0: g = cluster.value(pairs, beta); break;%0A//%09   %09%09case 1: g = cluster.value(id.atom1(), pairs, beta); break;%0A//%09   %09%09default: throw new RuntimeException();%0A//%09   %7D%0A//%09   ((PotentialCalculationEnergySum)pc).set(this).actionPerformed(phase);%0A//%09%7D//end of calculate%0A%0A%09
+%7D%0A%0A    public void setPhase(Phase phase) %7B%0A    %09phaseCluster = (PhaseCluster)phase;%0A    %7D%0A    %0A    public void setTemperature(double aTemperature) %7B%0A        temperature = aTemperature;%0A    %7D%0A%0A    /**%0A     * @deprecated use weight()%0A     */%0A    
 publ
@@ -1186,27 +1186,28 @@
  energy(
-Phase phase
+Atom%5B%5D atoms
 ) %7B%0A%09%09do
@@ -1215,17 +1215,13 @@
 ble 
-g
+w
  = p
-i((P
 hase
@@ -1231,53 +1231,98 @@
 ster
-)phase);%0A%09%09return (g==1.0) ? 0.0 : ( (g
+.getSampleCluster().value(phaseCluster.getCPairSet(), 1/temperature);%0A%09%09if (w 
 ==
+ 
 0) 
-?
+return
  Dou
@@ -1346,286 +1346,54 @@
 NITY
- : -phase.integrator().temperature()*Math.log(g) );%0A%09%7D%0A%0A%09/**%0A%09 * Returns true if the cluster value calculated during most recent call to%0A%09 * the pi method was positive; returns false otherwise.%0A%09 * @return boolean%0A%09 */%0A%09public boolean mostRecentIsPositive() %7Breturn isPositive;%7D%0A
+;%0A%09%09return -Math.log(w)*temperature;%0A    %7D%0A    
 %0A%7D%0A

update P0Cluster. add "weight" method, which returns sampling weight (actual value of sampling cluster). sampling cluster is obtained from PhaseCluster