<commit_before>package etomica.models.oneDHardRods;

import etomica.api.IAtomList;
import etomica.api.IAtomType;
import etomica.api.IBox;
import etomica.api.IRandom;
import etomica.api.IVectorMutable;
import etomica.box.Box;
import etomica.data.DataSourceScalar;
import etomica.data.meter.MeterPotentialEnergy;
import etomica.lattice.crystal.Basis;
import etomica.lattice.crystal.Primitive;
import etomica.nbr.list.PotentialMasterList;
import etomica.normalmode.CoordinateDefinition;
import etomica.normalmode.CoordinateDefinitionLeaf;
import etomica.normalmode.NormalModes;
import etomica.normalmode.NormalModes1DHR;
import etomica.normalmode.P2XOrder;
import etomica.normalmode.WaveVectorFactory;
import etomica.normalmode.CoordinateDefinition.BasisCell;
import etomica.potential.P2HardSphere;
import etomica.potential.Potential2;
import etomica.potential.Potential2HardSpherical;
import etomica.simulation.Simulation;
import etomica.space.Boundary;
import etomica.space.BoundaryRectangularPeriodic;
import etomica.species.SpeciesSpheresMono;
import etomica.units.Null;


/**
 * Uses a Widom-like insertion of a mode to calculate a probability.
 * Uses a different box than the main simulation, to assume a mode & rod 
 * are removed
 * 
 * @author cribbin
 *
 */
public class MeterDifferentImageSubtract extends DataSourceScalar {

    public int nInsert, counter;
    private MeterPotentialEnergy meterPE;
    private CoordinateDefinition cDef, simCDef;
    private int cDim, simCDim;
    private IVectorMutable[] waveVectors, simWaveVectors;
    private double[] simRealT, simImagT;
    protected double temperature;
    private double[] newU;
    private double[] wvCoeff, simWVCoeff, sqrtWVC;
    private double[][][] eigenVectors, simEigenVectors;
    private double[][] simOmegaSquared;

    protected final IRandom random;
    private IBox box;
    private int numAtoms;
    private Boundary bdry;
    private NormalModes nm;
    WaveVectorFactory waveVectorFactory;
    
    public MeterDifferentImageSubtract(String string, /*IPotentialMaster potentialMaster,*/ 
            int numSimAtoms, double density, Simulation sim,
            Primitive simPrimitive, Basis simBasis, CoordinateDefinition simCD,
            NormalModes simNM, double temp){
        super(string, Null.DIMENSION);
        this.random = sim.getRandom();
        this.temperature = temp;
        
        simWaveVectors = simNM.getWaveVectorFactory().getWaveVectors();
        this.simCDef = simCD;
        simCDim = simCD.getCoordinateDim();
        simEigenVectors = simNM.getEigenvectors();
        simWVCoeff = simNM.getWaveVectorFactory().getCoefficients();
        simRealT = new double[simCDim];
        simImagT = new double[simCDim];
        simOmegaSquared = simNM.getOmegaSquared();
        
        numAtoms = numSimAtoms - 1;
        box = new Box(sim.getSpace());
        sim.addBox(box);
        box.setNMolecules(sim.getSpecies(0), numAtoms);
        bdry = new BoundaryRectangularPeriodic(sim.getSpace(), numAtoms/density);
        box.setBoundary(bdry);
        
        int[] nCells = new int[]{numAtoms};
        cDef = new CoordinateDefinitionLeaf(box, simPrimitive, 
                simBasis, sim.getSpace());
        cDef.initializeCoordinates(nCells);
        cDim = cDef.getCoordinateDim();

        nm = new NormalModes1DHR(box.getBoundary(), numAtoms);
        nm.setHarmonicFudge(1.0);
        nm.setTemperature(temperature);
        nm.getOmegaSquared();
        waveVectorFactory = nm.getWaveVectorFactory();
        waveVectorFactory.makeWaveVectors(box);
        waveVectors = nm.getWaveVectorFactory().getWaveVectors();
        eigenVectors = nm.getEigenvectors();
        wvCoeff = nm.getWaveVectorFactory().getCoefficients();
        sqrtWVC = new double[wvCoeff.length];
        for (int i =0; i < wvCoeff.length; i++){
            sqrtWVC[i] = Math.sqrt(2*wvCoeff[i]);
        }
        
        PotentialMasterList potentialMaster = new PotentialMasterList(sim, sim.getSpace());
        Potential2 potential = new P2HardSphere(sim.getSpace(), 1.0, true);
        potential = new P2XOrder(sim.getSpace(), (Potential2HardSpherical)potential);
        potential.setBox(box);
        potentialMaster.addPotential(potential, new IAtomType[] {((SpeciesSpheresMono)sim.getSpecies(0)).getLeafType(), ((SpeciesSpheresMono)sim.getSpecies(0)).getLeafType()});
        double neighborRange = 1.01/density;
        potentialMaster.setRange(neighborRange);
        //find neighbors now.  Don't hook up NeighborListManager since the
        //  neighbors won't change
        potentialMaster.getNeighborManager(box).reset();
        
        meterPE = new MeterPotentialEnergy(potentialMaster);
        meterPE.setBox(box);
    }
    
