 package etomica.normalmode;

import etomica.lattice.crystal.Basis;
import etomica.lattice.crystal.BasisCubicFcc;
import etomica.lattice.crystal.Primitive;
import etomica.lattice.crystal.PrimitiveCubic;
import etomica.potential.P2SoftSphere;
import etomica.potential.P2SoftSphericalTruncated;
import etomica.potential.Potential2SoftSpherical;
import etomica.space.ISpace;
import etomica.space.Space;
import etomica.space3d.Space3D;

/**
 * Properties of a system of monatomic molecules occupying a lattice and interacting according
 * to a soft-sphere spherically-symmetric pair potential.  Properties are given by a lattice-dynamics treatment.
 * 
 * FCC Crystal Structure
 * 
 * @author kofke & Tai Tan
 *
 */
public class HarmonicCrystalSoftSphereFCC {

    public static void main(String[] args) {
        double rho = 1.256;
        int exponent = 12;
        int maxLatticeShell = 49;
        int nC =5;
        int numAtom = nC*nC*nC*4;
        
        double temperature = 0.001;
        String fileName = "inputSSDB"+ numAtom;
//        Primitive primitive = new PrimitiveFcc(Space3D.getInstance());
//        Basis basis = new BasisMonatomic(Space3D.getInstance());
        
        if (args.length > 0) {
            rho = Double.parseDouble(args[0]);
        }
        if (args.length > 1) {
            exponent = Integer.parseInt(args[1]);
        }
        if (args.length > 2) {
            nC = Integer.parseInt(args[2]);
        }
        if (args.length > 3) {
            temperature = Double.parseDouble(args[3]);
        }
        if (args.length > 4) {
        	fileName = args[4];
        }
        
        Primitive primitive = new PrimitiveCubic(Space3D.getInstance());
        Basis basisFCC = new BasisCubicFcc();
        
        Basis basis = new BasisBigCell(Space.getInstance(3), basisFCC, new int[]{nC,nC,nC});
        
        ISpace sp = Space3D.getInstance();
        final Potential2SoftSpherical potential = new P2SoftSphere(sp, 1.0, 1.0, exponent);
        double truncationRadius = 3.64137;
        P2SoftSphericalTruncated pTruncated = new P2SoftSphericalTruncated(sp, potential, truncationRadius);
        
        int[] nCells = new int[] {1, 1, 1};
        long startTime = System.currentTimeMillis();
        System.out.println("Start Time: " + startTime);
        
        System.out.println("Running lattice-dynamics derivatives-based FCC soft-sphere simulation");
        System.out.println("Temperature: " + temperature);
        System.out.println("Density: " + rho);
        System.out.println("Exponent: " + exponent +"\n");
        
        
        
        HarmonicCrystal harmonicCrystal = new HarmonicCrystal(nCells, primitive, basis, pTruncated, sp);
       
        harmonicCrystal.setCellDensity(rho/basis.getScaledCoordinates().length);

        double u = harmonicCrystal.getLatticeEnergy();
        harmonicCrystal.getNormalModes().setFileName(fileName);
        double a = harmonicCrystal.getHelmholtzFreeEnergy(temperature);
        
        System.out.println("\nLattice Energy: " + u);
        System.out.println("Helmholtz Free Energy at T"+temperature+ " is: "+a);
        System.out.println("Harmonic-reference free energy: "+ (a-u));
      
        System.out.println("\nCalcHarmonicA from file (Temperature-independent)");
        CalcHarmonicA.doit(harmonicCrystal.getNormalModes(), 3, temperature, basis.getScaledCoordinates().length);

        long endTime = System.currentTimeMillis();
        System.out.println("End Time: " + endTime);
        System.out.println("Time taken: " + (endTime - startTime));
        

    }
}

@@ -884,10 +884,9 @@
 l = 
-49
+2
 ;%0A  
@@ -899,16 +899,17 @@
 int nC =
+ 
 5;%0A     
@@ -2002,34 +2002,117 @@
 ble 
-truncationRadius = 3.64137
+rc = (maxLatticeShell*2)*nC*Math.pow(4.0/rho, 1.0/3.0)*0.495;%0A        System.out.println(%22truncation at %22+rc)
 ;%0A  
@@ -2203,24 +2203,10 @@
 al, 
-truncationRadius
+rc
 );%0A 
@@ -2662,41 +2662,588 @@
     
-%0A        %0A        HarmonicCrystal
+// we want an untruncated lattice energy, so extend lattice sum by 2 more cells%0A        double rcU = ((maxLatticeShell+2)*2)*nC*Math.pow(4.0/rho, 1.0/3.0)*0.495;%0A        pTruncated.setTruncationRadius(rcU);%0A        HarmonicCrystal harmonicCrystal = new HarmonicCrystal(nCells, primitive, basis, pTruncated, sp);%0A        %0A        harmonicCrystal.setCellDensity(rho/basis.getScaledCoordinates().length);%0A%0A        harmonicCrystal.setMaxLatticeShell(maxLatticeShell+2);%0A        double u = harmonicCrystal.getLatticeEnergy();%0A%0A%0A        // now actually calculate the free energy%0A       
  har
@@ -3413,35 +3413,24 @@
 h);%0A%0A       
- double u =
  harmonicCry
@@ -3434,33 +3434,50 @@
 Crystal.
-g
+s
 et
+Max
 Lattice
-Energy(
+Shell(maxLatticeShell
 );%0A     

use a more coherent truncation don't do maxLatticeShell=49 (crazy!) try to get untruncated lattice energy