<commit_before>package etomica.models.oneDHardRods;

import etomica.data.DataInfo;
import etomica.data.IData;
import etomica.data.IDataInfo;
import etomica.data.types.DataDouble;
import etomica.lattice.BravaisLatticeCrystal;
import etomica.lattice.LatticeSumCrystal;
import etomica.lattice.LatticeSumCrystal.DataGroupLSC;
import etomica.lattice.crystal.Basis;
import etomica.lattice.crystal.BasisCubicFcc;
import etomica.lattice.crystal.Primitive;
import etomica.lattice.crystal.PrimitiveCubic;
import etomica.normalmode.CalcHarmonicA;
import etomica.normalmode.NormalModesPotential;
import etomica.potential.P2SoftSphere;
import etomica.potential.P2SoftSphericalTruncated;
import etomica.potential.Potential2SoftSpherical;
import etomica.space.ISpace;
import etomica.space.Space;
import etomica.space3d.Space3D;
import etomica.space3d.Vector3D;
import etomica.units.Energy;
import etomica.util.FunctionGeneral;
import etomica.util.ParameterBase;

/**
 * Properties of a system of monatomic molecules occupying a lattice and 
 * interacting according to a spherically-symmetric pair potential.  Properties 
 * are given by a lattice-dynamics treatment.
 * 
 * set up to use truncated soft sphere potential and output files 
 * 
 * @author cribbin off HarmonicCrystal of kofke
 * 
 */
public class HarmonicCrystalSsFccNxy {

    public HarmonicCrystalSsFccNxy(int[] nCells, Primitive primitive, Basis 
            basis, Potential2SoftSpherical potential, ISpace _space) {
        this.potential = potential;
        this.nCells = (int[])nCells.clone();
        this.space = _space;
        lattice = new BravaisLatticeCrystal(primitive, basis);
        normalModes = new NormalModesPotential(nCells, primitive, basis, 
                potential, space);
    }

    public NormalModesPotential getNormalModes() {
        return normalModes;
    }

    public double getLatticeEnergy() {
        FunctionGeneral function = new FunctionGeneral() {
            public IData f(Object obj) {
                data.x = potential.u(((Vector3D)obj).squared());
                return data;
            }
            public IDataInfo getDataInfo() {
                return dataInfo;
            }
            final DataInfo dataInfo = new DataDouble.DataInfoDouble("Lattice " +
            		"energy", Energy.DIMENSION);
            final DataDouble data = new DataDouble();
        };
        LatticeSumCrystal summer = new LatticeSumCrystal(lattice);
        summer.setMaxLatticeShell(maxLatticeShell);
        summer.setK(lattice.getSpace().makeVector());
        double sum = 0;
        double basisDim = lattice.getBasis().getScaledCoordinates().length;
        DataGroupLSC data = (DataGroupLSC)summer.calculateSum(function);
        for(int j=0; j<basisDim; j++) {
            for(int jp=0; jp<basisDim; jp++) {
                sum += ((DataDouble)data.getDataReal(j,jp)).x; 
            }
        }
        return 0.5*sum/basisDim;
    }
    
    public double getHelmholtzFreeEnergy(double temperature) {

        int D = lattice.getSpace().D();
        int cellCount = 1;
        int differ = 1;
        for(int i=0; i<D; i++) {
            cellCount *= nCells[i];
            if(nCells[i] % 2 == 0) {
                differ *= 2;
            }
        }
        int basisDim = lattice.getBasis().getScaledCoordinates().length;
        int moleculeCount = cellCount*basisDim;
        double Acom = 0.5*D*Math.log(moleculeCount);

        double[][] omega2 = normalModes.getOmegaSquared();
        double[] coeffs = normalModes.getWaveVectorFactory().getCoefficients();
        double sumA = 0.0;
        double normalModeSum = 0.0;
        double omega2zeroCount = 0;
        for(int k=0; k<omega2.length; k++) {
            double coeff = coeffs[k];
            for(int i=0; i<omega2[k].length; i++) {
                if(!Double.isInfinite(omega2[k][i])) {
                    sumA += coeff*Math.log(omega2[k][i]/(2*temperature*Math.PI));
                    normalModeSum += coeff;
                } else {
                    omega2zeroCount++;
                }
            }
        }

        sumA -= Acom;
        sumA /= moleculeCount;
        sumA *= temperature;
        return sumA;
    }
    
    public void setCellDensity(double newDensity) {
        double oldVolume = lattice.getPrimitive().unitCell().getVolume();
        double scale = newDensity * oldVolume;
        Primitive primitive = lattice.getPrimitive();
        primitive.scaleSize(1.0/Math.pow(scale, 1.0/lattice.getSpace().D()));
    }
    
    public int getMaxLatticeShell() {
        return maxLatticeShell;
    }

    public void setMaxLatticeShell(int maxLatticeShell) {
        this.maxLatticeShell = maxLatticeShell;
        normalModes.setMaxLatticeShell(maxLatticeShell);
    }