    public double getDataAsScalar() {
        
        IAtomList atomlist = box.getLeafList();
//        for (int i = 0; i < atomlist.getAtomCount(); i++){
//            System.out.println("start i " + atomlist.getAtom(i).getPosition().getX(0));
//        }

        BasisCell[] simCells = simCDef.getBasisCells();
        BasisCell[] cells = cDef.getBasisCells();
        BasisCell cell = simCells[0];
        //nan this makes it 1D
        newU = new double[cDim];
        for(int i=0; i < newU.length; i++){
            newU[i] = 0.0;
        }
        
        //Calculate normal mode coordinates of simulation system.
        double[] realCoord = new double[simWaveVectors.length];
        double[] imagCoord = new double[simWaveVectors.length];
        for (int wvcount = 0; wvcount < simWaveVectors.length; wvcount++){
            simCDef.calcT(simWaveVectors[wvcount], simRealT, simImagT);
            realCoord[wvcount] = 0.0;
            imagCoord[wvcount] = 0.0;
            for (int i = 0; i < simCDim; i++){
                for (int j = 0; j < simCDim; j++){
                    realCoord[wvcount] += simEigenVectors[wvcount][i][j] * simRealT[j];
                    imagCoord[wvcount] += simEigenVectors[wvcount][i][j] * simImagT[j];
                }
            }
            //cleans up E-15 type numbers
            if(simWVCoeff[wvcount] != 1.0 ) {imagCoord[wvcount] = 0.0;}
            
            
        }
        
        //nan this will not work for more than 1D
        //Calculation of harmonic energy
        //The if statement determines whether the real coordinate or the 
        // imaginary coordinate is used, and zeroes out the last coordinate
        // calculated.
        double harmonic;
        int index = simWVCoeff.length-1;
        if(simWVCoeff[index] == 1.0){
            harmonic = simWVCoeff[index] * simOmegaSquared[index][0] * 
                    imagCoord[index] * imagCoord[index];
            imagCoord[index] = 0.0;
        } else {
            harmonic = simWVCoeff[index] * simOmegaSquared[index][0] * 
                    realCoord[index] * realCoord[index];
            realCoord[index] = 0.0;
        }
        
            
        //Calculate the positions for the meter's system
        for (int iCell = 0; iCell < cells.length; iCell++){
            cell = cells[iCell];
            for (int j = 0; j < cDim; j++) {
                newU[j] = 0.0;
            }
            for (int wvcount = 0; wvcount < waveVectors.length; wvcount++){
                //Calculate the change in positions.
                double kR = waveVectors[wvcount].dot(cell.cellPosition);
                double coskR = Math.cos(kR);
                double sinkR = Math.sin(kR);
                for (int i = 0; i < cDim; i++){
                    for (int j = 0; j < cDim; j++){
                       newU[j] += sqrtWVC[wvcount] * eigenVectors[wvcount][i][j] 
                            * (realCoord[wvcount] * coskR - imagCoord[wvcount] * sinkR);
                    }
                }
            }
            cDef.setToU(cells[iCell].molecules, newU);
        }
        
        return meterPE.getDataAsScalar() + harmonic;
    }

    
}<commit_msg>Added a commented-out comment for use in debugging.<commit_after>package etomica.models.oneDHardRods;

