package simulate;

import java.util.Random;

public class IntegratorGEMC extends Integrator {
    
    private final Random rand = new Random();
    public double maxRStep, maxVStep;
    public double pressure, betaMu, eBetaMu;
    private int freqDisplace, freqVolume, freqMolecule;
    private int iDisplace, iVolume, iMolecule, iTotal;
    public Phase secondPhase;
    private MCMove atomDisplace = new MCMoveAtom();
    
    public IntegratorGEMC() {
        super();
        freqDisplace = 1;
        freqVolume = freqMolecule = 1;
        maxRStep = 0.1;
        maxVStep = 0.1;
        phaseCountMax = 2;
//        super.phaseCountMax = 2;
        phase = new Phase[phaseCountMax];
        atomDisplace.setAdjustInterval(10000);
        atomDisplace.parentIntegrator = this;
    }
    
  public void registerPhase(Phase p) {
    super.registerPhase(p);
    if(phaseCount > phaseCountMax) {return;}
    if(phaseCount == 2) secondPhase = phase[1];
    System.out.println("phaseCount "+phaseCount);
  }
  
  //Performs basic GEMC trial
    public void doStep() {
        int i = (int)(rand.nextDouble()*iTotal);  //select type of trial to perform
        if(i < iDisplace) {                       //displace a randomly selected atom
            if(rand.nextDouble() < 0.5) {atomDisplace.doTrial(firstPhase);}
            else {atomDisplace.doTrial(secondPhase);}
        }
        else if(i < iVolume) {                    //volume exchange
            trialVolume();
        }
        else {                                    //molecule exchange
            if(rand.nextDouble() < 0.5) {trialExchange(firstPhase,secondPhase,Simulation.firstSpecies);}
            else {trialExchange(secondPhase,firstPhase,Simulation.firstSpecies);}
        }
    }
        
    private void trialVolume() {
        double v1Old = firstPhase.volume();
        double v2Old = secondPhase.volume();
        double hOld = firstPhase.energy.potential() + secondPhase.energy.potential();
        double vStep = (2.*rand.nextDouble()-1.)*maxVStep;
        double v1New = v1Old + vStep;
        double v2New = v2Old - vStep;
        double v1Scale = v1New/v1Old;
        double v2Scale = v2New/v2Old;
        double r1Scale = Math.pow(v1Scale,1.0/(double)Simulation.D);
        double r2Scale = Math.pow(v2Scale,1.0/(double)Simulation.D);
        firstPhase.inflate(r1Scale);
        secondPhase.inflate(r2Scale);
        double hNew = firstPhase.energy.potential() + secondPhase.energy.potential();
        if(hNew >= Double.MAX_VALUE ||
             Math.exp(-(hNew-hOld)/temperature+
                       firstPhase.moleculeCount*Math.log(v1Scale) +
                       secondPhase.moleculeCount*Math.log(v2Scale))
                < rand.nextDouble()) 
            {  //reject
              firstPhase.inflate(1.0/r1Scale);
              secondPhase.inflate(1.0/r2Scale);
            }
    }

    /**
     * Performs a GEMC molecule-exchange trial
     */
    private void trialExchange(Phase iPhase, Phase dPhase, Species species) {
        
        Species.Agent iSpecies = species.getAgent(iPhase);  //insertion-phase speciesAgent
        Species.Agent dSpecies = species.getAgent(dPhase);  //deletion-phase species Agent
        double uNew, uOld;
        
        if(dSpecies.nMolecules == 0) {return;}  //no molecules to delete; trial is over
        
        Molecule m = dSpecies.randomMolecule();  //select random molecule to delete
        uOld = dPhase.energy.meterPotential().currentValue(m); //get its contribution to energy
        m.displaceTo(iPhase.randomPosition());         //place at random in insertion phase
        uNew = iPhase.energy.meterPotential().currentValue(m); //get its new energy
        if(uNew == Double.MAX_VALUE) {  //overlap; reject
            m.replace();                //put it back 
            return;        
        }        
        double bFactor = dSpecies.nMolecules/dPhase.volume()  //acceptance probability
                         * iPhase.volume()/(iSpecies.nMolecules+1)
                         * Math.exp(-(uNew-uOld)/temperature);
        if(bFactor > 1.0 || bFactor > rand.nextDouble()) {  //accept
            iPhase.addMolecule(m,iSpecies);                 //place molecule in phase
        }
        else {              //reject
            m.replace();    //put it back
        }
    }
           
    public final int getFreqVolume() {return freqVolume;}
    public final void setFreqVolume(int f) {freqVolume = f;}
    public final int getFreqMolecule() {return freqMolecule;}
    public final void setFreqMolecule(int f) {freqMolecule = f;}
    public final int getFreqDisplace() {return freqDisplace;}
    public final void setFreqDisplace(int f) {freqDisplace = f;}
        
    public final double getMaxRStep() {return maxRStep;}
    public final void setMaxRStep(double s) {maxRStep = s;}
    
    public final double getMaxVStep() {return maxVStep;}
    public final void setMaxVStep(double s) {maxVStep = s;}
    
    public void initialize() {
        deployAgents();
        iDisplace = freqDisplace * firstPhase.moleculeCount;
        iVolume = iDisplace + freqVolume;
        iMolecule = iVolume + freqMolecule*firstPhase.moleculeCount;    
        iTotal = iMolecule;
    }
    
    public Integrator.Agent makeAgent(Atom a) {
        return new Agent(a);
    }
    
    public class Agent implements Integrator.Agent {
        public Atom atom;
        public Agent(Atom a) {atom = a;}
    }

}