    public static void main(String[] args) {
        Params params = new Params();
        
        double T = params.T;
        double rho = params.rho;
        int[] nCells = params.shape;
        int nA = (nCells[0] * nCells[1] * nCells[2] * 4);
        String filename = params.filename + nA;
        double rc = params.rc ;
        
        System.out.println("truncation " + rc);
        
        Space sp = Space3D.getInstance();
        Potential2SoftSpherical potentialBase = new P2SoftSphere(sp, 1.0,
                1.0, 12);
        P2SoftSphericalTruncated potential = new P2SoftSphericalTruncated(
                sp, potentialBase, rc);
        potential.setTruncationRadius(rc);
        
        Primitive primitive = new PrimitiveCubic(sp);
        Basis basis = new BasisCubicFcc();
        
        double L = Math.pow(basis.getScaledCoordinates().length / rho, 1.0 / 3.0);
        primitive.scaleSize(L);
        
        HarmonicCrystalSsFccNxy harmonicCrystal = new HarmonicCrystalSsFccNxy(
                nCells, primitive, basis, potential, sp);
        harmonicCrystal.setMaxLatticeShell(3);
        harmonicCrystal.getNormalModes().setFileName(filename);
        
        System.out.println("Density: " + rho);
        System.out.println("Temperature: " + T);
        double a = harmonicCrystal.getHelmholtzFreeEnergy(T);
        double u = harmonicCrystal.getLatticeEnergy();
                
        System.out.println("\nLattice Energy: " + u);
        System.out.println("Helmholtz Free Energy at T "+T+ " is: "+a);
        System.out.println("Harmonic-reference free energy: "+ (a-u));
      
        System.out.println("\nCalcHarmonicA from file (Temperature-independent)");
        CalcHarmonicA.doit(harmonicCrystal.getNormalModes(), 3, T, nA);

    }
    
    private NormalModesPotential normalModes;
    private BravaisLatticeCrystal lattice;
    private int[] nCells;
    private int maxLatticeShell;
    private Potential2SoftSpherical potential;
    private final ISpace space;
    private static final long serialVersionUID = 1L;
    
    
    public static class Params extends ParameterBase {
        public double T = 0.01;
        public double rho = 1.1964;
        public int[] shape = new int[] {2, 2, 4};
        public String filename = "inputSSDB_WV";
        public double rc = 1.4803453945760225 ;
    }
    
}<commit_msg>both truncation radii now listed<commit_after>package etomica.models.oneDHardRods;

import etomica.data.DataInfo;
import etomica.data.IData;
import etomica.data.IDataInfo;
import etomica.data.types.DataDouble;
import etomica.lattice.BravaisLatticeCrystal;
import etomica.lattice.LatticeSumCrystal;
import etomica.lattice.LatticeSumCrystal.DataGroupLSC;
import etomica.lattice.crystal.Basis;
import etomica.lattice.crystal.BasisCubicFcc;
import etomica.lattice.crystal.Primitive;
import etomica.lattice.crystal.PrimitiveCubic;
import etomica.normalmode.CalcHarmonicA;
import etomica.normalmode.NormalModesPotential;
import etomica.potential.P2SoftSphere;
import etomica.potential.P2SoftSphericalTruncated;
import etomica.potential.Potential2SoftSpherical;
import etomica.space.ISpace;
import etomica.space.Space;
import etomica.space3d.Space3D;
import etomica.space3d.Vector3D;
import etomica.units.Energy;
import etomica.util.FunctionGeneral;
import etomica.util.ParameterBase;

/**
 * Properties of a system of monatomic molecules occupying a lattice and 
 * interacting according to a spherically-symmetric pair potential.  Properties 
 * are given by a lattice-dynamics treatment.
 * 
 * set up to use truncated soft sphere potential and output files 
 * 
 * @author cribbin off HarmonicCrystal of kofke
 * 
 */
public class HarmonicCrystalSsFccNxy {

    public HarmonicCrystalSsFccNxy(int[] nCells, Primitive primitive, Basis 
            basis, Potential2SoftSpherical potential, ISpace _space) {
        this.potential = potential;
        this.nCells = (int[])nCells.clone();
        this.space = _space;
        lattice = new BravaisLatticeCrystal(primitive, basis);
        normalModes = new NormalModesPotential(nCells, primitive, basis, 
                potential, space);
    }

    public NormalModesPotential getNormalModes() {
        return normalModes;
    }

    public double getLatticeEnergy() {
        FunctionGeneral function = new FunctionGeneral() {
            public IData f(Object obj) {
                data.x = potential.u(((Vector3D)obj).squared());
                return data;
            }
            public IDataInfo getDataInfo() {
                return dataInfo;
            }
            final DataInfo dataInfo = new DataDouble.DataInfoDouble("Lattice " +
            		"energy", Energy.DIMENSION);
            final DataDouble data = new DataDouble();
        };
        LatticeSumCrystal summer = new LatticeSumCrystal(lattice);
        summer.setMaxLatticeShell(maxLatticeShell);
        summer.setK(lattice.getSpace().makeVector());
        double sum = 0;
        double basisDim = lattice.getBasis().getScaledCoordinates().length;
        DataGroupLSC data = (DataGroupLSC)summer.calculateSum(function);
        for(int j=0; j<basisDim; j++) {
            for(int jp=0; jp<basisDim; jp++) {
                sum += ((DataDouble)data.getDataReal(j,jp)).x; 
            }
        }
        return 0.5*sum/basisDim;
    }
    