import etomica.api.IAtomList;
import etomica.api.IAtomType;
import etomica.api.IBox;
import etomica.api.IRandom;
import etomica.api.IVectorMutable;
import etomica.box.Box;
import etomica.data.DataSourceScalar;
import etomica.data.meter.MeterPotentialEnergy;
import etomica.lattice.crystal.Basis;
import etomica.lattice.crystal.Primitive;
import etomica.nbr.list.PotentialMasterList;
import etomica.normalmode.CoordinateDefinition;
import etomica.normalmode.CoordinateDefinitionLeaf;
import etomica.normalmode.NormalModes;
import etomica.normalmode.NormalModes1DHR;
import etomica.normalmode.P2XOrder;
import etomica.normalmode.WaveVectorFactory;
import etomica.normalmode.CoordinateDefinition.BasisCell;
import etomica.potential.P2HardSphere;
import etomica.potential.Potential2;
import etomica.potential.Potential2HardSpherical;
import etomica.simulation.Simulation;
import etomica.space.Boundary;
import etomica.space.BoundaryRectangularPeriodic;
import etomica.species.SpeciesSpheresMono;
import etomica.units.Null;


/**
 * Uses a Widom-like insertion of a mode to calculate a probability.
 * Uses a different box than the main simulation, to assume a mode & rod 
 * are removed
 * 
 * @author cribbin
 *
 */
public class MeterDifferentImageSubtract extends DataSourceScalar {

    public int nInsert, counter;
    private MeterPotentialEnergy meterPE;
    private CoordinateDefinition cDef, simCDef;
    private int cDim, simCDim;
    private IVectorMutable[] waveVectors, simWaveVectors;
    private double[] simRealT, simImagT;
    protected double temperature;
    private double[] newU;
    private double[] wvCoeff, simWVCoeff, sqrtWVC;
    private double[][][] eigenVectors, simEigenVectors;
    private double[][] simOmegaSquared;

    protected final IRandom random;
    private IBox box;
    private int numAtoms;
    private Boundary bdry;
    private NormalModes nm;
    WaveVectorFactory waveVectorFactory;
    
    public MeterDifferentImageSubtract(String string, /*IPotentialMaster potentialMaster,*/ 
            int numSimAtoms, double density, Simulation sim,
            Primitive simPrimitive, Basis simBasis, CoordinateDefinition simCD,
            NormalModes simNM, double temp){
        super(string, Null.DIMENSION);
        this.random = sim.getRandom();
        this.temperature = temp;
        
        simWaveVectors = simNM.getWaveVectorFactory().getWaveVectors();
        this.simCDef = simCD;
        simCDim = simCD.getCoordinateDim();
        simEigenVectors = simNM.getEigenvectors();
        simWVCoeff = simNM.getWaveVectorFactory().getCoefficients();
        simRealT = new double[simCDim];
        simImagT = new double[simCDim];
        simOmegaSquared = simNM.getOmegaSquared();
        
        numAtoms = numSimAtoms - 1;
        box = new Box(sim.getSpace());
        sim.addBox(box);
        box.setNMolecules(sim.getSpecies(0), numAtoms);
        bdry = new BoundaryRectangularPeriodic(sim.getSpace(), numAtoms/density);
        box.setBoundary(bdry);
        
        int[] nCells = new int[]{numAtoms};
        cDef = new CoordinateDefinitionLeaf(box, simPrimitive, 
                simBasis, sim.getSpace());
        cDef.initializeCoordinates(nCells);
        cDim = cDef.getCoordinateDim();

        nm = new NormalModes1DHR(box.getBoundary(), numAtoms);
        nm.setHarmonicFudge(1.0);
        nm.setTemperature(temperature);
        nm.getOmegaSquared();
        waveVectorFactory = nm.getWaveVectorFactory();
        waveVectorFactory.makeWaveVectors(box);
        waveVectors = nm.getWaveVectorFactory().getWaveVectors();
        eigenVectors = nm.getEigenvectors();
        wvCoeff = nm.getWaveVectorFactory().getCoefficients();
        sqrtWVC = new double[wvCoeff.length];
        for (int i =0; i < wvCoeff.length; i++){
            sqrtWVC[i] = Math.sqrt(2*wvCoeff[i]);
        }
        
        PotentialMasterList potentialMaster = new PotentialMasterList(sim, sim.getSpace());
        Potential2 potential = new P2HardSphere(sim.getSpace(), 1.0, true);
        potential = new P2XOrder(sim.getSpace(), (Potential2HardSpherical)potential);
        potential.setBox(box);
        potentialMaster.addPotential(potential, new IAtomType[] {((SpeciesSpheresMono)sim.getSpecies(0)).getLeafType(), ((SpeciesSpheresMono)sim.getSpecies(0)).getLeafType()});
        double neighborRange = 1.01/density;
        potentialMaster.setRange(neighborRange);
        //find neighbors now.  Don't hook up NeighborListManager since the
        //  neighbors won't change
        potentialMaster.getNeighborManager(box).reset();
        
        meterPE = new MeterPotentialEnergy(potentialMaster);
        meterPE.setBox(box);
    }
    
