<commit_before>package etomica.models.oneDHardRods;

import etomica.api.IAtomList;
import etomica.api.IBox;
import etomica.api.IPotentialMaster;
import etomica.api.IRandom;
import etomica.api.IVectorMutable;
import etomica.atom.Atom;
import etomica.atom.iterator.AtomIterator;
import etomica.atom.iterator.AtomIteratorLeafAtoms;
import etomica.data.meter.MeterPotentialEnergy;
import etomica.integrator.mcmove.MCMoveBoxStep;
import etomica.normalmode.CoordinateDefinition;
import etomica.normalmode.CoordinateDefinition.BasisCell;

/**
 * A Monte Carlo move which compares several normal modes to harmonic normal
 * modes and then explores the phase space defined by the remaining normal
 * modes.
 * 
 * harmonicWVs are the wavevectors that are harmonic, and left out of the "change
 * a random mode" calculation of the doTrial() method.
 * 
 * comparedWVs are the wavevectors that are compared (removed).
 * 
 * @author cribbin
 * 
 */
public class MCMoveCompareMultipleWVLoop extends MCMoveBoxStep {

    private static final long serialVersionUID = 1L;
    protected CoordinateDefinition coordinateDefinition;
    private final AtomIteratorLeafAtoms iterator;
    protected double[][] uOld;
    protected double[] deltaU;
    protected final IRandom random;
    protected double energyOld,energyNew, energyEvenLater;
    protected final MeterPotentialEnergy energyMeter;
    private double[][][] eigenVectors;
    private IVectorMutable[] waveVectors;
    private double[] gaussian;
    protected double temperature;
    private double[][] stdDev;
    private double[] rRand, iRand, realT, imagT;
    private double[] waveVectorCoefficients;
    private double[][] omega2;
    double[] uNow;
    int changedWV, howManyChangesToHardRodModes;
    int[] comparedWVs, harmonicWVs;
    
    public MCMoveCompareMultipleWVLoop(IPotentialMaster potentialMaster,
            IRandom random) {
        super(potentialMaster);
        this.random = random;
        iterator = new AtomIteratorLeafAtoms();
        energyMeter = new MeterPotentialEnergy(potentialMaster);
        gaussian = new double[2];
        howManyChangesToHardRodModes = 1;
    }
    
