/* This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

package etomica.virial;

import etomica.atom.IAtomList;
import etomica.integrator.mcmove.MCMoveAtom;
import etomica.space.Space;
import etomica.space.Vector;
import etomica.util.random.IRandom;

/**
 * Class that samples molecule positions based on a chain of hard spheres of
 * diameter sigma.  The sequence of spheres is chosen randomly.
 *
 * @author Andrew
 */
public class MCMoveClusterAtomHSChain extends MCMoveAtom {

    protected final double sigma;
    protected int[] seq;

    public MCMoveClusterAtomHSChain(IRandom random, Space _space, double sigma) {
        super(random, null, _space);
        this.sigma = sigma;
    }

    public boolean doTrial() {

        IAtomList leafAtoms = box.getLeafList();
        int n = leafAtoms.size();
        if (seq == null) {
            seq = new int[n];
        }
        for (int i = 0; i<n; i++) {
            seq[i] = i;
        }
        for (int i = 0; i<n; i++) {
            int j = i+random.nextInt(n-i);
            int k = seq[j];
            seq[j] = seq[i];
            seq[i] = k;
        }
        leafAtoms.get(seq[0]).getPosition().E(0);

        for (int i = 1; i<n; i++) {
            Vector pos = leafAtoms.get(seq[i]).getPosition();

            pos.setRandomInSphere(random);
            double sig = getSigma(seq[i - 1], seq[i]);
            if (sig < 0) {
                // we want to force the position to be in the well (between 1 and sigma)
                sig = -sig;
                while (pos.squared() < 1 / (sig * sig)) {
                    pos.setRandomInSphere(random);
                }
            }
            pos.TE(sig);
            pos.PE(leafAtoms.get(seq[i-1]).getPosition());
        }

        ((BoxCluster)box).trialNotify();
        return true;
    }

    protected double getSigma(int i, int j) {
        return sigma;
    }

    public double getChi(double temperature) {
        return 1;
    }

    public void rejectNotify() {
        throw new RuntimeException("nope");
    }

    public void acceptNotify() {
        ((BoxCluster)box).acceptNotify();
    }
}

@@ -1409,17 +1409,91 @@
 tion();%0A
-%0A
+            boolean insideBox = false;%0A            while(!insideBox) %7B%0A    
         
@@ -1539,16 +1539,20 @@
         
+ 
+   
 double s
@@ -1594,24 +1594,28 @@
             
+    
 if (sig %3C 0)
@@ -1613,24 +1613,28 @@
 (sig %3C 0) %7B%0A
+    
             
@@ -1710,16 +1710,20 @@
  sigma)%0A
+    
         
@@ -1754,24 +1754,28 @@
             
+    
 while (pos.s
@@ -1800,24 +1800,28 @@
 g * sig)) %7B%0A
+    
             
@@ -1867,34 +1867,41 @@
                 
-%7D%0A
+    %7D%0A   
             %7D%0A  
@@ -1888,34 +1888,39 @@
 %0A               
+ 
 %7D%0A
+    
             pos.
@@ -1928,16 +1928,352 @@
 E(sig);%0A
+                insideBox = true;%0A                for (int j = 0; j %3C pos.getD(); j++) %7B%0A                    if (Math.abs(pos.getX(j)) %3E box.getBoundary().getBoxSize().getX(j)/2) %7B%0A                        insideBox = false;%0A                        break;%0A                    %7D%0A                %7D%0A                if (insideBox) %7B%0A        
         
@@ -2302,17 +2302,19 @@
 et(seq%5Bi
--
+ - 
 1%5D).getP
@@ -2325,16 +2325,48 @@
 ion());%0A
+                %7D%0A            %7D%0A
         

Handle case where chain is out of box