    public double getDataAsScalar() {
        
        IAtomList atomlist = box.getLeafList();
//        for (int i = 0; i < atomlist.getAtomCount(); i++){
//            System.out.println("start i " + atomlist.getAtom(i).getPosition().getX(0));
//        }

        BasisCell[] simCells = simCDef.getBasisCells();
        BasisCell[] cells = cDef.getBasisCells();
        BasisCell cell = simCells[0];
        //nan this makes it 1D
        newU = new double[cDim];
        for(int i=0; i < newU.length; i++){
            newU[i] = 0.0;
        }
        
        //Calculate normal mode coordinates of simulation system.
        double[] realCoord = new double[simWaveVectors.length];
        double[] imagCoord = new double[simWaveVectors.length];
        for (int wvcount = 0; wvcount < simWaveVectors.length; wvcount++){
            simCDef.calcT(simWaveVectors[wvcount], simRealT, simImagT);
            realCoord[wvcount] = 0.0;
            imagCoord[wvcount] = 0.0;
            for (int i = 0; i < simCDim; i++){
                for (int j = 0; j < simCDim; j++){
                    realCoord[wvcount] += simEigenVectors[wvcount][i][j] * simRealT[j];
                    imagCoord[wvcount] += simEigenVectors[wvcount][i][j] * simImagT[j];
                }
            }
            //cleans up E-15 type numbers
            if(simWVCoeff[wvcount] != 1.0 ) {imagCoord[wvcount] = 0.0;}
            
            
        }
        
        //nan this will not work for more than 1D
        //Calculation of harmonic energy
        //The if statement determines whether the real coordinate or the 
        // imaginary coordinate is used, and zeroes out the last coordinate
        // calculated.
        double harmonic;
        int index = simWVCoeff.length-1;
        if(simWVCoeff[index] == 1.0){
            harmonic = simWVCoeff[index] * simOmegaSquared[index][0] * 
                    imagCoord[index] * imagCoord[index];
            imagCoord[index] = 0.0;
        } else {
            harmonic = simWVCoeff[index] * simOmegaSquared[index][0] * 
                    realCoord[index] * realCoord[index];
            realCoord[index] = 0.0;
        }
        
            
        //Calculate the positions for the meter's system
        for (int iCell = 0; iCell < cells.length; iCell++){
            cell = cells[iCell];
            for (int j = 0; j < cDim; j++) {
                newU[j] = 0.0;
            }
            for (int wvcount = 0; wvcount < waveVectors.length; wvcount++){
                //Calculate the change in positions.
                double kR = waveVectors[wvcount].dot(cell.cellPosition);
                double coskR = Math.cos(kR);
                double sinkR = Math.sin(kR);
                for (int i = 0; i < cDim; i++){
                    for (int j = 0; j < cDim; j++){
                       newU[j] += sqrtWVC[wvcount] * eigenVectors[wvcount][i][j] 
                            * (realCoord[wvcount] * coskR - imagCoord[wvcount] * sinkR);
                    }
                }
            }
            cDef.setToU(cells[iCell].molecules, newU);
        }
        
//        return Math.exp(-1*(meterPE.getDataAsScalar()+harmonic));
        return meterPE.getDataAsScalar() + harmonic;
    }

    
}
{'fraction_non_alphanumeric': 0.10613266583229036, 'fraction_numerical': 0.007571964956195244, 'mean_word_length': 2.4885396201702683, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5529390', 'n_tokens_mistral': 5191, 'n_tokens_neox': 4899, 'n_words': 1037}