    public boolean doTrial() {
        int coordinateDim = coordinateDefinition.getCoordinateDim();
        BasisCell[] cells = coordinateDefinition.getBasisCells();
        rRand = new double[coordinateDim];
        iRand = new double[coordinateDim];
        realT = new double[coordinateDim];
        imagT = new double[coordinateDim];
        
        BasisCell cell = cells[0];
        uOld = new double[cells.length][coordinateDim];
        double normalization = 1 / Math.sqrt(cells.length);
        int numWV = comparedWVs.length;
        
        /*
         * This loop looks at each wavevector, asks if that wavevector is
         * removed, and calculates what happens if it is.
         */
// ZERO OUT A NORMAL MODE.
        
        //Store old positions
        for(int iCell = 0; iCell < cells.length; iCell++){
            //store old positions.
            uNow = coordinateDefinition.calcU(cells[iCell].molecules);
            System.arraycopy(uNow, 0, uOld[iCell], 0, coordinateDim);
        }
        for (int wvCount = 0; wvCount < numWV; wvCount++) {
            int comparedwv = comparedWVs[wvCount];
            
            // Get normal mode coordinate information
            coordinateDefinition.calcT(waveVectors[comparedwv], realT, imagT);
            
            for (int iCell = 0; iCell < cells.length; iCell++) {
                cell = cells[iCell];
                uNow = coordinateDefinition.calcU(cells[iCell].molecules);
                // rezero deltaU
                for (int j = 0; j < coordinateDim; j++) {
                    deltaU[j] = 0.0;
                }

                // Calculate the contributions to the current position of
                // the zeroed mode, and subtract it from the overall position
                double kR = waveVectors[comparedwv].dot(cell.cellPosition);
                double coskR = Math.cos(kR);
                double sinkR = Math.sin(kR);
                for (int i = 0; i < coordinateDim; i++) {
                    // Calculate the current coordinate:
                    double realCoord = 0, imagCoord = 0;
                    for (int j = 0; j < coordinateDim; j++) {
                        realCoord += eigenVectors[comparedwv][i][j] * realT[j];
                        imagCoord += eigenVectors[comparedwv][i][j] * imagT[j];
                    }
                    for (int j = 0; j < coordinateDim; j++) {
                        deltaU[j] -= waveVectorCoefficients[comparedwv] 
                                * eigenVectors[comparedwv][i][j] * 2.0
                                * (realCoord * coskR - imagCoord * sinkR);
                    }
                }
                for (int i = 0; i < coordinateDim; i++) {
                    deltaU[i] *= normalization;
                }

                for (int i = 0; i < coordinateDim; i++) {
                    uNow[i] += deltaU[i];
                }
                coordinateDefinition.setToU(cells[iCell].molecules, uNow);
            }//end of cell loop
        }// end of wvCount loop
        
        energyOld = energyMeter.getDataAsScalar();
        if (Double.isInfinite(energyOld)) {
            int limit = coordinateDefinition.getBox().getLeafList().getAtomCount();
            for (int k = 0; k < limit; k++) {
                System.out.println(k + " " + 
                        coordinateDefinition.getBox().getLeafList().getAtom(k)
                        .getPosition());
            }
            throw new IllegalStateException(
                    "Overlap after the removal of a mode!");
        }

// MOVE SOME NUMBER OF RANDOM HARD ROD POTENTIAL MODES, AND MEASURE energyNew
        // equivalent to MCMoveChangeMode for several modes.
//        int changedWV;
        for(int wvCount = 0; wvCount < howManyChangesToHardRodModes; wvCount++){
            // Select the wave vector whose eigenvectors will be changed.
            // The zero wavevector is center of mass motion, and is rejected as
            // a possibility.
            boolean isAccepted = true;
            do{
                isAccepted = true;
                changedWV = random.nextInt(waveVectorCoefficients.length);
                for(int i = 0; i < harmonicWVs.length; i++){
                    if (changedWV == harmonicWVs[i]) {
                        isAccepted = false;
                    }
                }
            } while (!isAccepted);
//            System.out.println(changedWV);
            
            // calculate the new positions of the atoms.
            // loop over cells
            double[] delta = new double[coordinateDim*2];<commit_msg>minor cleanup of useless commented out line<commit_after>package etomica.models.oneDHardRods;

import etomica.api.IAtomList;
import etomica.api.IBox;
import etomica.api.IPotentialMaster;
import etomica.api.IRandom;
import etomica.api.IVectorMutable;
import etomica.atom.Atom;
import etomica.atom.iterator.AtomIterator;
import etomica.atom.iterator.AtomIteratorLeafAtoms;
import etomica.data.meter.MeterPotentialEnergy;
import etomica.integrator.mcmove.MCMoveBoxStep;
import etomica.normalmode.CoordinateDefinition;
import etomica.normalmode.CoordinateDefinition.BasisCell;

/**
 * A Monte Carlo move which compares several normal modes to harmonic normal
 * modes and then explores the phase space defined by the remaining normal
 * modes.
 * 
 * harmonicWVs are the wavevectors that are harmonic, and left out of the "change
 * a random mode" calculation of the doTrial() method.
 * 
 * comparedWVs are the wavevectors that are compared (removed).
 * 
 * @author cribbin
 * 
 */
public class MCMoveCompareMultipleWVLoop extends MCMoveBoxStep {

    private static final long serialVersionUID = 1L;
    protected CoordinateDefinition coordinateDefinition;
    private final AtomIteratorLeafAtoms iterator;
    protected double[][] uOld;
    protected double[] deltaU;
    protected final IRandom random;
    protected double energyOld,energyNew, energyEvenLater;
    protected final MeterPotentialEnergy energyMeter;
    private double[][][] eigenVectors;
    private IVectorMutable[] waveVectors;
    private double[] gaussian;
    protected double temperature;
    private double[][] stdDev;
    private double[] rRand, iRand, realT, imagT;
    private double[] waveVectorCoefficients;
    private double[][] omega2;
    double[] uNow;
    int changedWV, howManyChangesToHardRodModes;
    int[] comparedWVs, harmonicWVs;
    
    public MCMoveCompareMultipleWVLoop(IPotentialMaster potentialMaster,
            IRandom random) {
        super(potentialMaster);
        this.random = random;
        iterator = new AtomIteratorLeafAtoms();
        energyMeter = new MeterPotentialEnergy(potentialMaster);
        gaussian = new double[2];
        howManyChangesToHardRodModes = 1;
    }
    
    public boolean doTrial() {
        int coordinateDim = coordinateDefinition.getCoordinateDim();
        BasisCell[] cells = coordinateDefinition.getBasisCells();
        rRand = new double[coordinateDim];
        iRand = new double[coordinateDim];
        realT = new double[coordinateDim];
        imagT = new double[coordinateDim];
        
        BasisCell cell = cells[0];
        uOld = new double[cells.length][coordinateDim];
        double normalization = 1 / Math.sqrt(cells.length);
        int numWV = comparedWVs.length;
        
        /*
         * This loop looks at each wavevector, asks if that wavevector is
         * removed, and calculates what happens if it is.
         */
// ZERO OUT A NORMAL MODE.
        
        //Store old positions
        for(int iCell = 0; iCell < cells.length; iCell++){
            //store old positions.
            uNow = coordinateDefinition.calcU(cells[iCell].molecules);
            System.arraycopy(uNow, 0, uOld[iCell], 0, coordinateDim);
        }
        for (int wvCount = 0; wvCount < numWV; wvCount++) {
            int comparedwv = comparedWVs[wvCount];
            
            // Get normal mode coordinate information
            coordinateDefinition.calcT(waveVectors[comparedwv], realT, imagT);
            
            for (int iCell = 0; iCell < cells.length; iCell++) {
                cell = cells[iCell];
                uNow = coordinateDefinition.calcU(cells[iCell].molecules);
                // rezero deltaU
                for (int j = 0; j < coordinateDim; j++) {
                    deltaU[j] = 0.0;
                }

                // Calculate the contributions to the current position of
                // the zeroed mode, and subtract it from the overall position
                double kR = waveVectors[comparedwv].dot(cell.cellPosition);
                double coskR = Math.cos(kR);
                double sinkR = Math.sin(kR);
                for (int i = 0; i < coordinateDim; i++) {
                    // Calculate the current coordinate:
                    double realCoord = 0, imagCoord = 0;
                    for (int j = 0; j < coordinateDim; j++) {
                        realCoord += eigenVectors[comparedwv][i][j] * realT[j];
                        imagCoord += eigenVectors[comparedwv][i][j] * imagT[j];
                    }
                    for (int j = 0; j < coordinateDim; j++) {
                        deltaU[j] -= waveVectorCoefficients[comparedwv] 
                                * eigenVectors[comparedwv][i][j] * 2.0
                                * (realCoord * coskR - imagCoord * sinkR);
                    }
                }
                for (int i = 0; i < coordinateDim; i++) {
                    deltaU[i] *= normalization;
                }

                for (int i = 0; i < coordinateDim; i++) {
                    uNow[i] += deltaU[i];
                }
                coordinateDefinition.setToU(cells[iCell].molecules, uNow);
            }//end of cell loop
        }// end of wvCount loop
        
        energyOld = energyMeter.getDataAsScalar();
        if (Double.isInfinite(energyOld)) {
            int limit = coordinateDefinition.getBox().getLeafList().getAtomCount();
            for (int k = 0; k < limit; k++) {
                System.out.println(k + " " + 
                        coordinateDefinition.getBox().getLeafList().getAtom(k)
                        .getPosition());
            }
            throw new IllegalStateException(
                    "Overlap after the removal of a mode!");
        }

// MOVE SOME NUMBER OF RANDOM HARD ROD POTENTIAL MODES, AND MEASURE energyNew
        // equivalent to MCMoveChangeMode for several modes.
        for(int wvCount = 0; wvCount < howManyChangesToHardRodModes; wvCount++){
            // Select the wave vector whose eigenvectors will be changed.
            // The zero wavevector is center of mass motion, and is rejected as
            // a possibility.
            boolean isAccepted = true;
            do{
                isAccepted = true;
                changedWV = random.nextInt(waveVectorCoefficients.length);
                for(int i = 0; i < harmonicWVs.length; i++){
                    if (changedWV == harmonicWVs[i]) {
                        isAccepted = false;
                    }
                }
            } while (!isAccepted);
//            System.out.println(changedWV);
            
            // calculate the new positions of the atoms.
            // loop over cells
            double[] delta = new double[coordinateDim*2];
{'fraction_non_alphanumeric': 0.09043927648578812, 'fraction_numerical': 0.003986710963455149, 'mean_word_length': 1.9712656284272867, 'pattern_counts': {'":': 0, '<': 27, '<?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': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3404765', 'n_tokens_mistral': 3868, 'n_tokens_neox': 3579, 'n_words': 1003}