    public double getHelmholtzFreeEnergy(double temperature) {

        int D = lattice.getSpace().D();
        int cellCount = 1;
        int differ = 1;
        for(int i=0; i<D; i++) {
            cellCount *= nCells[i];
            if(nCells[i] % 2 == 0) {
                differ *= 2;
            }
        }
        int basisDim = lattice.getBasis().getScaledCoordinates().length;
        int moleculeCount = cellCount*basisDim;
        double Acom = 0.5*D*Math.log(moleculeCount);

        double[][] omega2 = normalModes.getOmegaSquared();
        double[] coeffs = normalModes.getWaveVectorFactory().getCoefficients();
        double sumA = 0.0;
        double normalModeSum = 0.0;
        double omega2zeroCount = 0;
        for(int k=0; k<omega2.length; k++) {
            double coeff = coeffs[k];
            for(int i=0; i<omega2[k].length; i++) {
                if(!Double.isInfinite(omega2[k][i])) {
                    sumA += coeff*Math.log(omega2[k][i]/(2*temperature*Math.PI));
                    normalModeSum += coeff;
                } else {
                    omega2zeroCount++;
                }
            }
        }

        sumA -= Acom;
        sumA /= moleculeCount;
        sumA *= temperature;
        return sumA;
    }
    
    public void setCellDensity(double newDensity) {
        double oldVolume = lattice.getPrimitive().unitCell().getVolume();
        double scale = newDensity * oldVolume;
        Primitive primitive = lattice.getPrimitive();
        primitive.scaleSize(1.0/Math.pow(scale, 1.0/lattice.getSpace().D()));
    }
    
    public int getMaxLatticeShell() {
        return maxLatticeShell;
    }

    public void setMaxLatticeShell(int maxLatticeShell) {
        this.maxLatticeShell = maxLatticeShell;
        normalModes.setMaxLatticeShell(maxLatticeShell);
    }

    public static void main(String[] args) {
        Params params = new Params();
        
        double T = params.T;
        double rho = params.rho;
        int[] nCells = params.shape;
        int nA = (nCells[0] * nCells[1] * nCells[2] * 4);
        String filename = params.filename + nA;
        double rc = params.rc ;
        
        System.out.println("truncation " + rc);
        
        Space sp = Space3D.getInstance();
        Potential2SoftSpherical potentialBase = new P2SoftSphere(sp, 1.0,
                1.0, 12);
        P2SoftSphericalTruncated potential = new P2SoftSphericalTruncated(
                sp, potentialBase, rc);
        potential.setTruncationRadius(rc);
        
        Primitive primitive = new PrimitiveCubic(sp);
        Basis basis = new BasisCubicFcc();
        
        double L = Math.pow(basis.getScaledCoordinates().length / rho, 1.0 / 3.0);
        primitive.scaleSize(L);
        
        HarmonicCrystalSsFccNxy harmonicCrystal = new HarmonicCrystalSsFccNxy(
                nCells, primitive, basis, potential, sp);
        harmonicCrystal.setMaxLatticeShell(3);
        harmonicCrystal.getNormalModes().setFileName(filename);
        
        System.out.println("Density: " + rho);
        System.out.println("Temperature: " + T);
        double a = harmonicCrystal.getHelmholtzFreeEnergy(T);
        double u = harmonicCrystal.getLatticeEnergy();
                
        System.out.println("\nLattice Energy: " + u);
        System.out.println("Helmholtz Free Energy at T "+T+ " is: "+a);
        System.out.println("Harmonic-reference free energy: "+ (a-u));
      
        System.out.println("\nCalcHarmonicA from file (Temperature-independent)");
        CalcHarmonicA.doit(harmonicCrystal.getNormalModes(), 3, T, nA);

    }
    
    private NormalModesPotential normalModes;
    private BravaisLatticeCrystal lattice;
    private int[] nCells;
    private int maxLatticeShell;
    private Potential2SoftSpherical potential;
    private final ISpace space;
    private static final long serialVersionUID = 1L;
    
    
    public static class Params extends ParameterBase {
        public double T = 0.1;
        public double rho = 1.1964;
        public int[] shape = new int[] {16, 16, 16};
        public String filename = "inputSSDB_WV";
       // public double rc = 1.4803453945760225 ;
        public double rc = 2.2;
    }
    
}
{'fraction_non_alphanumeric': 0.10159682030541803, 'fraction_numerical': 0.013109267136182972, 'mean_word_length': 2.476848484848485, 'pattern_counts': {'":': 0, '<': 13, '<?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': '6988388', 'n_tokens_mistral': 4601, 'n_tokens_neox': 4296, 'n_words': 912}